Constraints in methanol
This is an exercise to see the effect of constraints and timesteps on simulations.
This exercise uses a system of 216 methanol molecules and an system of 1 methanol molecule in a cubic box with edge 2.4 nm. The methanol model using a unit-atom CH3 group, mainly to speed up the simulations for the exercise, but also to avoid the use of virtual interaction sites when we want to remove the angle vibrations involving hydrogens. The latter can now be achieved by constraining the angles involving hydrogens, as this does not affect other degrees of freedom (there are none left in this molecule in that case). The setup is an NVT ensemble with a thermostat for the 216 molecule system and a stochastic dynamics integrator for the single molecule to ensure ergodicity.
The input can be found here.
The idea is to run a combinations of settings:
time steps 0.5, 1, 2 and 5 fs
constraints: none, H-bonds, all-bonds and H-angles
Some settings will give, rightly so, a warning in grompp. You can force grompp to produce the run input file anyhow with the option: -maxwarn 1. Some combinations might lead to instability.
Will we calculate the heat of vaporization to compare the quality of the simulations. The heat of vaporatization is the net interaction energy between the molecules. This is calculated as the potential energy of a single molecule in the gas phase minus the average potential energy per molecule in the liquid. To obtain this you need to run for each setting a single molecule and a liquid simulation and then get the potential energy, which is most easily obtained by doing:
echo pot | g_energy
You can think about the results you obtained, but you might not be able to explain all differences you observe.
You can also look at the energy conservation by plotting the conserved energy.