Currently there are no CHARMM / NAMD users in dev team, and script for CHARMM/NAMD is currently not very well maintained. If you found any strange / suspicious behavior please let us know!
Download CHARMM force field parameters from http://mackerell.umaryland.edu/CHARMM_ff_params.html, and copy both CHARMM22 force field file (top_all22_prot.inp) and stream parameters for ETOH (stream/toppar_all22_prot_model.str) to working directory. Run VMD (http://www.ks.uiuc.edu/Research/vmd/), and type the following into VMD’s Tk console (accessible via Extensions -> Tk console).
topology top_all22_prot.inp
topology toppar_all22_prot_model.str
resetpsf
segment ETOH { residue 1 ETOH }
coord ETOH 1 C1 { 0 0 0 }
coord ETOH 1 C2 { 1.52 0 0 }
guesscoord
writepsf etoh.psf
writepsf etoh.pdb
solvate etoh.psf etoh.pdb –o solution -t 11.838
The resulting structure will appear as solution.psf and solution.pdb, with 500 water molecules (1509 atoms in total). If the number of atoms does not match, try tweaking –t argument. Also, etoh.psf and etoh.pdb are produced, which can be used to specify the solute structure.
Similarly, create a water box for the reference solvent system. The system should contain exactly 500 molecules of waters. The resulting structure should appear as solvent.psf and solvent.pdb.
::::tcl
solvate -o solvent -minmax {{0 0 0} {0 0 0}} -t 12.888
To obtain the solvation free energy, user runs the solution system of interest, the pure solvent system, and the isolated solute in vacuum.
The simulation of isolated solute is not necessary if the solute is treated as rigid.
The solution system is first minimized and equilibrated in NVT and is then equilibrated in NPT.
The pure solvent system is minimized and equilibrated in NPT.
All bond lengths may be constrained for speeding up calculation.
For detailed instruction consult the NAMD manual.
Example configuration files are available on https://sourceforge.net/projects/ermod/files/data_example/.
MD production runs are performed in NPT condition, over 100 ps for the solution system to take 10000 snapshots with 10 fs interval, and over 50 ps for the pure solvent to take 500 snapshots with 100 fs interval. An MD of isolated solute is done over 50 ns in addition, to take 500000 snapshots with 100 fs interval.
In the following example, the trajectory files for solution system, pure solvent system, isolated solute system are named solution_run.dcd, solvent_run.dcd, and solute_run.dcd, respectively. The log files are named similarly.
Input configuration files for ermod and slvfe will be generated via a few helper scripts. First, execute the following command to generate an input file containing the system information:
::::sh
(ERmod directory)/share/ermod/tools/NAMD/gen_structure --psf solution.psf --log solution_run.log
where (ERmod directory) is the directory where the ERmod programs are installed; the directory was specified with the option of --prefix of ./configure described in General build guide.
The log file specified as the last argument can be any log file that uses solution.psf (is only used to find appropriate parameter files).
The options --psf and --log can be simplified to -t and -l, respectively.
path/to/ermod should be replaced with an appropriate path from the current working directory. Upon execution, user specifies the solute segment name.
In the example case, the ethanol molecule is the solute; type “ETOH” (do not include quotation marks) to proceed.
After the execution, two directories are created, namely soln and refs.
Inside these directories, soln/MDinfo
, soln/MolPrm1
, soln/SltInfo
, refs/MDinfo
, refs/MolPrm1
, and refs/SltInfo
files are generated.
See Parameter files for ermod for detailed description of SltInfo, MolPrm1, and so on.
Caution: When the number of solute molecules in the solution system of interest is more than unity, user needs to modify MDinfo, MolPrmX, and SltInfo manually.
See Finite concentration of solute: when more than one molecules are present for the solute species.
The manual modification of MDinfo, MolPrmX, and SltInfo in the case of more than one solute molecules is necessary only for NAMD (and CHARMM).
Caution: When ERmod is of ver 0.3.5 or earlier, the solution MD needs to be arranged so that the solute molecule appears as the first molecule in the coordinate files. This restriction is true only in the case of NAMD (and CHARMM). Further, even when the solution MD is not arranged so, one can run ermod correctly when MDinfo, MolPrmX, and SltInfo are modified manually. Still, this modification needs the knowledge about the structure of the input files for ermod, it is advised that the solute molecule be placed as the first in the coordinate files when the solution MD is conducted with NAMD (or CHARMM). Note that the above restriction is removed as of ver 0.3.6.
soln/MDinfo should read as
::::text
FRAME 2
1 500
9 3
The second column of the first line indicates that there are 2 types of molecular species.
The second line states that the system contains 1 ETOH molecule and 500 water molecules.
The last line states that 9 atoms are present within a single ethanol molecule and that 3 atoms are present within a single water molecule.
In the reference system (refs directory,) MDinfo file is also generated.
But the number of molecular species in refs/MDinfo is 1, since MDinfo in this case carries the information of the pure solvent system.
Both in the soln and refs directories, MolPrm1 and SltInfo should be found.
MolPrm1 is the parameter file for solvent molecule, and SltInfo is the parameter file for solute. Both files contain the Lennard-Jones parameters and charge information, which are necessary in the calculation. See Parameter files for ermod for detailed description of MDInfo, SltInfo, MolPrm1, and so on.
Now the parameters for the solution system are prepared.
Move to soln directory and run the helper script:
::::sh
cd soln
(ERmod directory)/share/ermod/tools/NAMD/gen_input --dcd ../solution_run.dcd --log ../solution_run.log
The trajectory (dcd file) and log files are fed to the script.
The options --dcd and --log can be simplified to -x and -l, respectively.
After running the script, a symbolic link file, namely HISTORY, is to be created.
A file named parameters_er is also to be created; this file controls how the ermod program runs.
MDinfo is further to be updated; “FRAME” on the first line is replaced by the number of snapshots in trajectory.
Start the solute-solvent interaction calculation in the solution system with:
::::sh
(ERmod directory)/bin/ermod
If the program runs successfully, there should be files named engsln.xx and aveuv.tt, where xx is a two-digit integer, besides several other outputs.
The former is the histogram of the solute-solvent pair interaction energy, while the latter is the averaged sum of the solute-solvent interaction energy.
The two-digit numbers indicate which fraction of the trajectory the file corresponds to.
For example, the trajectory is divided into 10 blocks.
engsln.01 is then the histogram obtained from the first one of the 10 blocks, and engsln.02 is from the second one of the 10 blocks.
Correspondingly, the output distribution functions are obtained from each of the 10 blocks.
The division into blocks is done to examine the convergence of the solvation free energy, as described in the end of Quick Start Guide ---- Running slvfe: Getting Final output
Then move on to refs directory to start calculating the reference part.
::::sh
cd ../refs
(ERmod directory)/share/ermod/tools/NAMD/gen_input --dcd ../solvent_run.dcd --log ../solvent_run.log --flexible ../solute_run.dcd
(ERmod directory)/bin/ermod
The first two arguments are to be taken similarly to the case of the solution (of course, the dcd and log files are changed to those of pure solvent).
The third argument, --flexible, specifies that the test-particle insertion calculation is performed with flexible structures of the solute.
The option --flexible can be simplified to -s.
In addition to a symbolic link named HISTORY, a symbolic link file named SltConf is generated; HISTORY refers to the trajectory of the solvent, solvent_run.dcd, and SltConf is linked to the trajectory file of the isolated solute, solute_run.dcd.
After running ermod, user obtains the energy distributions for the insertion operation.
engref.xx is the histogram of the solute-solvent pair interaction energy at insertions in reference system, and corref.xx is the correlation matrix.
User may notice that there are 5 outputs of energy distribution functions.
This is because the reference trajectory is divided into 5 sections,
and the output distribution functions are obtained from each of the 5 blocks.
When two ermod runs are finished in the soln and refs directories, the solvation free energy can be calculated from the distribution functions obtained.
Go to Quick Start Guide ---- Running slvfe: Getting Final Output.
A sample output is prepared as slvfe-samle.log.
The trajectory is divided into a set of blocks both for soln and refs.
The default value of the number of divisions is 10 for the solution system and is 5 for the reference.
The number of divisions can be modified with the numdiv parameter in the parameters_er file when the ermod program is run.
SeeParameter files for ermod for the description of the parameters which can be specified in parameters_er.
If the solute is rigid (its intramolecular structure is fixed), the procedure is simplified.
The scheme for the solution part is the same as above.
Note, of course, that the solution MD is supposed to be done using a rigid solute.
Further, it is not necessary to perform a simulation of isolated solute when the solute is rigid.
The use of rigid solute simplifies the operation of the reference part.
In the refs directory, run the helper script:
::::sh
cd ../refs
(ERmod directory)/share/ermod/tools/NAMD/gen_input --dcd ../solvent_run.dcd --log ../solvent_run.log --rigid ../etoh.pdb
(ERmod directory)/bin/ermod
The difference from the flexible solute case is the third argument.
The third argument, --rigid, specifies that the test-particle insertion calculation is performed with fixed structure of the solute.
The option --rigid can be simplified to -r.
The structure of the solute is taken from the SltInfo file when the solute is rigid.
See Parameter files for ermod for the format of the SltInfo file in the case of rigid solute.
When the refs calculation is done, the calculation of solvation free energy of rigid solute can be done with the same procedure as that in the flexible case.
Wiki: FixedLigandBind
Wiki: MoreThanOneSolute
Wiki: QuickStartGuide
Wiki: build-Guide