QuickStartGuide-LAMMPS

Quick Start Guide: LAMMPS

Generate input configuration for solvation free energy calculation

Input configuration file for erdst and slvfe will be generated via helper scripts.
First, execute the following command to generate system information:

(ERmod directory)/share/ermod/tools/LAMMPS/gen_structure --data solution.data

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 argument for --data refers to the data file for the solution system containing the solute species.
The option --data can be simplified to -d, and takes a LAMMPS data file as an argument.
Upon execution, user is asked to state which molecule is the solute.
For example, user is asked

Molecule types in topology file: C2H6O H2O
Which molecules are solutes? (For multiple choice please specify as comma-separated list) 

When C2H6O is the solute, user types “C2H6O” (do not include quotation marks) to proceed.
After the execution, two directories are created, namely soln and refs.
The soln directory has MDinfo, MolPrmX (X = 1, 2, ...) , SltInfo, and .extraparam,
and the refs directory has MDinfo, MolPrmX (X = 1, 2, ...) , and SltInfo.
See Parameter files for erdst for detailed description of SltInfo, MolPrm1, and so on.

soln/MDinfo should read, for example, as

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 molecule of C2H6O type and 500 molecules of H2O type (C2H6O and H2O are the molecular types in the topology file).
The last line states that 9 atoms are present within a single C2H6O molecule and that 3 atoms are present within a single H2O molecule.

In the reference system (refs directory,) MDinfo file is also generated.
But the number of molecular species in refs/MDinfo is smaller by 1 than that in soln/MDinfo, since refs/MDinfo carries only the information of the solvent system without solute.

Both in soln and refs directories, MolPrmX (X = 1, 2, ...) and SltInfo should be found.
MolPrmX is the parameter file for solvent species, 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 erdst for detailed description of MDInfo, SltInfo, MolPrmX, and so on.

Now the parameters for the solution system are prepared.
Move to soln directory and run the helper script:

cd soln
(ERmod directory)/share/ermod/tools/LAMMPS/gen_input --traj ../solution.xtc --log ../solution.log --input ../solution.in

The trajectory file (.xtc), the log file of MD (.log), and the LAMMPS input script (.in) are fed to the script.
The options --traj, --log, and --input can be simplified to -t, -l, and -i, 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 erdst program runs.
MDinfo is further to be updated; “FRAME” on the first line is replaced by the number of snapshots in trajectory.
Between the .log and .in files, the information within the former is of higher priority.

Start the solute-solvent interaction calculation in the solution system with:

(ERmod directory)/bin/erdst

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.

cd ../refs
(ERmod directory)/share/ermod/tools/LAMMPS/gen_input --traj ../solvent.xtc --log ../solvent.log --input ../solvent.in --flexible ../solute.xtc
(ERmod directory)/bin/erdst

The first three arguments are to be taken similarly to the case of the solution (of course, the .xtc, .log, and .in files are changed to those of the solvent system without solute). The fourth 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.xtc, and SltConf is linked to the trajectory file of the isolated solute, solute.xtc.

After running erdst, 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.
The two-digit numbers indicate which fraction of the trajectory the file corresponds to.
For example, the trajectory is divided into 5 blocks.
engref.01 and corref.01 are then obtained from the first one of the 5 blocks, and engref.02 and corref.02 are from the second one.
Correspondingly, the output distribution functions are obtained from each of the 5 blocks.

When two erdst 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.

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 erdst program is run.
See Parameter files for erdst 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 in this case.
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:

cd ../refs
(ERmod directory)/share/ermod/tools/LAMMPS/gen_input --traj ../solvent.xtc --log ../solvent.log --input ../solvent.in --rigid ../solute.data
(ERmod directory)/bin/erdst

The difference from the flexible solute case is the fourth argument.
The fourth 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 erdst 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.

The option of --help or -h shows the help message for gen_structure and gen_input.

At the moment, gen_input cannot read parameters from variable in the LAMMPS input script for

fix npt
fix nvt
fix temp/*
pair_style 
dump
kspace_style

For the above, variable should not be used in an LAMMPS input script that is to be fed to gen_input.
For example,

variable        temp equal 300             # Final temperature (K)
variable        taut equal 1000            # Coupling time of thermostat (fs)
variable        taup equal 2000            # Coupling time of barostat (fs)
fix             1 all npt temp ${temp} ${temp} ${taut} iso 1 1 ${taup}

are valid in an LAMMPS input script, but are not for gen_input.
To read these info with gen_input, they need to be modified into

fix     1 all npt temp 300 300 1000 iso 1 1 2000