Umbrella sampling
The tcsh shell script sample_us in samples/umbrella_sampling invokes five ChemShell runs to sample the free-energy barrier of the rotation of ethanol around its C–C bond:
#!/bin/tcsh -f
# prepare a directory for the sampling results
rm -rf umb_files
mkdir umb_files
# run the simulations
foreach bias ( 60 30 0 -30 -60 )
set var = `echo $bias\*3.1415927/180.| bc -l`
echo running window for $bias degrees
sed -e "s/VAR_bias/$bias/g" ethanol_umbrella.chm > onewindow.chm
chemsh onewindow.chm > log
grep "Restraint 1 tors : value" log | cut -c26-40 > torsion
sed -e "s/en/0.03 $var/g" torsion > umb_files/window$bias
cp torsion torsion$bias
end
# analyse
umbrella_integration.x -ss 500 -d umb_files -min -1.7 -max 1.7 -v
The system is defined in the file ethanol_umbrella.chm. In this example, a force field description is used. The restraint for umbrella sampling is included in the dynamics command by
dynamics dyn1 coords=c \
...
restraints = { { torsion 2 1 6 9 VAR_bias 0.03 } } \
where VAR_bias is replaced by the actual value by a sed script.
Five windows with the bias cantered at -60, -30, 0, 30, and 60 degrees, respectively, are simulated. In this example, only 3 ps are simulated in each window. To achieve proper sampling, one want to equilibrate much better and sample longer in real applications. One also might use more windows. Note that the value of the bias has to be input in radians.
After each window was run, the trajectory of the reaction coordinate is extracted into a file called torsion with the value of the bias appended.
The results are analysed with the tool umbrella_integration.x the source code of which can be found in the directory contrib. Try umbrella_integration.x -h for an explanation of the command-line arguments.
The detailed results (free-energy surfaces) are written to a file fe_ui.xy for the umbrella integration analysis, and a file fe_wham.xy for wham analysis. One can see that WHAM analysis produces more noise. The main results (energy differences between stationary points) of umbrella integration are printed to standard output. The difference between the minima (about 0.4 milli-Hartree, with an error bar of 0.7 mH) provides an estimate of the accuracy of the sampling. The barrier hight results in about 4.6 mH p/m 0.5 mH. This of course depends on the force field (and no effort was made here to use a realistic force field).
Plotting the global histogram of the sampling (the file global_histogram.xy) shows that, although the whole range of the reaction coordinate within the boundaries had been sampled, the system mainly resides in the minima.