Adding new potential energy surfaces

Adding new features to dynamiq always involves two main parts: adding in the code, and adding support for the new feature in the input handling.

A new potential generally requires one pair of files in the potentials/ directory. Parameters are generally implemented as gsl_vectors or gsl_matrices in the .h file.

You'll have to implement several functions. First, there are two functions which take a pointer to an md_sys and return a double. These are for the potential and kinetic energies. While the kinetic energy is typically the standard kinetic energy, the real purpose of it in the calculations is that the difference of these terms should be the Lagrangian of the system. This is particularly important when writing up potentials that use a "classical electronic state" model of nonadiabatic dynamics (i.e., the MMST Hamiltonian).

There should also be two functions which take a pointer to an md_sys and return a pointer to a gsl_vector. These are the qdot and pdot functions, given by Hamilton's equations.

There are also two (gsl_)matrices which need to be returned (also from functions which take a pointer to md_sys). The first is the Hessian, the second is the inverse mass matrix (formally, the second derivative of the Hamiltonian with respect to momentum -- for cases of nontrivial kinetic energy expressions). If needed (again, MMST Hamiltonian being an example) one can also define functions to return the mixed derivatives of the Hamiltonian with respect to position and momentum.

There are some cases where the sampling is done in a different coordinate system than the dynamics. If this is your case, you should also write an "mc2dyn" function to convert before each trajectory.

All of these functions are made visible to the program through an initialization function, which is a void that takes pointers to an md_sys and an md_pot, and the integer number of dimensions. The md_pot input should be initialized so that all the functions described above are assigned to the correct elements of that structure.

Once the code for the potential energy surface is done, you need to write the support for the input. Start with the file lexing/md_lexer.l. Add a regex that returns a token for your potential in POTENTIALS section of that file.

In lexing/md_parser.y, add your token to the "potentials" tokens. Under the POTENTIALS part of the objtok heading, push your token onto the stack, and set thisCalc->pot->init to your initialization function. Then go to the "PARAMETERS FOR POTENTIALS" section. You'll want to start off with the line printf0("\n* "); for output purposes. Then every potential must parse the keywords "DIM" and "M" into the integer sys->ndim and the gsl_vector sys->m (respectively). This is also where you should parse any system-specific parameters you'd like to add.

Recompile (using make clean && make) and you should be ready to go!