[gmx-users] Automated simulations with Gromacs: MD, BAR, Umbrella, Force Field conversion.
Andrey Frolov
andrey.i.frolov at mail.ru
Tue Jul 8 00:18:38 CEST 2014
Dear GMX-users,
I would like to share some scripts which help to automate MD and free energy
calculations with Gromacs. They are written in Bash and thus can be easily
customized for one's specific needs. At present everything is designed to run
absolute solvation free energy calculations of small organic molecules/ligands in
arbitrary (mixed) solvent. The following scripts are available:
(http://frolov-pchem.wikispaces.com/Downloads)
1) Pipeline for molecular dynamics simulations with Gromacs: run_md.sh.
2) Pipeline for BAR absolute solvation free energy calculations: bennett_pipe.sh
3) Pipeline for 1 collective variable umbrella sampling for Gromacs + PLUMED:
umbrella_1CV.sh
4) Force field converter: ffconv.py
5) Python unit converter: unitconv.py
All scripts come with real life examples how to use them
(http://frolov-pchem.wikispaces.com/Tutorials), where you can find a one-click
RUN_ALL.sh script, which performs all the computations in one go. If you
have difficulties, please, let me know, I can assist.
The pipeline scripts are intended to save user's time on writing a code for
automation of calculations. Below a more detailed description of the scripts
is given.
--------------------------------------------------------------------------------
1. Pipeline for molecular dynamics simulations with Gromacs: run_md.sh.
--------------------------------------------------------------------------------
This script can perform in one go or separately the following steps: building
the topology and coordinates of the system, performing energy minimization, NVT
and NPT equilibration, production run MD:
EXAMPLE USAGE.
paracetamol dissolved in supercritical CO2 with addition of ethanol at 250 bar and 308.15K
run_md.sh \
--lbuildcoor yes \
--lbuildtop yes \
--FF oplsaa \
--box 3.5,3.5,3.5 \
--components paracetamol,ethanol,CO2 \
-N 1,32,485 \
--toppath $p/top \
--mdppath $p/mdp \
--mdseq em,nvteq,npteq,run \
--mdpkeys ref_t=308.15,gen_temp=308.15,ref_p=250.0
--emmdpkeys nsteps=200 \
--nvteqmdpkeys nsteps=5000 \
--npteqmdpkeys nsteps=5000,nstxtcout=1000 \
--runmdpkeys nsteps=10000,nstxtcout=100,nstenergy=1000
see: http://frolov-pchem.wikispaces.com/Tutorials#MD
--------------------------------------------------------------------------------
2. Pipeline for BAR absolute solvation free energy calculations: bennett_pipe.sh
--------------------------------------------------------------------------------
The script is an environment to automatize free energy calculations with
Gromacs (www.gromacs.org). This script uses "run_md.sh" for running MD simulations,
therefore it also inherits all the features of that script. It was designed
to calculate solvation free energies and herefore operates with lambda_coul
and lambda_vdw scaling parameters (see Gromacs manual). For each pair of lambda
parameters it prepares the simulation setups. Then one can submit for each pair
of lambda values equilibration and/or production simulation to a queue system.
bennett_pipe.sh --couplemoltype paracetamol \
--lcouplegro yes \
--vdwlambdas 0.0__0.0__0.0__0.2__0.3__0.4__0.5__0.6__0.7__0.8__1.0 \
--coullambdas 0.0__0.5__1.0__1.0__1.0__1.0__1.0__1.0__1.0__1.0__1.0 \
--grofile `pwd`/initdata/system.gro \
--runmdkeys "--mdseq em,nvteq,npteq,run --lbuildtop no --lbuildcoor no --lbuildmdp yes --topfile `pwd`/initdata/system.top --mdppath `pwd`/mdp --nvteqmdpkeys nsteps=5000 --npteqmdpkeys nsteps=5000 --runmdpkeys nsteps=5000 --mdruncmd mdrun_d"
See: http://frolov-pchem.wikispaces.com/Tutorials#BAR
--------------------------------------------------------------------------------
3. Pipeline for 1 collective variable umbrella sampling for Gromacs + PLUMED
--------------------------------------------------------------------------------
The script is an environment to calculate potential of mean force (PMF) for a
single collection variable (CV) using PLUMED (www.plumed-code.org) modified
Gromacs package (www.gromacs.org).
umbrella_1CV.sh --CV 20:20:180 \
--CVunits degree \
--ks 300.0 \
--ksunits kJ/mol/rad^2 \
--leql yes \
--lrun yes
See: http://frolov-pchem.wikispaces.com/Tutorials#Umbrella
--------------------------------------------------------------------------------
4. Force field converter: ffconv.py
--------------------------------------------------------------------------------
ffconv.py is a python program which performs conversion between different file
formats of classical force field for small molecules. Supplementary scripts check
if the conversion was correct by comparing energies calculated before and after
conversion.
Now ffconv can do the following:
to read:
1) OPLS-AA: Schrodinger ffld_server output
2) AMBER topology file CHARMM force field file
3) CHARMM molecule stream file
4) PDB coordinate file
to write:
1) Gromacs topology itp files
2) RISM-MOL topology file
3) Amber topology file - only for single point calculations (e.g. with
3D-RISM)
to test conversion to/from:
1) ffld_server (runs "ffld_server"),
2) AMBER (runs NAMD... yes, NAMD because AMBER is not free and AmberTools
utilities (e.g. NAB) are not able to understand all the fields
of Amber topology file)
3) CHARMM (runs charmm)
to store the topology files into local topology database, which is compatible with
"run_md.sh" pipeline system building capabilities.
EXAMPLE USAGE (convertion of ffld_server generated topology to Gromacs format)
### to GROMACS
# 1) convert
$FFCONVPATH/ffconv.py --fromprog FFLDSERVER \
--toprog GROMACS \
--ffname oplsaa \
--moltopfn paracetamol.ffld \
--ffconvpath $FFCONVPATH
# 2) check conversion
$FFCONVPATH/scripts/check_conversion.sh \
--mol paracetamol \
--fromprog FFLDSERVER \
--toprog GROMACS \
--confpath `pwd` \
--ldbstore1 yes \
--ldbstore2 yes \
--dbpath $FFCONVPATH/top \
--ffname oplsaa
STRUCTURE
The code is an objective oriented python code. First, the input topology files
are read and stored in a molecule and/or ForceField objects. The information in
these objects is used to write out the topology files. Therefore, to incorporate
a new force field format one has to write a function which fills in the information
to the objects and/or function which writes the objects information in the required
format.
molecule_class.py - contains the classes definition and internal functions.
func_*_io.py - I/O functions
So-far, unique functionality of ffconv.py (with respect to "amb2gmx.pl" and "charmm_to_gromacs") is:
1) the ability to have automated OPLS-AA ligand topology creation for Gromacs
2) checking if the conversion was correct in each particular case
3) ability to be easily extended to other topology file formats
4) topology conversion into RISM-MOL and AMBER formats for single molecules RISM and 3D-RISM calculations.
see: http://frolov-pchem.wikispaces.com/Downloads#ffconv.py
--------------------------------------------------------------------------------
5. Python unit converter: unitconv.py
--------------------------------------------------------------------------------
The script takes two command line arguments and returns the factor to
convert units.
EXAMPLE USAGES:
unitconv.py kJ/mol/rad^2 kcal/mol/rad^2
unitconv.py mol*m^-3 nm^-3
unitconv.py kcal/mol toGMX
Also, one can load the module in his/her own python code.
see: http://frolov-pchem.wikispaces.com/Downloads#unitconv.py
Please, give your feedback to: andrey.i.frolov at mail.ru
Please note, that it was observed that wikispaces.com does not always work well with Mozilla Firefox (does not load pages, crashes).
With recent Google Chrome no problems were observed.
Hope this will be of help.
With kind regards,
Andrey
(Former position)
Junior researcher
Institute of Solution Chemistry
Russian Academy of Sciences
Akademicheskaya St. 1
153045 Ivanovo, Russia
(Current position)
Sanofi-aventis R&D,
Lead Generation to Compound Realisation,
Structure Design & Informatics,
94403 Vitry, France
More information about the gromacs.org_gmx-users
mailing list