dynamic is able to perform molecular dynamics simulations using a number of different potential energy surfaces. NVE, NVT and NpT ensembles can be simulated, both periodic and nonperiodic calculations are possible.

List of examples:

• 1: Gluatathione Dihedral Angle (MD with a single QMDFF, coordinate analysis): dynamic/glutathione
• 2: Ionic Liquid with GFN2-xTB (MD of a ionic liquid pair, with GFN2-xTB): dynamic/ionic_xtb
• 3: Ferrocene with pGFN-FF (How to couple GULP for pGFN-FF calculations): dynamic/ferrocene_pgfnff
• 4: H+H2 with external program (How to connect an external PES routine to Caracal): dynamic/h3_external
• 5: Hydrogen shift on Cu-surface (Calculation of periodic systems (QMDFF)): dynamic/hydrogen_shift
• 6: Disulfide Mechanochemistry (External force MD trajectory): dynamic/mechanochemistry
• 7: Ethanol Box (Periodic system with polymerized QMDFF): dynamic/ethanol_box
• 8: H on Ni100 (How to use a neural network potential): dynamic/h_on_ni_ann
• 8: O4 on custom PES (How to couple an analytical PES with gradients to the code): dynamic/o4_custom
• 9: PO2 on custom PES (How to couple an analytical PES without gradients): dynamic/po2_custom

Folder: dynamic/glutathione

In this calculation, a single glutathione molecule is sampled in vacuum and its central peptide dihedral angle is monitored, such that the rotation between the two stable configurations can be seen as a rare event. The molecule is shown below, with the atoms along the monitored peptide angle marked by circles.

The main input file is glutatione.key, where all needed commands for the calculation are listed:

pes qmdff

qmdff {
   ffname glutathione.qmdff
   eshift -1403.0282384116763
}

xyzstart start.xyz
rpmd_beads 1
ensemble nvt 
steps 1000000
deltat 1.0
tdump 500

nvt {
   temp 500
   thermostat nose-hoover
}

eval_coord 500

The pes qmdff keyword states that a single QMDFF is used to describe the molecule. It was generated from a PBE calculation using orca (see section qmdffgen). During the QMDFF generation, the torsional force constant of the molecule are optimized using simple EHT scans of the coordinates. Therefore, the torsional landscape of glutathione can be sampled qualitatively. qmdffname gives the name of the QMDFF file present in the folder, eshift has no direct influence on the results, since it only results in a global shift in total energy. With this, the total energy is one the same scale as for the PBE reference.

The start structure start.xyz invoked by xyzstart must be present as well and contains the inital structure of the sampled molecule, in the simple xyz format. Here a classical calculation is ordered via rpmd_beads 1, you can do the calculations also with larger number of RPMD beads. The nvt ensemble is used as standard case for nonperiodic calculations (assuming infinite volume).

Since the torsional rotation barrier around the peptide bonds is quite large, roughly 1,000,000 MD steps are needed to see some of them, ordered with steps 1000000. The time step is 1 fs (deltat 1.0). Each 500 fs, a structure is written out (tdump 500).

In the section nvt, the settings for the NVT ensemble are defined: A temperature of 500 Kelvin (temp 500), and the deterministic Nose-Hoover chain thermostat for sampling of the system at the desired temperature (thermostat nose-hoover).

Finally, the keyword eval_coord 500 states that the values of certain coordinates shall be calculated and printed out every 500 fs. Which coordinates are selected for this analysis can be listed in the file coord_eval.inp, which has the following entries in this example:

8 15
2 15 8 27

This means that the bond between the atoms 8 and 15 and the peptide angle along the atoms 2, 15, 8 and 27 shall be calculated. If the position (x,y and z coordinates) of a single atom shall be printed, only one number (its index) should be given, three numbers indicate the bond between three atoms.

The calculation can be invoked by calling:

dynamic.x glutathione.key 

If everything works as it should, a long list of lines will be printed to the terminal (every 500 fs, i.e., 500 MD steps), at the end, the following lines will appear:

...
Step:   998500  --  pot. energy:  -1403.02274 Hartree  --  temperature:   427.4313 K
Step:   999000  --  pot. energy:  -1403.04119 Hartree  --  temperature:   482.8053 K
Step:   999500  --  pot. energy:  -1403.01065 Hartree  --  temperature:   574.1100 K
Step:  1000000  --  pot. energy:  -1402.99316 Hartree  --  temperature:   475.2845 K
The averaged measured temperature was:   500.48480 K.

The calculation needed a time of      60.186 seconds.
 Trajectory written to 'trajectory.xyz', temperatures to 'temperature.dat'.

.. .----. -- -.. --- -. . 
Dynamic calculation successfully finished!

Note that the exact values will differ for each run, since the initial velocities are set by a random number generator obeying the Maxwell-Boltzmann velocity distribution of the desired temperature. In this case, the desired temperature was 500 K, which was archived quite will (average temperature: 500.48480 K).

As can be seen, the trajectory of the system is stored in trajectory.xyz and might be openened with vmd or a similar program. The actual temperature of the system is stored in temperature.dat, it fluctuates around the desired average of 500 K. The fluctuations are quite large, which is due to the small system. A large box of molecules would show much smaller fluctuations.

The peptide dihedral angle was written to the file coord_eval.dat. Since it was defined in the second line of the coord_eval.inp file, it is printed here in the third column (first column: MD step). It is plotted in the figured below (in your case, the angle will probably switch on other times, or maybe not at all, depending on the Maxwell-Boltzmann initialization).

Folder: dynamic/ionic_xtb

Since version 1.1, it is possible to choose one of the GFN-xTB semiempirical methods (originating from the density functional tight binding concept) as a PES for all calculations in Caracal. The implementation is based on the tblite program (see on github), which is (with some minor changes) now included completely into Caracal. Three different Hamiltonians can be used: GFN1-xTB, GFN2-xTB and IPEA1-xTB. Periodic and nonperiodic calculations can be done, and two different implicit solvation schemes (Analytically-Linearized Poisson-Boltzmann (ALPB) and Conductor-like Polarizable Continuum Model (CPCM)) are available.

In this example, a simple short MD trajectory of a ionic liquid pair ([5-oxo-C6C1Im][NTf2]) soluted into ethanol shall be run. Main input file is dynamic.key:

pes gfn-xtb

gfn-xtb {
   hamiltonian gfn2-xtb
   solv_model alpb
   solv_spec  ethanol
  # solv_epsilon 70.0
}

xyzstart start.xyz

ensemble nvt
steps 1000
deltat 1
tdump 10

nvt {
   temp 400
   thermostat nose-hoover
}

Besides the usual dynamic keywords already introduced in the first section, the pes gfn-xtb keyword activates the family of GFN-xTB methods as PES description. Further specifications can/need to be done in the respective gfn-xtb{ section. Here, it must be stated, which of the three available Hamiltonians shall be used (GFN2-xTB) in thise case. All other keywords are optional. In this example, the ALPB implicit solvation scheme is activated, and ethanol is chosen as solvent. Further, the explicit definition of the dielectric constant of a solvent is given as commented alternative. The value 70.0 is rather close to that of water. In order to use this explicit constant, the solv_spec ethanol keyword must be commented out and solv_epsilon 70.0 be commented in. Then, the modification of the system's behavior upon change of the solvent can be studied (A complete list of all explicitly defined solvents can be seen [here](link to be added)).

Since the tblite program part is paralellized with OpenMP, it makes sense to explicitly state the number of threads before starting the calculation (e.g., 4 threads):

 export OMP_NUM_THREADS=4 
 dynamic.x dynamic.key

During the calculation (which might take roughly 5 minutes), detailed information about the SCF cycles and general timing of the xTB method is written to gfn_xtb.log. Finally, the trajectory can be inspected by opening trajectory.xyz with a suitable program.

Folder: dynamic/ionic_xtb

Since version 1.1, it is possible to couple Caracal to an existing installation of the General Utility Lattice Program (GULP) by Julian Gale (weblink). With this coupling, several other PES methods implemented in GULP can be used by Caracal. The present simple example uses the pGFN-FF force field to setup a PES for a ferrocene molecule and run a short MD trajectory with 8 RPMD beads. Since pGFN-FF has its main advantage in simulating periodic systems, this example serves mainly as description how to connect GULP to Caracal. The actual calculation setup is simple.

1: Download GULP

Go to https://gulp.curtin.edu.au/download.html and choose the latest version, 6.2 onwards (tgz file). Enter your email adress (only free for academic use, use your university domain!) and proceed. Some seconds later, you should recieve a mail with the download link. The download starts directly after clicking on it.

2: Install GULP

Extract and go into the gulp folder. Look into the README file if something has changed or should be adapted to your operating system. Then go into the Src subfolder and execute:

./mkgulp -t lib

After finishing this compilation, the pGFN-FF library needed for external calls must be built as well. For this, go back into the GULP main folder and then into the subfolder Utils/pGFNFF/Src and execute:

make

3: Modify Caracal Makefile

Now go back to the Caracal repo directory and open the Makefile in the src subfolder. There, edit the following lines:

#   Location of the GULP package, if compiled with Caracal
#   Choose either Linux (for Linux) or Darwin (for Mac)
#OSTYPE = Linux
#   The absolute or relative path to the GULP main directory
#GULPDIR =  ../../gulp-6.2
#   Comment in for the compiler flags 
#GULPFF = -DGULP=def -I${GULPDIR}/Src/${OSTYPE}
#   Comment in for the linking flags
#GULPLINK = ${GULPDIR}/Src/${OSTYPE}/libgulp.a  ${GULPDIR}/Utils/pGFNFF/Src/libpGFNFF.a 

Comment in the lines starting with OSTYPE, GULPDIR, GULPFF and GULPLINK. Then, edit the GULPDIR line and replace the dummy path there (../../gulp-6.2) by the absolute or relative path to your new GULP installation.

Then clean your current Caracal installation:

make clean

and remake everything:

make

Now, Caracal should have been patched with GULP!

4: Run the calculation

Go back into the directory of this example. There, two files exist: dynamic.key containing all commands and ferrocene.xyz containing the initial structure of the molecule. dynamic.key has the usual content:

pes pgfn-ff

pgfn-ff {
   struc_init ferrocene.xyz
}

xyzstart ferrocene.xyz

rpmd_beads 8
steps 10000
deltat 1
tdump 100
ensemble nvt

nvt {
   temp 298
   thermostat nose-hoover
}

The pGFN-FF is invoked by the pes pgfn-ff main keyword, the respective pgfn-ff section only contains the initial structure: struc_init ferrocene.xyz. The structure given with this keyword is used for the pGFN-FF initialization (buildup of bonding matrix etc.). The actual start file for the MD trajectory is given by xyzstart ferrocene.xyz. Although both files are identical here, another structure could be given as well. Note that at least the bonding configuration should be kept the same for a useful description, since this is not updated during the dynamics. The remaining keywords are the same as usual. Now start the MD trajectory with:

dynamic.x dynamic.key

The positions of all beads in each frame are written to trajectory.xyz, those of the centroids to traj_centroid.xyz.

Folder: dynamic/h3_external

This example mainly serves as a manual to connect an external PES calculation program to Caracal. With the presented interface it is possible to connect every available program that calculates the energy and the cartesian gradient of a structure. Note, however, that this method of external calling is quite inefficient. Especially if very fast programs or routines shall be connected, the needed time for file IO will be much longer than the actual external calculation time! If, e.g., a publicly available analytical PES shall be used with Caracal that is not already implemented, please contact me! I might add this routine to the program package itself.

This example case simply calculates a short trajectory of the H+H2 reaction, starting at the TS. It uses the well known H+H2 function by Truhlar and Horowitz in an updated version (DOI).
Since this surface is already implemented within Caracal (ANA_H3), this is rather pointless and shall mainly serve as a tutorial how to connect an external program/subroutine to Caracal.

The dynamic.key file looks as follows:

pes external

xyzstart start.xyz
rpmd_beads 32

ensemble nvt
steps 1000

deltat 0.1
tdump 1.0

nvt {
   temp 200
   thermostat andersen
}

external {
   symlink ./external.x 
}

The pes external keyword activates the external call. After the usual dynamic keywords for a short unbiased NVT trajectory, the external section with the symlink keyword needs to be added. In this case, the external program is called with the command ./external.x.

This program is a short Fortran90 program, serving as handler for the H+H2 PES, which is only a subroutine. The Fortran file (call_external.f90) has the following lines:

program call_external
implicit none
real(kind=8),allocatable::xyz(:,:),grad(:,:)
integer::i,j,natoms
real(kind=8)::energy
character(len=3)::adum
integer::info

open(unit=37,file="coords.xyz",status="old")  
read(37,*) natoms

allocate(xyz(3,natoms))
allocate(grad(3,natoms))

read(37,*)
do i=1,natoms
   read(37,*) adum, xyz(:,i)
end do
close(37)
xyz=xyz/0.52917721092d0


call potential(xyz,natoms,energy,grad,info)

open(unit=38,file="egrad_out.dat",status="replace")
write(38,*) energy
do i=1,natoms
   write(38,*) grad(:,i)
end do

end program call_external

The procedure works as follows:

1. Caracal (dynamic.x) writes a coords.xyz file in the usual xyz format with the coordinates in Angstrom.
2. Caracal calls the external program (external.x in this case).
3. Caracal reads the results file egrad_out.dat, where the energy (Hartree) is printed in the first line and the gradient components are written: one atom per line, x-,y- and z- components in Hartree/bohr.
4. Caracal uses the obtained energy/gradient for its purpose: in this case, 32 calls are performed for each MD step (one for each bead), then the next time step is calculated.

The external program can be written with every programming language (or simply be a shell script), as long as it reads in the coords.xyz file and writes the calculated energy and gradient to the egrad_out.dat file!

In the current example, the external program (call_external.f90 + h3_external.f) must be compiled first, e.g., with gfortran:

gfortran -o external.x call_external.f90 h3_external.f

Then, the dynamic.x program can be executed as usual:

dynamic.x dynamic.key

The calculation will take some seconds (which is some orders of magnitudes slower than the internal ANA_H3 PES!) and print out the usual output (trajectory.xyz etc.).

Folder: dynamic/hydrogen-shift

This example covers the calculation of a larger system: A cutout of a copper surface, on which a hydrogen atom is placed, containing 81 atoms in total. A simple dE-EVB-QMDFF was set up, where the two QMDFFs each describe the hydrogen being adsorbed on a hollow position on the surface. The EVB is parametrized to enable a correct description of the transition between those to (neighbored) hollow positions. This subdivided description is shown in the figure below. (QMDFF1: left minimum, QMDFF2: right minimum, dE-EVB: transition region)

As usual, the main input file is the .key file, shift.keyin this example:

pes de_evb
qmdffnames min1.qmdff min2.qmdff
2evb  -3850.13597864118219    -3850.13586849674994

de_evb {
   function 1g
}

xyzstart ts.xyz
rpmd_beads 64
ensemble nvt
steps 2000
deltat 0.5
tdump 10.0

nvt {
   temp 50
   thermostat nose-hoover
}

fix_atoms fix_atoms.dat
eval_coord 10

In the first section, the potential energy description is set up. A preoptimized dE-EVB is the model (pes de_evb), with the QMDFFs min1.qmdff and min2.qmdff building the diabatic surfaces to be coupled. The absolute reference energies (resulting from a CP2K calculation) of both QMDFFs are given with the 2evb keyword, and are now of central importance (at least their relative values).

More details concerning the dE-EVB coupling term are given in the de_evb { section, which here only contains the definition of the simple 1g coupling function (function 1g) (a Gaussian function, dependent on the local energy difference of both QMDFFs). The parameters defining the coupling term are stored in the file evb_pars.dat, which was generated from evbobt. It contains just two parameters: the prefactor of the Gaussian (a) and the exponential prefactor of the Gaussian (b).

 ** This file contains parameters for a EVB-dE-calculation **
   2.4754341517693738E-002
   512.90654128480833  

The following keywords are the usual dynamic keywords, similar to the first example. The main difference is that now 64 RPMD beads are used (bead_number 64). This allows for the description of dynamic quantum effects, since especially the hydrogen atom will be modeled not as a single point particle, but as a somewhat spatially smeared entity. A low temperature of just 50 K was chosen in order to make the quantum effects visible (more delocalization at low temperatures).

Since we model a surface slab, which is thought to be attached to a metal bulk beneath it, it makes sense to somehow fix the lower metal atoms. This is done via the fix_atoms keyword, which freezes the copper atoms of the lower two layers. For this, the numbers of all atoms to be fixed must be stored into a separate file: fix_atoms.dat. This file contains one atom index per line (no ascending order required):

1
2 
3
4
5
6
... 
36
37
38
39
40

For the purpose of subsequent analysis, the keyword eval_coord is added again. The respective coord_eval.inp file now contains only one single number:

81

With this, the position of the hydrogen atom (being the last one) is monitored, x-, y- and z-coordinates are written to the file coord_eval.dat every 10 fs. It should be noted, that only the position of the centroid, i.e., the average value of all RPMD beads is monitored by eval_coord, since only the centroid has a direct connection to experimental observables.

If we look at a snapshot of the trajectory during the dynamics, it will be seen that the hydrogen atom is delocalized significantly, probably even over the barrier into both minima at once, whereas the 64 RPMD beads of the much heavier Cu atoms still cover almost the same positions:

We can further inspect the position of the hydrogen centroid over the whole trajectory. This can be done by plotting the cartesian positions in coord_eval.dat in a 3D plot, with x and y as independent coordinates:

The area in the top left has a slightly higher z coordinate, corresponding to the valley betweeen the hollow positions, where the hydrogen atom is located at the beginning (ts.xyz). Probably due to the smearing of its position into both valleys at a time, it remains some 100 fs in the TS region, until it finally falls into the minimum in the bottom right, described by QMDFF2.

Folder: dynamic/glutathione

The dynamic program is able to apply forces on one or two atoms during MD simulations. One special mode covers the simulation of a atomic force microscope (AFM) experiment. In this, the molecule of interest is clamped between a carrier surface, which stays at constant height during the experiment and the AFM tip, which moves one part of the molecule away from the surface. At some time, the strain will be so strong that the molecule will rip apart.

In our example, we have a disulfide molecule which has the special feature of ripping apart in two parts:

By pulling on the red atom on the right, the strain will first be focused on the central disulfide bridge, which is weaker than the C-C bonds and will therefore rip apart first. Then, the remaining chain will be strained as well, and at some stage a C-C bond will break, resulting in complete destruction of the molecule.

In order to carry out such a simulation, the following commands as listed in the mechano.key file are needed:

 pes qmdff

 qmdff {
    ffnames  freq.qmdff
    eshift -2377.3292550041192
 }

 xyzstart min.xyz
 rpmd_beads 1
 steps 500000
 deltat 0.5
 tdump 50
 ensemble nvt

 nvt {
   temp 200
   thermostat  nose-hoover
 }

 force {
   afm_run
   afm_fix 6
   afm_move 70 35.0 25.2965 0.02 -1.171
   afm_avg 20
   afm_second 100000
 }

A single QMDFF is a sufficient PES description of the whole process, since it is semireactive and allows for the breaking of covalent bonds, such as the disulfide or one of the C-C bonds. Only if a further reaction involving bridging of new bonds should be described as well, an EVB coupling term would be needed. We will perform a classical MD (bead_number 1), within a NVT ensemble, taking 500000 time steps, in order to get a smooth description of the process.

Besides the already known nvt section specifying the properties of the ensemble, a new section named force must be added in order to set up the AFM simulation. The keyword afm_run activates this special kind of force application. afm_fix 6 specifies that the 6th atom (being the green marked in the figure above) shall be hold fixed. Atom 70 (the red marked) is moved with the keyword afm_move 70 35.0 25.2965 0.02 -1.171. Here, 70 specifies the index of the atom, 35.0 says that it will be moved 35 Angstroms in total along the vector 25.2965 0.02 -1.171. No force must be given, since the atom will be moved straightforward, with an "infinite force". In order to make the system able to react smoothly to this somewhat artificial modification, a long trajectory is recommended if a large molecule like this disulfide shall be simulated in the AFM mode.

The remaining two keywords affect only the evaluation of the calculation. afm_avg 20 orders the program to average the force imposed by the molecule on the moved atom (70) every 20 steps in order to reduce thermal noise. This averaged force (in nN) will be written to file afm_averages.dat, where also the averaged length of the molecule (distance between atoms 6 and 70 in this case) will be written every 20 time steps. After the simulation, the maximum applied force during the simulation will be written to the screen, together with the molecule length at this point. If one is interested into a two-step reaction (like in this case), the keyword afm_second which a number of time steps orders the program to do this analysis before and after the proposed number of time steps, such that two reactive events are monitored independently. This analysis can of course also be done manually by inspecting the file afm_averages.dat, such that this is only a matter of convenience.

After running the simulation by invoking

dynamic.x mechano.key

the following information (or similar) will be written to the screen:

... 
Step:   499900  --  pot. energy:  -2377.00520 Hartree  --  temperature:   191.5480 K
Step:   500000  --  pot. energy:  -2377.01925 Hartree  --  temperature:   196.6160 K
Averaged values for AFM forces and distances were
written to file 'afm_averages.dat'.

The first reaction is located at     10.3283731 nN at     25.1564797 Angstrom.
The second reaction is located at     14.2123746 nN at     38.6247096 Angstrom.
The averaged measured temperature was:   199.86418 K.

The calculation needed a time of      77.917 seconds.
 Trajectory written to 'trajectory.xyz', temperatures to 'temperature.dat'.

.. .----. -- -.. --- -. . 
Dynamic calculation successfully finished!

Therefore, it required a force of roughly 10 nN to break the disulfide bond (where the molecule had a length of 25 Angstroms). The second bond breaking, a C-C bond, required a larger force of roughly 14 nN, the molecule was elongated to 39 Angstroms at this point. If we plot the file afm_averages.dat, the time-dependent force can be seen explicitly:

Besides the AFM mode, dynamic is also able to apply force vectors in two atoms in the system. In this case, no atom is moved with a constant speed, but they are torn apart with a certain strength, where it is not certain that it will result in the desired bond rupture. A input file for such a calculation is placed in the folder as well. It is named forcevec.key and has the following entries:

pes qmdff

qmdff {
   ffnames  freq.qmdff
   eshift -2377.3292550041192
}

xyzstart min.xyz
rpmd_beads 1
steps 1000000
deltat 0.5
tdump 100
ensemble nvt

nvt {
  temp 335
  thermostat  nose-hoover
}

force {
  vec1 6 2.520E-9 35 25.2965 0.02 -1.171
  vec2 70 2.520E-9   -25.2965 -0.02 1.171
}

eval_coord 50

Most keywords are left unchanged. The number of time steps was doubled and the temperature raised such that a rupture event can be seen. In the force section, both applied forces are specified separately. Their structure is as follows: Index of atom, on which force shall be applied; magnitude of the force (in Newton), and the direction of the force vector. Here, the atoms 6 and 70, serving as anchors in the AFM simulation, were chosen again, forces of 2.52 nN are applied in both directions, summing up to 5.04 nN. This is only half the maximum value measured in the AFM experiment, a smaller value is needed here since the force is applied during longer times and the the atoms are accelerated at the beginning, leading to a larger effective force acting on the bonds when the molecule is completely stretched first. The force vectors are the same as in the afm_move keyword, but now in opposite directions, avoiding a movement of the molecules center of mass. After executing the dynamic.x program, the following output will be written:

...
Step:   999800  --  pot. energy:  -2377.05097 Hartree  --  temperature:   319.7660 K
Step:  1000000  --  pot. energy:  -2377.06109 Hartree  --  temperature:   331.4498 K
Mechanochemistry results:
 Length of the molecule at the beginning:                      23.090453 Angstrom.
 Length of the molecule at the end:            33.865335 Angstrom.
The averaged measured temperature was:   335.24351 K.

The calculation needed a time of     157.234 seconds.
 Trajectory written to 'trajectory.xyz', temperatures to 'temperature.dat'.

.. .----. -- -.. --- -. . 
Dynamic calculation successfully finished!

If we look at the lengths of the molecule (defined again as the distances between both atoms on which the forces are applied), it can directly be seen that the disulfide bond was ripped apart, leading to an elongation of around 10 Angstroms. Further insight can be gained by plotting the coord_eval.dat file, which was written here by the eval_coord keyword. The file coord_eval.inp contains just one line:

6 70

It thus measures the distance between both active atoms. The plot will look something like this:

The rupture of the disulfide bond happens shortly after 0.4 ns. When (and if) the rupture happens during the 1 ns simulation time is different for each run, since the initial velocities of the atoms are initialized randomically. Raising the applied force will lead to almost direct rupture, since the effective force is especially large during the first fs, when the chain is somewhat elongated from its initial state and an additional velocity acts on the disufide bond. Raising the temperature will also lead to faster rupture, since more kinetic energy is available within the degrees of freedom to be added to the external force.

Not that also much smaller applied forces will lead to ruptures, but those might happen after several 100s of ns or even later, since the activation barrier grows exponentially with declining force.

Folder: dynamic/ethanol_box

In this example, an ethanol box built by poly_qmdff in advance (see tutorial) is equilibrated by an NpT trajectory in a classical calculation (one RPMD bead) and subsequently sampled with NVT and 16 RPMD beads.

Both QMDFF file box.qmdff and initial structure box.xyz are taken directly from the poly_qmdff calculation. For both the equilibration and the sampling, keyfiles are located in the folder as well.

In the equilibration.key has the following content:

 pes qmdff

 qmdff {
    ffnames box.qmdff
    eshift 0.0
 }

 xyzstart box.xyz
 rpmd_beads 1
 steps 10000
 deltat 0.5
 tdump 10.0

 ensemble npt

 npt {
    temp 200
    pres 1.0
    thermostat nose-hoover
    barostat nose-hoover
    baro_damp 10000
    periodic 27.0 27.0 27.0  
 }

The equilibration takes 10000 time steps (steps 10000) which is rather short and might be raised if enough time is available (e.g., to 40000). All other settings are similar to that shown above.

Now, however, the NpT ensemble is used (ensemble npt). The respective settings for this ensemble are located in the npt section. Besides the temperature control, additional keywords concerning the pressure are present: pres 1.0 sets the pressure to one atmosphere, barostat nose-hoover sets the Nose-Hoover chain barostat, baro_damp sets the barostat damping parameter (how fast the reaction of a difference between actual and desired pressure will be). Too small values might lead to explosions of the system, wheras too large values raise the needed number of time steps. The keyword periodic 27.0 27.0 27.0 activates a periodic calculation, which is always needed for a NpT calculation but is optional for NVT calculations (assuming an infinite volume if not given). The dimensions are also taken from poly_qmdff which prints the dimensions of the generated box.xyz file to the command line.

Invoke the calculation by typing:

dynamic.x 

Since tdump 10 is set, an information line will be printed every 10 time steps to the command line, containing the potential energy, actual temperature, volume, pressure and density:

 Performing MD ....
   Step    Epot(Hartree)    T(K)      Volume(A^3)     Press.(atm)   density(g/cm^3)
      20    -0.16509       128.9973    19685.3932      32.54611        0.485476
      40    -0.00554       103.3309    19670.0749     102.48462        0.485855
 ....
    9980    -1.00975       206.6086    10904.1019    -168.26176        0.876440
   10000    -0.98194       201.3728    10911.1905    -190.58679        0.875871
 The averaged measured temperature was:   199.09688 K.

Folder: dynamic/o4_custom

In this example, it will be shown how you can manually integrate your own (Fortran) potential energy surface subroutine to Caracal. I tried to make this whole process as simple and straightforward as possible. Nevertheless, this example is much more involved than all others shown here, and of course each routine might have its own features that must be considered (units of input/output variables, initialization needed or not, etc.)

Here, a PES by Paukku et al. (see paper), obtained from the POTLIB page describing the adiabatic ground state of the O2+O4 quintet state shall be coupled to Caracal and a simple unbiased dynamical trajectory shall be calculated as a benchmark.

The following steps need to be taken:

1: Add the PES source code file to the Caracal code

The file o4_quintet.f need to be copied into the Caracal src folder. Then, it must be added to the Makefile such that it can be compiled. Open it and modify the OBJS section by adding the entry o4_quintet.o, for example like this:

   ...
   lm_qmdff.o water_init.o egrad_water.o ev_coord_init.o \
   external_grad.o custom_init.o custom_grad.o \
   \
   egrad_ch4oh.o util_ch4oh.o egrad_h3.o \
   ...

to

   ...
   lm_qmdff.o water_init.o egrad_water.o ev_coord_init.o \
   external_grad.o custom_init.o custom_grad.o \
   o4_quintet.o \
   \
   egrad_ch4oh.o util_ch4oh.o egrad_h3.o \
   ...

2: Implement an energy/gradient call to the routine

For this, stay in the Caracal src folder and open the file custom_grad.f90. As you can see, it is possible to add several (in principle infinitely many) external subroutines to the Caracal code! Which one is used for a calculation is then defined by the cust_number variable (see below).

Since our routine calculates analytical gradients directly, we can ignore the num_grad section. Instead, you need to remove the comment before the gradient call line in the cust_number .eq. 1 section:

!
!     If the external routine is able to calculate the energy and analytical 
!     (or numerical gradients), activate this call here and adapt it to your 
!     routine!
!
!    call your_routine(xyz2,e_evb,g_evb)

to

!
!     If the external routine is able to calculate the energy and analytical 
!     (or numerical gradients), activate this call here and adapt it to your 
!     routine!
!
call your_routine(xyz2,e_evb,g_evb)

Next, open the subroutine file o4_quintet.f. Here, the sole relevant part is the first subroutine:

C   Input: X(4),Y(4),Z(4)               in unit of bohr
C   Output: E                           in unit of hartree
C   Output: dEdX(4),dEdY(4),dEdZ(4)     hartree/bohr
C
      subroutine pot(X,Y,Z,E,dEdX,dEdY,dEdZ)

      implicit double precision (a-h,o-z)

C
C Convert factors (Use conversion factors compiled on Dec. 5)
C Cconv: bohr to Angstrom 
C        1 bohr = 0.52917721092 angstrom
C Econv: kcal/mol to hartree 
C        1 kcal/mol = 0.159360144 * 10^-2 hartree
C Gconv: kcal/(mol*Angstrom) to hartree/bohr
C        1 kcal mol^-1 angstrom^-1 = 0.843297564 * 10^-3 hartree/bohr
C
      double precision Cconv
      double precision Econv
      double precision Gconv
      parameter(Cconv=0.52917721092d0)
      parameter(Econv=0.159360144d-2)
      parameter(Gconv=0.843297564d-3)
C
C Reference energy of infinitely separated O2 + O2 in hartree (taken
C from
C calculations)
C
      double precision Eref
      parameter(Eref=-299.842500d0)

      integer i
      double precision E,V
      double precision X(4),Y(4),Z(4),dEdX(4),dEdY(4),dEdZ(4)
      double precision Xcart(12),dVdX(12)

C Convert to local variables
      do i=1,4
        Xcart(3*i-2)=X(i)*Cconv
        Xcart(3*i-1)=Y(i)*Cconv
        Xcart(3*i)=Z(i)*Cconv
      enddo

      call o4pes(Xcart,v,dVdX,1)

C Convert local output to the ones ANT wants

      E = v*Econv + Eref
      do i=1,4
        dEdX(i)=dVdX(3*i-2)*Gconv
        dEdY(i)=dVdX(3*i-1)*Gconv
        dEdZ(i)=dVdX(3*i)*Gconv
      enddo

      end

The first three parameters are the input x-, y- and z coordinates of the four atoms as separate vectors, followed by the energy and the x-,y- and z coordinates of the resulting gradient as output. All parameters need to be in atomic units, which is also the case for Caracal. If your subroutine needs other units, e.g., Angstrom for the input coordinates, you can use one of the global unit conversion parameters explained in the header section of custom_grad.

Therefore, the connection is quite straightforward: Simply submit/receive the x-,y- and z coordinate columns of the coordinates and gradient separately. The call in custom_grad.f90 should then look like this:

call pot(xyz2(1,:),xyz2(2,:),xyz2(3,:),e_evb,g_evb(1,:),g_evb(2,:),g_evb(3,:))

Optionally, it might be wise to rename the subroutine, since pot is not quite a unique name... Especially in the case that you want to add several subroutine at once.

Some subroutine require a separate initialization routine (e.g., for read in of external files with parameters), which is executed only once after calling the program. Open the file custom_init.f90 to add the routine to Caracal (in the section of the same cust_number!)

3: Recompile Caracal

After all changes are implemented, recompiling the code and making the executables should be no problem:

make 

4: Run the calculation

Now switch back to the folder of this example. Besides the source code for the subroutine, a starting structure start.xyz and a keyfile benchmark.key are located there. start.xyz contains two oxygen molecules, separated by 4 Angstroms. The keyfile has the following content:

pes custom

xyzstart start.xyz
rpmd_beads 4

ensemble nvt
steps 10000
deltat 0.5
tdump 50.0

nvt {
   temp 300
   thermostat andersen
}

custom {
   pes_number 1 
}

Caracal is ordered to use a manually added PES by the general PES definition pes custom. In the custom section at the end of the file, the specific PES is chosen with the pes_number keyword (currently 1, since only one PES has been added in this example).

Start the dynamics trajectory with

dynamic.x bechmark.key

and the calculation should run through in less than a second, emitting the usual output files (see previous sections).

Folder: dynamic/po2_custom

This example is essentially a continuation of the previous one, therefore, the explanations will be somewhat shorter. Here, another external subroutine shall be added to Caracal, but now, the task is somewhat more challenging, since the PES does not provide an analytical gradient (This example can be combined with example 6, such that two custom PESs are added to Caracal. Then, everything still works as noted here, but the gradient calls must be added to the else if (cust_number .eq. 2) then section into custom_grad.f90. Simply copy and pase the numerical gradient code from the first section and modify it as described below.) The provided PES by Chen el al. (see paper) describes the ground state of the PO2 molecule and its decomposition pathways (to P + O2 and to O + PO)

The following steps need to be taken:

1: Add the PES source code file to the Caracal code

Copy the file pes_po2.f to the src folder of Caracal and add it to the Makefile.

   ...
   lm_qmdff.o water_init.o egrad_water.o ev_coord_init.o \
   external_grad.o custom_init.o custom_grad.o \
   \
   egrad_ch4oh.o util_ch4oh.o egrad_h3.o \
   ...

to

   ...
   lm_qmdff.o water_init.o egrad_water.o ev_coord_init.o \
   external_grad.o custom_init.o custom_grad.o \
   pes_po2.o \
   \
   egrad_ch4oh.o util_ch4oh.o egrad_h3.o \
   ...

2: Implement an energy/gradient call to the routine

Open the pes_po2.f file. The main subroutine is placed directly at the top:

      SUBROUTINE PES_PO2(R1,R2,R3,ENER)
c       R1=roo;R2=rpo;R3=rpo
      IMPLICIT  REAL * 8 (A-H,O-Z)
      DOUBLE PRECISION :: R1,R2,R3,V3
      DOUBLE PRECISION, DIMENSION(3) :: R
      DOUBLE PRECISION :: E12,E13,E23


      R(1)=R1
      R(2)=R2
      R(3)=R3
      CALL CHIPR_O2(R1,E12)
      CALL CHIPR_PO(R2,E13)
      CALL CHIPR_PO(R3,E23)
      CALL TRIAT(R1,R2,R3,V3)
      ENER=E12+E13+E23+V3
      RETURN
      END

It can be seen that the subroutine requires three interatomic distances as input (R1, R2 and R3) and gives the current energy as output (ENER). Thus, there is no gradient calculated. Further, it does not accept a xyz structure as provided by Caracal as input! Therefore, the subroutine need to be changed somewhat into the following (one possibility):

      SUBROUTINE PES_PO2(xyz,ENER)
c       R1=roo;R2=rpo;R3=rpo
      IMPLICIT  REAL * 8 (A-H,O-Z)
      DOUBLE PRECISION :: R1,R2,R3,V3
      DOUBLE PRECISION, DIMENSION(3) :: R
      DOUBLE PRECISION, DIMENSION(3,3) :: xyz
      DOUBLE PRECISION :: E12,E13,E23

      R1 = sqrt(dot_product(xyz(:,2)-xyz(:,3),xyz(:,2)-xyz(:,3)))
      R2 = sqrt(dot_product(xyz(:,1)-xyz(:,2),xyz(:,1)-xyz(:,2)))
      R3 = sqrt(dot_product(xyz(:,1)-xyz(:,3),xyz(:,1)-xyz(:,3)))

      R(1)=R1
      R(2)=R2
      R(3)=R3
      CALL CHIPR_O2(R1,E12)
      CALL CHIPR_PO(R2,E13)
      CALL CHIPR_PO(R3,E23)
      CALL TRIAT(R1,R2,R3,V3)
      ENER=E12+E13+E23+V3
      RETURN
      END

First, the arguments of the subroutine were changed, R1, R2, R3 were replaced by xyz. This array is then defined in the declaration section: DOUBLE PRECISION, DIMENSION(3,3) :: xyz. The first dimension always are the x-, y- and z coordinates of the current atom, the second dimension is the index of the atom. Now, the needed interatomic distances must be calculated from the given structure. This is done directly after the variable declaration. According to the given comment describing the distances (with P being the first atom), the O-O distance is calculated by: R1 = sqrt(dot_product(xyz(:,2)-xyz(:,3),xyz(:,2)-xyz(:,3))) and so on.

Then, open the file custom_grad.f90 and add the gradient routine. There, add the subroutine to the code snippet for numerical gradients:

!   call your_routine(xyz2,e_evb)

   step=num_grad_step
   do i=1,natoms
      do j=1,3
         do k=1,2
            if (k.eq.1) then
               xyz2(j,i)=xyz2(j,i)+step
            else if (k.eq.2) then
               xyz2(j,i)=xyz2(j,i)-2*step
            end if
!
!     Activate this call here and adapt it to your routine!
!
!               call your_routine(xyz2,e_tmp)


            if (k.eq.1) then
               e_upper=e_tmp
            else if (k.eq.2) then
               e_lower=e_tmp
            end if
         end do
         g_evb(j,i)=(e_upper-e_lower)/(2*step)
         g_evb(j,i)=g_evb(j,i)*hartree*bohr
         xyz2(j,i)=xyz2(j,i)+step
      end do
   end do
end if

Two calls can be seen here. The first call calculates the energy of the undistorted structure (e_evb). It must be changed accordingly:

call pes_po2(xyz2,e_evb)

The second call inside the do-loop calculates the actual energy for each geometry distortion and stores it into e_tmp. Those values are then postprocessed to obtain the gradient components. The call need to be changed to:

call pes_po2(xyz2,e_tmp)

The subroutine works with atomic units, therefore, no unit conversions need to be done here.

3: Recompile Caracal

After all changes are implemented, recompiling the code and making the executables should be no problem:

make 

4: Run the calculation

Again, the input for a short MD trajectory is given in the folder of this example. Here, the start structure is given in start.xyz and the keywords for dynamic are given in benchmark.key:

 pes custom
 num_grad

 xyzstart start.xyz
 rpmd_beads 4

 ensemble nvt
 steps 10000
 deltat 0.5
 tdump 20.0

 nvt {
    temp 500
    thermostat andersen
 }

 custom {
    pes_number 1
 }

Again, the newly implemented PES is activated by the keyword pes custom, and the first implemented PES is chosen by pes_number 1 in the custom section. (If you already added the PES from example 6, the keyword should now be pes_number 2)

Further, it is important to add the keyword num_grad if the gradient needs to be calculated numerically! Else, the trajectory would be useless since no gradient would be calculated at all.

Start the calculation as usual:

dynamic.x benchmark.key

and the usual output files are obtained. If everything worked correctly, the PO2 molecule should stay stable and vibrating a bit, with the single RPMD beads being almost completely merged. Since the gradient is calculated numerically, the calculation takes somewhat longer, approximately 10 seconds.