The following tutorial explains how to compute potentials of mean force (PMF) or the free-energy profiles using HFB method as implemented in CHARMM as the driver and NAMD as the MD engine. We will compute the PMF for the [NaCl] solution of 1.0 M using the Cartesian gradient-augmented HFB method so that we do not need to worry about the Jacobian correction ever. We use this particular example because we have performed extensive benchmarks using two independent conventional techniques for this case. All the benchmarks have been performed using CHARMM for both the driver and the MD engine. As in the original publication we will be using NPT ensemble at T=298K and p=1atm with truncated octahedron box of 40Angstroms that contains 31 NaCl ion pairs and 1668 water molecules. Nonbonded interactions will be described using the following NAMD parameters (these are rather expensive, and perhaps unnecessarily so, so feel free to make any compromises you are used to).
wrapAll on
wrapNearest on
PME on
PMEGridsizeX 40
PMEGridsizeY 40
PMEGridsizeZ 40
margin 5
switching on
switchDist 10.0
cutoff 12.0
pairListDist 16.0
nonbondedFreq 1
fullElectFrequency 2
stepspercycle 20
For the integration of MD steps we will use
timestep 2
rigidBonds all
rigidTolerance 0.00000001
Although the validity of the Na+ and Cl- parameters as well as the CHARMM version of the TIP3P water (that has non-zero vdW parameters on its hydrogen atoms and has been estimated to have dielectric constant epsilon of 92+-5) could be debated, these parameters have been used by us to establish what we believe is high quality free-energy benchmarks. We also use the CHARMM benchmarks to establish the validity of the developed hybrid CHARMM/NAMD methodology (as described in details below) that employs NAMD code instead of CHARMM for the MD engine. In this respect we greately appreciate the NAMD development for allowing us to use CHARMM style parameter (PRM) and toplogy (PSF) files.
Because we use NAMD as the MD engine, we must be able to impose positional restraints on the RCS atoms using NAMD itself. In this particular case, the RCS space consists of one Na+ and one Cl- ion, leaving 1668 water molecules and 30 NaCl ions pairs in the SCS space (refer to the HFB documentation for definitions of the RCS and SCS spaces). Since we have two ions in the RCS space we will be restraining 6 Cartesian degrees of freedom, namely (x1,y1,z1) and (x2,y2,z2) for Na+ and Cl-, respectively. The off-the-shelf NAMD code provides two options to restrain the absolute positions of atoms to their reference positions in Cartesian space, namely, the "PDB-style" restraints and the "TCL-style" urestraint function with "posi" feature. For the purpose of the present tutorial we will use the former, as the most simple in terms of the NAMD/CHARMM hybridization.
We wish to use a truncated octahedron box with NAMD, just like we did with CHARMM. Note the difference between representations of the truncated octahedron in NAMD and CHARMM. NAMD uses a very simple and general triclinc representation of the unit cell, whereas CHARMM uses yet another equivalent representation. Although otherwise equivalent the two boxes do not look alike when visualized with VMD - so don't be surprised. To proceed with our tutorial we need the NAMD-style truncated octahedron box. If you are lazy, you can use the CHARMM-style box to start with and equilibrate it using NAMD (make sure to setup a proper triclinic representation of the truncated octahedron in the NAMD input file). It takes only 100 ps or so to morph one representation into the other by running an NPT MD simulation. If you want to be smarter than that, consider writing a coordinate transformation between the two representations.
The next thing we need to setup is the reference beads. We will use 32 beads to span distance range from 1.67 to 16 Angstroms and to ensure proper integration of the PMF. We will use CHARMM program to setup the reference beads. We choose to separate the ions along the x-axis from the reactant configuration Na+(1.01212,0.0,0.0); Cl-(-0.65638,0.0,0.0) to the product configuration Na+(9.70565,0.0,0.0); Cl-(-6.29423,0.0,0.0) using linear interpolation within TReK/HFB module (note that center of mass of the ion pair is at the origin).
Let us setup a psf file for the ion pair, and prepare two charmm coordinate files (chr) with the specified reactant and product structures. Now we use HFB to generate all the missing beads between the reactant and the product at equal intervals. The following charmm input script takes care of all these steps(provide a hyperlink to the file here).
At this point we should have 32 reference structures in charmm format. However these only include a single ion pair and do not have the other ions or water molecules. Using a psf and a chr files for the full system, the following charmm input script prepares complete pdb-style reference beads to be used with NAMD(provide a hyperlink to the file here).
Next step is to equilibrate all the beads to the obtained reference structures using NAMD. The corresponding NAMD input script is provided here(provide a hyperlink to the file here). To avoid overwhelming file editing, we use the following perl-script that takes care of all the tasks mentioned above(provide a hyperlink to the file here).
(Sorry the tutorial is currently under construction. More information will become availabel as we go.)