[gmx-users] Comparing Charmm 27 energies from GROMACS and NAMD
Reza
resal81 at gmail.com
Mon May 13 21:58:43 CEST 2013
Here is the script:
#!/bin/bash
# note: you must have "par_all27_prot_lipid.prm" in the starting directory
set -e
PROGRAMS="/g1/home/resal/Programs"
GMX_BIN_457="$PROGRAMS/gromacs/4.5.7/thread/bin"
GMX_BIN_461="$PROGRAMS/gromacs/4.6.1/thread/bin"
NMD_BIN="$PROGRAMS/namd/NAMD_2.9_Linux-x86_64-multicore"
VMD_BIN="$PROGRAMS/vmd/1.9.1/bin"
if [ -d prep ]; then
rm -rf prep
fi
if [ -d gromacs-4.5.7 ]; then
rm -rf gromacs-4.5.7
fi
if [ -d gromacs-4.6.1 ]; then
rm -rf gromacs-4.6.1
fi
if [ -d namd-2.9 ]; then
rm -rf namd-2.9
fi
mkdir prep
cd prep
cat > psfgen << EOF
# create the membrane
package require membrane
membrane -l POPC -x 40 -y 40 -o membrane
mol delete all
# delete water molecules
mol load psf membrane.psf pdb membrane.pdb
package require psfgen
readpsf membrane.psf
coordpdb membrane.pdb
set wat_sel [atomselect top "water"]
set wat_segs [lsort -unique [\$wat_sel get segid]]
foreach wseg \$wat_segs {
set sel [atomselect top "segid \$wseg and water"]
set wat_res [lsort -unique [\$sel get resid]]
foreach res \$wat_res {
delatom \$wseg \$res
}
}
writepsf memb_nowat.psf
writepdb memb_nowat.pdb
quit
EOF
$VMD_BIN/vmd -dispdev none -e psfgen
cd ..
mkdir gromacs-4.5.7
cd gromacs-4.5.7
cat > pdb2gmx.in << EOF
8
6
EOF
cat > mdpfile.mdp << EOF
integrator = md
nsteps = 0
nstlog = 1
nstlist = 0
ns_type = simple
rlist = 0
coulombtype = cut-off
rcoulomb = 0
rvdw = 0
pbc = no
EOF
GMXBIN=$GMX_BIN_457
source $GMXBIN/GMXRC
$GMXBIN/pdb2gmx -f ../prep/memb_nowat.pdb # << pdb2gmx.in # choose CHARMM27 for ff, None for water
$GMXBIN/grompp -f mdpfile.mdp -p topol.top -c ../prep/memb_nowat.pdb -o topol.tpr
$GMXBIN/mdrun -nt 1 -s topol.tpr -g gromacs.log
cd ..
mkdir gromacs-4.6.1
cd gromacs-4.6.1
GMXBIN=$GMX_BIN_461
source $GMXBIN/GMXRC
$GMXBIN/pdb2gmx -f ../prep/memb_nowat.pdb # << ../gromacs-4.5.7/pdb2gmx.in
$GMXBIN/grompp -f ../gromacs-4.5.7/mdpfile.mdp -p topol.top -c ../prep/memb_nowat.pdb -o topol.tpr
$GMXBIN/mdrun -nt 1 -s topol.tpr -g gromacs.log
cd ..
mkdir namd-2.9
cd namd-2.9
cat > conf << EOF
structure ../prep/memb_nowat.psf
coordinates ../prep/memb_nowat.pdb
paratypecharmm on
parameters ../par_all27_prot_lipid.prm
exclude scaled1-4
1-4scaling 1.0
switching off
cutoff 1000
pairlistdist 1000
timestep 1.0
outputenergies 1
outputtiming 1
binaryoutput no
outputname namd
dcdfreq 1
temperature 300
run 0
EOF
$NMD_BIN/namd2 conf > namd.out
cd ..
cat > analyze_energies.py << EOF
import sys
lines = open(sys.argv[1]).readlines()
if sys.argv[1].endswith('log'):
match = lambda line, fields: all([f in line for f in fields])
for i, line in enumerate(lines):
if match(line, ('Bond', 'U-B', 'LJ-14')):
B, A, D, I, L1 = map(float, lines[i+1].split())
elif match(line, ('Coulomb-14','LJ (SR)')):
C1,L2, C2,P, K = map(float, lines[i+1].split())
elif match(line, ('Total Energy', 'Temperature', 'Pressure (bar)')):
break
print sys.argv[1]
print '%10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f' % (
B/4.184, A/4.184, D/4.184, I/4.184,
(C1+C2)/4.184, (L1+L2)/4.184, P/4.184)
elif sys.argv[1].endswith('out'):
for line in lines:
if line.startswith('ENERGY:'):
f = line.split()
if f[1] == '0':
E = map(float, f[2:])
B, A, D, I, C, L = E[0:6]
P = E[11]
break
print sys.argv[1]
print '%10.3f %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f' % (
B, A, D, I, C, L, P)
EOF
echo
echo "==============================="
echo " BOND ANGLE DIH IMP COUL LJ POT"
python analyze_energies.py gromacs-4.5.7/gromacs.log
python analyze_energies.py gromacs-4.6.1/gromacs.log
python analyze_energies.py namd-2.9/namd.out
echo "==============================="
echo
On May 13, 2013, at 3:57 PM, Reza <resal81 at gmail.com> wrote:
> So I think I figured out what was causing the discrepancy of Charmm27 energies between gromacs and NAMD. It appears that it's related to the gromacs version: energies from 4.5.7 and NAMD match very well while 4.6.1 gives energies that are different from both 4.5.7 and NAMD.
>
> Here are the results form my tests on POPC membrane:
>
> ===============================
> BOND ANGLE DIH IMP COUL LJ POT
> gromacs-4.5.7/gromacs.log
> 1539.120 3111.902 1250.540 16.284 -1705.306 -1219.345 2993.188
>
> gromacs-4.6.1/gromacs.log
> 1539.120 3111.902 1250.540 16.284 -5997.992 350.084 269.945
>
> namd-2.9/namd.out
> 1539.116 3111.920 1250.542 16.284 -1705.307 -1219.343 2993.211
> ===============================
>
>
> I think this could be due to either 1) the simulation setup should be different between 4.5.7 and 4.6.1 when running single point calculations, or 2) my 4.6.1 compilation is not correct (although all regression tests passed), or 3) there is a bug somewhere.
>
> I'll attach a script that automates the comparison between GROMACS and NAMD in a separate post (it will make this post > 50kB). If you run it and got similar/different results as mine, I'll be happy to know. Also any comment on the simulations setup is highly appreciated.
>
> Reza
>
>
>
> On May 3, 2013, at 11:57 AM, Reza <resal81 at gmail.com> wrote:
>
>>
>> Actually pdb which is similarly doesn't have velocity information. In this case however I'm mainly interested to see how the potential terms compare between the two packages with Charmm27 ff.
>>
>> Reza
>>
>>
>>>
>>>
>>> On 5/3/13 11:33 AM, Reza wrote:
>>>> Thanks Mark!
>>>>
>>>>> No, they weren't. See
>>>>> http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy
>>>>> You cannot hope to reproduce accurately the energy of a configuration if
>>>>> you let the coordinates be manipulated.
>>>>>
>>>>
>>>>
>>>> The "rerun" option is interesting - in my case however the potential energy terms stayed identical but the kinetic term became zero.
>>>>
>>>
>>> If you are using an .xtc file for the rerun, this makes sense since the .xtc does not store velocities and hence kinetic energy cannot be calculated.
>>>
>>> -Justin
>>>
>>>>>
>>>>> Something you think is equivalent is not :-) Move to testing a system with
>>>>> two lipids. Inspect all the logfile outputs very carefully for clues.
>>>>
>>>>
>>>> I totally agree :) So far I found out that for "no cut-off" simulation in Gromacs, rather that specifying a large cut-off, it needs rlist=rvdw=rcoulomb=0 and pbc=no along with ns_type=simple and nstlist=0 (according to the manual). I am running various tests and will update if I find out what is causing the discrepancy.
>>>>
>>>> Reza
>>>>
>>>> On May 1, 2013, at 5:46 PM, Mark Abraham <mark.j.abraham at gmail.com> wrote:
>>>>
>>>>> On Wed, May 1, 2013 at 2:32 PM, Reza Salari <resal81 at gmail.com> wrote:
>>>>>
>>>>>> Hi Justin,
>>>>>>
>>>>>> I actually did :) but it ended up being bigger than 50 kb so it needed
>>>>>> moderator approval to show up. I was hoping it would've been released by
>>>>>> now. I'll attach a the details below.
>>>>>>
>>>>>> Any help/hint is highly appreciated.
>>>>>>
>>>>>> Reza
>>>>>>
>>>>>> Details:
>>>>>>
>>>>>> *1)* Both systems were prepared using VMD "membrane" package and then
>>>>>> waters were removed. I used Gromacs 4.6.1 and NAMD 2.9.
>>>>>>
>>>>>> *2)* Simulations were run in vacuum as a single-point energy calculations
>>>>>> (0 step).
>>>>>
>>>>>
>>>>> No, they weren't. See
>>>>> http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy
>>>>> You cannot hope to reproduce accurately the energy of a configuration if
>>>>> you let the coordinates be manipulated.
>>>>>
>>>>> PME was not used.
>>>>>>
>>>>>> *3) *For Gromacs runs, pdb2gmx was used to prepare the system and the
>>>>>> output was saved as the pdb format. The mdp file:
>>>>>>
>>>>>> integrator = md
>>>>>> nsteps = 0
>>>>>> nstlog = 1
>>>>>> nstlist = 1
>>>>>> ns_type = grid
>>>>>> rlist = 100.0
>>>>>> coulombtype = cut-off
>>>>>> rcoulomb = 100.0
>>>>>> rvdw = 100.0
>>>>>> pbc = no
>>>>>>
>>>>>> *4) *NAMD input file:
>>>>>>
>>>>>> structure ../0_prep/memb_nowat.psf
>>>>>> paratypecharmm on
>>>>>> parameters par_all27_prot_lipid.prm
>>>>>> exclude scaled1-4
>>>>>> 1-4scaling 1.0
>>>>>> switching off
>>>>>> switchdist 8
>>>>>> cutoff 1000
>>>>>> pairlistdist 1000
>>>>>> margin 0
>>>>>> timestep 1.0
>>>>>> outputenergies 1
>>>>>> outputtiming 1
>>>>>> binaryoutput no
>>>>>> coordinates ../0_prep/memb_nowat.pdb
>>>>>> outputname out
>>>>>> dcdfreq 10
>>>>>> temperature 300
>>>>>> run 0
>>>>>>
>>>>>> *5)* Energies:
>>>>>>
>>>>>> For Single POPC (kcal/mol)
>>>>>>
>>>>>>
>>>>>> Gromacs NAMD Diff
>>>>>>
>>>>>> Bond 43.0022 43.0015 -0.0007
>>>>>> Angle 80.6568 80.6571 0.0003
>>>>>> Dih 29.8083 29.8083 0.0000
>>>>>> Imp 0.8452 0.8452 0.0000
>>>>>>
>>>>>> Coul -17.2983 -17.2983 0.0000
>>>>>> LJ -7.0798 -7.0798 0.0000
>>>>>>
>>>>>> Pot 129.9343 129.9340 -0.0003
>>>>>>
>>>>>>
>>>>> The intra-molecule terms look fine. Since this is a lipid, there are
>>>>> non-bonded interactions that are intra-molecule, so the non-bondeds also
>>>>> seem fine.
>>>>>
>>>>> For POPC Memb (kcal/mol)
>>>>>>
>>>>>> Gromacs NAMD Diff
>>>>>>
>>>>>> Bond 1539.1181 1539.1162 -0.0019
>>>>>> Angle 3111.9264 3111.9197 -0.0067
>>>>>> Dih 1250.5425 1250.5421 -0.0004
>>>>>> Imp 16.2837 16.2837 0.0000
>>>>>>
>>>>>> Coul -1837.8585 -1705.3075 132.5510
>>>>>> LJ -995.0311 -1219.3432 -224.3121
>>>>>>
>>>>>> Pot 3084.9904 2993.2110 -91.7794
>>>>>>
>>>>>
>>>>> Something you think is equivalent is not :-) Move to testing a system with
>>>>> two lipids. Inspect all the logfile outputs very carefully for clues.
>>>>>
>>>>> Mark
>>>>>
>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Tue, Apr 30, 2013 at 8:33 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On 4/30/13 4:19 PM, Reza Salari wrote:
>>>>>>>
>>>>>>>> Hi
>>>>>>>>
>>>>>>>> I have set up two small systems, one with a single POPC lipid, and
>>>>>> another
>>>>>>>> system with 23 POPC's arranged as a lipid bilayer. I am hoping to do a
>>>>>>>> similar comparison as in Table 1 in the Par Bjelkmar et al paper
>>>>>> (porting
>>>>>>>> charmm ff to gromacs) for my systems. My main question is that for the
>>>>>>>> single POPC, all the potential energy terms match very well, but for the
>>>>>>>> membrane system the non-bonding terms differ significantly.
>>>>>>>>
>>>>>>>> I am providing the full details below and greatly appreciate any hint
>>>>>> for
>>>>>>>> better comparison of the energies.
>>>>>>>>
>>>>>>>>
>>>>>>>> Thanks,
>>>>>>>> Reza Salari
>>>>>>>>
>>>>>>>> Details:
>>>>>>>>
>>>>>>>> 1) Both systems were prepared using VMD "membrane" package:
>>>>>>>>
>>>>>>>> 2)
>>>>>>>>
>>>>>>>>
>>>>>>> It appears you hit "send" too early. Please provide the entirety of the
>>>>>>> details you intended. Complete .mdp files and actual quantification of
>>>>>> the
>>>>>>> differences you are observing are also very important.
>>>>>>>
>>>>>>> -Justin
>>>>>>>
>>>>>>> --
>>>>>>> ==============================**==========
>>>>>>>
>>>>>>> Justin A. Lemkul, Ph.D.
>>>>>>> Research Scientist
>>>>>>> Department of Biochemistry
>>>>>>> Virginia Tech
>>>>>>> Blacksburg, VA
>>>>>>> jalemkul[at]vt.edu | (540) 231-9080
>>>>>>> http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin<
>>>>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin>
>>>>>>>
>>>>>>> ==============================**==========
>>>>>>> --
>>>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>>>> http://lists.gromacs.org/**mailman/listinfo/gmx-users<
>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users>
>>>>>>> * Please search the archive at http://www.gromacs.org/**
>>>>>>> Support/Mailing_Lists/Search<
>>>>>> http://www.gromacs.org/Support/Mailing_Lists/Search>before posting!
>>>>>>> * Please don't post (un)subscribe requests to the list. Use the www
>>>>>>> interface or send it to gmx-users-request at gromacs.org.
>>>>>>> * Can't post? Read http://www.gromacs.org/**Support/Mailing_Lists<
>>>>>> http://www.gromacs.org/Support/Mailing_Lists>
>>>>>>>
>>>>>> --
>>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>>> * Please search the archive at
>>>>>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>>>> * Please don't post (un)subscribe requests to the list. Use the
>>>>>> www interface or send it to gmx-users-request at gromacs.org.
>>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>>>
>>>>> --
>>>>> gmx-users mailing list gmx-users at gromacs.org
>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>>>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>>>> * Please don't post (un)subscribe requests to the list. Use the
>>>>> www interface or send it to gmx-users-request at gromacs.org.
>>>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>>>
>>>
>>> --
>>> ========================================
>>>
>>> Justin A. Lemkul, Ph.D.
>>> Research Scientist
>>> Department of Biochemistry
>>> Virginia Tech
>>> Blacksburg, VA
>>> jalemkul[at]vt.edu | (540) 231-9080
>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>>
>>> ========================================
>>> --
>>> gmx-users mailing list gmx-users at gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>>> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>>> * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-request at gromacs.org.
>>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>
More information about the gromacs.org_gmx-users
mailing list