MD Tutorial
In this tutorial you will learn how to perform a molecular dynamics simulation in CAST and how to deal with the different options for this task.
External programs you need: VMD
- First example: NVE MD
- Temperature stuff
- Velocity scaling
- Additional Potentials
- Trackoffset
- Snapshot buffer
- Optimization of Snapshots
- Pre-Optimization
- Refinement
- Rattle
- Restart broken simulations
- Pressure control
- Conclusion
Files needed for this example:
- structure file: tridecaalanin.arc
- forcefield parameter file: oplsaa.prm (download from tinker homepage)
- CAST inputfile: CAST.txt
This is your CAST.txt file:
#### GENERAL INFORMATION ####
verbosity 4 # how much information is printed
name tridecaalanin.arc # name of the structure file
outname output # name of the output file(s)
inputtype TINKER # inputfile type (in this case: tinker)
#### TASK ####
task MD # task: MD simulation
#### OPTIONS FOR ENERGY INTERFACE ###
interface OPLSAA # energy interface: OPLSAA forcefield
paramfile oplsaa.prm # name of forcefield parameter file
cutoff 12.0 # cutoff for non-bonded interactions
switchdist 10.0 # additional option for cutoff
#### OPTIONS FOR MD TASK ####
MDsteps 50000 # number of MD steps
MDintegrator 0 # which integrator? (0 = velocity verlet)
MDtimestep 0.001 # time of one MD step (in ps)
MDtrack 1 # activate tracking (otherwise no outputfiles)
MDsnap 1000 # number of snapshots
Attention! The stuff behind the # is only for explanatory reasons. Please remove it before starting the calculation!
As outputfiles you get:
- output_MD_TRACE.csv: table with timestep, temperature, pressure (not calculated), kinetic, potential, total energy and snapshot number, can be opened for example with LibreOfficeCalc
- output_MD_SNP.txt: file with all structures that are snapshots, can be viewed in VMD
What we have done here is to take a molecule, don't give it any temperature (i. e. none of the atoms has any velocity) and then let it move at a constant energy level. This means that in the beginning the molecule has only potential energy, no kinetic energy. During the simulation the potential energy gets lower (as the systems moves downwards on the potential energy surface) and the kinetic energy gets higher whereas the total energy remains constant. This can be seen nicely in "output_MD_TRACE.csv".
Now we want the same behaviour as before (so we want a NVE ensemble where the total energy is constant) but we don't want to start at 0 K but at another temperature. So we add the following lines to our inputfile:
MDtemp_control 0
MDheat 0 300
With the first line we switch off temperature control (if we wouldn't CAST would keep the temperature at 300 K for the whole simulation time, see next section). The second line gives the starting temperature (300 K at step 0, i. e. at the beginning). If we take a look at "output_MD_TRACE.csv" we see that the simulation starts at around 300 K but then the temperature changes quite rapidly. The total energy however remains constant.
In order to activate temperature control add the following lines to your CAST.txt file:
MDtemp_control 1
MDheat 1000 50
MDheat 2000 300
MDheat 3000 100
With the first line we activate temperature control, i. e. the temperatures we give below are not just for the first step, but for the whole simulation.
This simulation starts with 50K, keeps its temperature constant till 1000 steps, then rises it linearly to 300K at 2000 steps, decreases it back to 100K at 3000 steps and keeps it constant at 100K till the end of the simulation. This temperature behaviour can be seen very nicely in the file "output_MD_TRACE.csv" as it is done by direct velocity scaling which means that the temperature is scaled exactly to the desired temperature in every step. Of cause this means that energy is put it on taken out of the system so the total energy is not constant anymore.
For switching on the Nosé-Hoover thermostat add the following line to your inputfile:
MDthermostat 1
Now the same temperature evolution is done as in the section before but the temperature control is done by Nosé-Hoover thermostat instead of direct velocity scaling. As a consequence of this the temperature varies a bit around the desired value but it is more like a real physical behavior.
With the option MDnosehoover_Q you can manipulate the Nosé-Hoover thermostat. The default value for it is 0.1. In order to show the effect of Q we start two simulations. For the first one we just make our former simulation shorter and let the temperature raise from 0 to 300 K over the simulation time of 1000 steps:
MDsteps 1000
MDheat 0 0
MDheat 1000 300
For the second we just add one more line:
MDnosehoover_Q 0.5
Now we plot the temperature evolution during the simulations for the two different Q's:
So when Q is higher the temperature seems to react more slowly.
In this section we will see what the option MDveloscale does. For this purpose we use a structure of ethane that looks like that:
8 (oplsaa.prm Parameter)
1 C 8.071103 -2.620014 -3.714228 80 2 3 4 5
2 H 8.147952 -1.583782 -3.321000 85 1
3 H 9.086855 -3.054936 -3.829940 85 1
4 H 7.574121 -2.549369 -4.705345 85 1
5 C 7.133741 -3.345059 -2.763358 80 1 6 7 8
6 H 6.871486 -4.375160 -3.086501 85 5
7 H 6.162287 -2.808674 -2.709798 85 5
8 H 7.547320 -3.311831 -1.732754 85 5
Just copy the above into a textfile and rename it "ethane.arc".
This is our CAST.txt file:
verbosity 4
name ethane.arc
outname no_veloscale
inputtype TINKER
task MD
interface OPLSAA
paramfile oplsaa.prm
cutoff 12.0
switchdist 10.0
MDsteps 5000
MDintegrator 0
MDtimestep 0.001
MDtrack 1
MDsnap 1000
MDthermostat 1
MDheat 0 300
MDveloscale 0
Notice that the option MDveloscale is still switched off. Now we do the same calculation with this option switched on. So we change the name of the outputfile (outname) to veloscale and the option MDveloscale to 1 and run CAST again. Now we compare the trajectories "no_veloscale_MD_SNAP.arc" and "veloscale_MD_SNAP.arc" by opening both of them in VMD.
We see that in the trajectory without veloscale the molecule moves outside of the window whereas in the trajectory with velocity scaling it stays in the same place. This is because the option MDveloscale removes the total translation and rotation of the molecule, both before the simulation starts and after every MD step. As normally this is what we want it is switched on by default if you don't give this option.
Of course you can also use the general bias potentials that are available for every task of CAST, i. e. spherical and cubic potentials and potentials on individual bonds, angles and torsions. For more information see our manual (LaTeX file). Here only those biases will be described that are only for task MD.
Now we will add a spherical potential to our MD. As example we use a cluster of 27 argon atoms. The tinkerstructure of the system ("argon.arc") looks like that:
27 (oplsaa.prm Parameter)
1 Ar 0.000000 0.000000 0.000000 60
2 Ar 0.000000 0.000000 3.000000 60
3 Ar 0.000000 0.000000 6.000000 60
4 Ar 0.000000 3.000000 0.000000 60
5 Ar 0.000000 3.000000 3.000000 60
6 Ar 0.000000 3.000000 6.000000 60
7 Ar 0.000000 6.000000 0.000000 60
8 Ar 0.000000 6.000000 3.000000 60
9 Ar 0.000000 6.000000 6.000000 60
10 Ar 3.000000 0.000000 0.000000 60
11 Ar 3.000000 0.000000 3.000000 60
12 Ar 3.000000 0.000000 6.000000 60
13 Ar 3.000000 3.000000 0.000000 60
14 Ar 3.000000 3.000000 3.000000 60
15 Ar 3.000000 3.000000 6.000000 60
16 Ar 3.000000 6.000000 0.000000 60
17 Ar 3.000000 6.000000 3.000000 60
18 Ar 3.000000 6.000000 6.000000 60
19 Ar 6.000000 0.000000 0.000000 60
20 Ar 6.000000 0.000000 3.000000 60
21 Ar 6.000000 0.000000 6.000000 60
22 Ar 6.000000 3.000000 0.000000 60
23 Ar 6.000000 3.000000 3.000000 60
24 Ar 6.000000 3.000000 6.000000 60
25 Ar 6.000000 6.000000 0.000000 60
26 Ar 6.000000 6.000000 3.000000 60
27 Ar 6.000000 6.000000 6.000000 60
First we perform an MD with the spherical potential switched off:
verbosity 4
name argon.arc
outname no_spherical
inputtype TINKER
task MD
interface OPLSAA
paramfile oplsaa.prm
cutoff 12.0
switchdist 10.0
MDsteps 50000
MDintegrator 0
MDtimestep 0.001
MDtrack 1
MDsnap 1000
MDspherical 0 10.0 15.0 10.0 20.0 2.0 4.0
Then we switch it on by changing the first number inMDspherical to 1. We also change the outname to spherical and run CAST again. Then we compare the 2 trajectories "no_spherical_MD_SNAP.arc" and "spherical_MD_SNAP.arc" in VMD.
Whereas the atoms constantly move away from each other without the spherical potential they are kept close to each other when the spherical potential is moved on.
MDspherical even applies two additional spherical potentials which you can specify by 7 parameters:
- <0/1> switches on and off spherical potential
- radius where the first potential begins (here: 10.0)
- radius where the second potential begins (here: 15.0)
- force constant of first potential (here: 10.0)
- force constant of second potential (here: 20.0)
- exponent of first potential (here: 2.0 - "normal" harmonic potential)
- exponent of second potential (here: 4.0 - more steep potential)
In order to avoid problems with velocities that are too high the radius of the inner potential should be chosen in a ways that none of the atoms in the starting structure feel the additional potential.
By the way... You might wonder why there is a spherical potential for task MD when there is also a spherical potential that is available in all tasks of CAST, i. e. also in MD simulations. This is a good question. I guess it has some historical reasons.
Another possibility to add an additional potential to an MD simulation is the option MD_biased_potential. There the potential is not centered around the geometrical center of the molecule but you can define an active site by choosing some atoms. The potential then depends on the distance to the active site:
- if the distance is smaller than the inner cutoff (first number in
MDcutoff) you have "normal" movement - if the distance is bigger than the outer cutoff (second number in
MDcutoff) you have no movement - between these extremes a gradual scaling of the velocities is applied
But let's see an example. For this we use the structure ubq.arc and the following CAST.txt:
verbosity 3
name ubq.arc
outname not_adjust
inputtype TINKER
task MD
interface OPLSAA
paramfile oplsaa.prm
cutoff 12.0
switchdist 10.0
MDsteps 5000
MDintegrator 0
MDtimestep 0.001
MDtrack 1
MDsnap 1000
MDheat 0 0.001
MDheat 500 300
MDveloscale 0 # switch off veloscale
#### OPTIONS FOR BIASED POTENTIAL ####
MDbiased_potential 1 # switch on and off
MDactive_site 1008 # define active site
MDactive_site 1017
MDactive_site 1018
MDactive_site 1019
MDcutoff 10 15 # inner and outer cutoff
MDadjust_by_step 0 # adjust active site every step
When we open the trajectory in VMD we see that only a part of the molecule moves, that part that is next to the NH3 group that defines the active site.
Attention! You should always switch off MDveloscale if you use this option!
There is one more option for the biased potential: MD_adjust_by_step. If this is set to 0 the position of the active site and the distances of every atom to it are only calculated once at the beginning of the simulation. If it's set to 1 the active site and the distances are updated every MD step. To see the effect you can perform another simulation where you switch on this option (rename your outputfile to adjust) and then compare the two trajectories. Especially take a look at the phenyle ring with atom indices 62 to 67 (VMD numbering). In the not_adjust trajectory it moves the whole time, in the adjust trajectory it only moves when the NH3 group that defines the active site points towards it.
By the way... If you only want to fix certain atoms without scaling down their velocity you might rather do this with the general options FIXrange (give atom indices of atoms you want to fix) or FIXsphere (fix all atoms outside of a sphere of a given atom).
With MDtrackoffset you can determine how many information is written to the outputfile "output_MD_TRACE.csv". The trajectory "output_MD_SNAP.arc" is not affected. To test this option just use the very first inputfile (or any other) and add the following line to it:
MDtrackoffset 100
When opening the file "output_MD_TRACE.csv" you see that only every 100th frame is written into it.
With the option MDsnap_buffer you decide how many snapshots are saved before they are written into a file. To test it add the following line to you inputfile:
MDsnap 50000 # make a snapshot at every step
MDsnap_buffer 100
The start CAST and break the run after some time. Look at how many snapshots are saved (by opening "output_MD_SNAP.arc" in VMD). There will always be multiples of 100 snapshots (i. e. 100, 200, 300, ...) because only when 100 snapshots are "full" they are written into file.
When switching on the option MDsnap_opt every snapshot is optimized by a local optimization. To test it use the following options:
verbosity 4 # needed to compare output
MDsnap 10 # only a small number of snapshots
MDsnap_opt 1 # optimize snapshots
Compare the console output at the end to the one you get if you just let the task LOCOPT run on your structure and you'll see that 10 local optimizations are performed, one for each snapshot.
If you switch on MDpre_optimize a local optimization is performed before starting the simulation. For seeing this it might be a good idea to use a very short simulation so change the following:
MDsteps 10 # short simulation
MDsnap_opt 0 # do not optimize snapshot (might be confusing)
MDpre_optimize 1 # perform preoptimization
Again you can see the local optimization in the console output.
Refinement means that the pairlist for non-bonded interactions is newly built. This might be necessary if you use a cutoff because then atompairs whose distance was longer than the cutoff might move in a way that they are now nearer to each other or the other way round. If you're not sure what you're doing just set it to the recommended value of 100 (default is that there is no refinement). In order to see the refinement (as alert Refining structure/nonbondeds in console output, not as any effect) use the following inputfile:
#### GENERAL INFORMATION ####
verbosity 4
name tridecaalanin.arc
outname output
inputtype TINKER
#### TASK ####
task MD
#### OPTIONS FOR ENERGY INTERFACE ###
interface OPLSAA
paramfile oplsaa.prm
cutoff 12.0
switchdist 10.0
#### OPTIONS FOR MD TASK ####
MDsteps 50000
MDintegrator 0
MDtimestep 0.001
MDveloscale 0
MDsnap 1000
MDrefine_offset 10 # refinement every 10 steps
With the rattle algorithm atom distances can be kept constant during the MD simulation. Mostly this algorithm is used to keep constant all bonds to any hydrogen atom. In order to do this add the following lines to the inputfile of the section before (you might set the MDrefine_offset to 100).
MDrattle 1 # keep all hydrogen bonds constant with rattle
MDrattle_use_paramfile 1 # get distances from parameterfile
With the first line you tell CAST to keep all bonds constant where one bonding partner is a hydrogen atom. Then CAST needs to know which is the distance to which the bonds are constraint. The second line tells CAST to take them from the bond parameter of the parameter file. When you keep all bonds with a hydrogen constant you have to use this option. So we run CAST and look at the trajectory with VMD. There we can see that the length of all bonds with a hydrogen atom remains the same during the whole simulation.
It is also possible only to constrain certain bond length. In the next example we want the bond distance of the C-terminal carbonyle bonds to remain constant. So we set MDrattle to 2 (which means: only constrain certain bonds) and add the following lines to define our constraint bonds:
MDrattlebond 125 126
MDrattlebond 125 133
In the trajectory of this simulation the bonds with hydrogen again change a bit whereas the bonds of the carboxyle group don't change at all. The distance to which they are constraint is still taken from the parameterfile.
But we don't have a take the distances from the parameterfile we can also define them by ourselves. In order to do that we set MDrattle_use_paramfile to 0 and add the following lines:
MDrattledist 1.25
MDrattledist 1.5
So we constrain our first rattlebond (125, 126) to 1.25 angstrom, the second one (125, 133) to 1.5 angstrom. There must always be the same number of line with the option MDrattledist as with MDrattlebond as we need a distance for every bond we want to restrain.If we run CAST again we see that the carboxyle group is asymmetric now.
When getting the distances not from parameterfile but giving them ourselves we can also restrain distances between atoms that are not connected by a bond. We can for example keep the C-teminal carboxyle C atom and the N-terminal N atom at a distance of 3 angstrom by the following lines:
MDrattlebond 1 125
MDrattledist 3
In some cases MD simulations break, often this is the case if velocities get too big. Normally then CAST gives a lot of zeros and/or NaNs and finally the program breaks.
If you activate the option MDrestart_if_broken CAST tests after every MD step if all bonds are still intact, i.e. if the bond length of every bond is reasonable. If it recognizes broken bonds the coordinates of all atoms are set to the original coordinates at the beginning of the simulation and new random velocities are assigned (according to the currently desired temperature, i.,e. if you've applied a gradient the gradient
won't start again but the temperature will continue to act how you would expect it if the simulation had run normally).
In order to see this in action use the following inputfile:
#### GENERAL INFORMATION ####
verbosity 3
name tridecaalanin.arc
outname output
inputtype TINKER
#### TASK ####
task MD
#### OPTIONS FOR ENERGY INTERFACE ###
interface OPLSAA
paramfile oplsaa.prm
cutoff 12.0
switchdist 10.0
#### OPTIONS FOR MD TASK ####
MDsteps 1000
MDintegrator 0
MDtimestep 0.001
MDveloscale 0
MDsnap 1000
MDrefine_offset 100
MDheat 0 300
MDheat 1000 30000 # use a very high temperature to cause breaking
MDrestart_if_broken 1
In the console output you see that there is a lot bond breaking and starting again (in my case the first time between 700 and 800 steps) but in the outputfile "output_MD_TRACE.csv" you see that the temperature keeps rising normally (well, as normal as we told it). If we look at the trajectory with VMD we can see the structure unfolding and then jumping back to the starting structure when CAST recognizes it's broken.
It is also possible to perform a NPT simulation by controlling the pressure via Behrendsen barostat. It is important to know that it is only possible to use pressure control if periodic boundaries are activated because the pressure control is done by modification of the periodic box size. You can see this in VMD if you plot the distance between two of the box atoms against the simulation timestep.
So our example for this is a water box, water.arc. This is our inputfile:
#### GENERAL INFORMATION ####
verbosity 3
name water.arc
outname output
inputtype TINKER
#### TASK ####
task MD
#### OPTIONS FOR ENERGY INTERFACE ###
interface OPLSAA
paramfile oplsaa.prm
cutoff 12.0
switchdist 10.0
#### OPTIONS FOR PERIODIC BOUNDARIES ####
Periodics 1 21.0 21.0 21.0
Periodicp 1
#### OPTIONS FOR MD TASK ####
MDsteps 1000
MDintegrator 0
MDtimestep 0.001
MDsnap 1000
MDrefine_offset 100
MDheat 0 300
MDpress 1
MDptarget 1.0
MDpcompress 0.000046
MDpdelay 2.0
With the option Periodics we activate periodic boundaries. The first number is switching them on (1) or off (0). By setting Periodicp to 1 we tell CAST to print dummy atoms into the outputfile that mark the corners of the box. So after running this simulation we can look at the trajectory with VMD and see a lot of water atoms moving inside a box that is marked by 8 dummy atoms. Sometimes we see a water molecule disappear at one end of the box and entering it at the other side. This is caused by the periodicity as we could imagine a lot of these boxes next to each other and it doesn't matter in which one of them a molecule is.
The last four line control the Behrendsen barostat. MDpress just switches it on. The MDptarget option gives the target pressure in atmosphere, so in our case it is 1 atm. You can see the current pressure in the file "press_MD_TRACE.txt" in the third column. The values are quite far from the target pressure and
variy quite strongly. However it is comparable to the results of a tinker calculation performed with similar parameters.
MDpcompress is the isothermal compressibility of the substance (the default value 0.000046 is for water) and MDpdelay barostat delay in picoseconds. It is needed to calculate the scaling factor for the box. When in doubt, just leave the default value of 2.0 ps.
Now you have seen some examples of the currently working options for task MD in CAST. Two options that have not been mentioned are MDrestart_offset and MDresume. These are for restarting a simulation that has broken for some reason at exactly the positions and velocities where it stopped. However saving and restarting a simulation does currently not work so I can't show an example for that. What is also not shown are the options for analyzing MD simulations MDana_pair, MDanalyze_zones and MDzonewidth. They are for plotting some data while doing an MD simulation. For details about them (and many of the other options please) see our manual (LaTeX file).
Date: 17 Jun 2019