AMBER Archive (2003)

Subject: RE: AMBER: PMEMD and EAMBER

From: Kristina Furse (kristina.e.furse_at_vanderbilt.edu)
Date: Fri Aug 08 2003 - 13:33:28 CDT

Great--thanks! I'll check it out and let you know.

The constraint is definitely holding like it should with pmemd (pattern of
torsion vs. time very similar to sander7 test), and the energy profile is
similar as well, though the average potential energy over 10ps is slightly
lower with pmemd (-351752.7754 pmemd vs. -351716.9462 sander7). And I am very
very glad you re-included this option--even though I'm on IRIX not LINUX,
you're still saving me a week per month for running this monster!

Kristina

>===== Original Message From "Robert Duke" <rduke_at_email.unc.edu> =====
>Kristina -
>
>I got a chance to look at your example runs and the pmemd code this morning,
>and it is actually really a simple problem. I reordered some code in
>force_cit.f90 in an effort to make things a little more clear, and screwed
>up in regard to when ene(1) gets updated if nmropt is nonzero. The simple
>one line fix to src.pmemd/pmemd/force_cit.f90 is:
>664a665
>> ene(1) = ene(1) + enmr(1) + enmr(2) + enmr(3)
>
>I have also attached the modified version of force_cit.f90 for your
>convenience, though that may give folks on the amber reflector grief.
>
>I am virtually certain this fixes the problem, but if you would rebuild,
>test it, and let me know, I would appreciate it. Normally, I like to test
>stuff myself, but this one was dead easy, if you don't mind rebuilding and
>testing yourself, I think that would be best. Sorry for messing it up. At
>one point in time, I had actually removed support for nmropt == 1 from
>pmemd, since I was not sure how the coordinate and forces distribution
>issues were going to play out there. I ended up re-including it when I
>determined that it could be done relatively easily, and our temperature
>scaling tests worked. If everyone would check nmropt== 1 results carefully,
>though, I would appreciate it. I simply don't have test coverage of the
>gazillion different possible options. I would note that in my perusal of
>the original amber 6 code, I noticed that pme breaks nmropt==1- based
>modification of cut, so I don't even trust that all of this stuff still
>works in amber 6 (could be thats the only problem, but if you don't have a
>complete test suite, you DON'T KNOW).
>
>I expect that in the next week or so I will include this fix in a new drop
>of pmemd that also includes support for the various alphaservers and the
>cray t3e at PSC.
>
>Regards - Bob
>
>----- Original Message -----
>From: "Kristina Furse" <kristina.e.furse_at_vanderbilt.edu>
>To: "Robert Duke" <rduke_at_email.unc.edu>
>Sent: Friday, August 08, 2003 1:21 PM
>Subject: RE: AMBER: PMEMD and EAMBER
>
>
>> Thank you so much! My system is pretty big--104,234 atoms. If that sounds
>like
>> too much trouble to work with, I could do a quick test with a smaller
>system
>> and send that. If not, I'll ship the big guy.
>>
>> Kristina
>>
>> >===== Original Message From "Robert Duke" <rduke_at_email.unc.edu> =====
>> >Kristina -
>> >It is very possible that something (like an enmr) got whacked in the
>pmemd
>> >code work without getting noticed. The nmropt == 1 stuff is not
>extensively
>> >tested; we mostly use the temperature stepping capabilities in our group.
>> >Non-nmr constraints were tested and do work (hb benchmark), producing the
>> >same results for sander 6, 7, pmemd, and pmemd 7-compatibility mode. If
>you
>> >will send me your input (gzipped, please) I will be glad to track down
>and
>> >fix the bug. I'll do this for any nmropt == 1 stuff folks turn up; I have
>> >tried to not touch the nmropt == 1 code because of the testing impact,
>but
>> >it is always possible that something obscure got missed.
>> >Regards - Bob
>> >----- Original Message -----
>> >From: "Kristina Furse" <kristina.e.furse_at_vanderbilt.edu>
>> >To: <amber_at_scripps.edu>
>> >Sent: Thursday, August 07, 2003 4:32 PM
>> >Subject: AMBER: PMEMD and EAMBER
>> >
>> >
>> >> Hello-
>> >>
>> >> We've been trying out PMEMD, using it for some potential of mean force
>> >> calculations (using NMR-style torsional constraints to bias the
>sampling,
>> >then
>> >> a WHAM method to get the PMF curve).
>> >>
>> >> I ran some initial test runs to compare pmemd with sander in AMBER7
>> >(running
>> >> on an SGI origin). While comparing the test runs, I noticed something
>> >about
>> >> the computation of EAMBER in minimization and dynamics output (below).
>> >Looking
>> >> at the first step of minimization output for the very same input
>> >coordinates
>> >> run in sander (amber7) and pmemd (amber7_compat=1), the energies are
>> >identical
>> >> (well, to be exact, EEL is different by 0.0001) except for EAMBER,
>which
>> >is
>> >> off by 2.5778kcal, which is exactly the restraint energy minus the
>0.0001.
>> >I
>> >> summed the energies by hand and it looks like PMEMD may be subtracting
>the
>> >> restraint energy twice, or else not including it in the total potential
>> >energy
>> >> calculation then still subtracting it to get EAMBER?
>> >>
>> >> I took a look at some dynamics output (same input coordinates as
>before)
>> >from
>> >> sander7 and pmemd, and I think the same thing may be going on there.
>For
>> >> pmemd, it seems like the number for EPtot already does not include
>> >restraint
>> >> energy, then EAMBER has restraint energy subtracted again.
>> >>
>> >> Has anyone else noticed this? Did I miss something somewhere, or am I
>just
>> >> adding wrong...
>> >>
>> >> Thanks!
>> >>
>> >> Kristina
>> >>
>> >>
>> >>
>> >> Minimization output:
>> >>
>> >> Sander (amber 7):
>> >>
>> >> NSTEP ENERGY RMS GMAX NAME
>NUMBER
>> >> 1 -3.5193E+05 1.6203E+01 1.3187E+02 CD
>17121
>> >>
>> >> BOND = 3402.1061 ANGLE = 9311.0284 DIHED =
>> >11245.0567
>> >> VDWAALS = 48314.8917 EEL = -472531.9934 HBOND =
>> >0.0000
>> >> 1-4 VDW = 4345.3096 1-4 EEL = 43982.3959 RESTRAINT =
>> >2.5779
>> >> EAMBER = -351931.2050
>> >>
>>
>>===========================================================================
>=
>> >==
>> >> =
>> >> NMR restraints for step 1
>> >> Energy (this step): Bond = 2.578 Angle = 0.000 Torsion =
>> >0.000
>> >> Energy (tot. run) : Bond = 2.578 Angle = 0.000 Torsion =
>> >0.000
>> >>
>> >> DEVIATIONS: Target=(r2+r3)/2 Target = closer of
>r2/r3
>> >> This step Entire run This step Entire
>> >run
>> >> ave. rms ave. rms ave. rms ave.
>> >rms
>> >> Bond 0.253 0.254 0.253 0.254 0.253 0.254 0.253
>> >0.254
>> >> Angle 0.000 0.000 0.000 0.000 0.000 0.000 0.000
>> >0.000
>> >> Torsion 0.000 0.000 0.000 0.000 0.000 0.000 0.000
>> >0.000
>> >>
>>
>>===========================================================================
>=
>> >==
>> >> =
>> >>
>> >>
>> >>
>> >>
>> >> PMEMD (amber7_compat=1):
>> >>
>> >> NSTEP ENERGY RMS GMAX NAME
>NUMBER
>> >> 1 -3.5193E+05 1.6203E+01 1.3187E+02 CD
>17121
>> >>
>> >> BOND = 3402.1061 ANGLE = 9311.0284 DIHED =
>> >11245.0567
>> >> VDWAALS = 48314.8917 EEL = -472531.9933 HBOND =
>> >0.0000
>> >> 1-4 VDW = 4345.3096 1-4 EEL = 43982.3959 RESTRAINT =
>> >2.5779
>> >> EAMBER = -351933.7828
>> >>
>>
>>===========================================================================
>=
>> >==
>> >> =
>> >> NMR restraints for step 1
>> >> Energy (this step): Bond = 2.578 Angle = 0.000 Torsion =
>> >0.000
>> >> Energy (tot. run) : Bond = 2.578 Angle = 0.000 Torsion =
>> >0.000
>> >>
>> >> DEVIATIONS: Target=(r2+r3)/2 Target = closer of
>r2/r3
>> >> This step Entire run This step Entire
>> >run
>> >> ave. rms ave. rms ave. rms ave.
>> >rms
>> >> Bond 0.253 0.254 0.253 0.254 0.253 0.254 0.253
>> >0.254
>> >> Angle 0.000 0.000 0.000 0.000 0.000 0.000 0.000
>> >0.000
>> >> Torsion 0.000 0.000 0.000 0.000 0.000 0.000 0.000
>> >0.000
>> >>
>>
>>===========================================================================
>=
>> >==
>> >> =
>> >>
>> >>
>> >>
>> >> Molecular Dynamics Output:
>> >>
>> >> Sander (amber7)
>> >>
>> >> NSTEP = 100 TIME(PS) = 3150.100 TEMP(K) = 296.72 PRESS =
>32.9
>> >> Etot = -287279.7774 EKtot = 64194.7299 EPtot
> = -351474.5072
>> >> BOND = 3498.2008 ANGLE = 9330.0143 DIHED =
>11174.5893
>> >> 1-4 NB = 4342.5808 1-4 EEL = 44008.7977 VDWAALS =
>48274.1435
>> >> EELEC = -472109.4246 EHBOND = 0.0000 RESTRAINT =
>6.5909
>> >> EAMBER (non-restraint) = -351481.0982
>> >> EKCMT = 25330.3166 VIRIAL = 24605.6087 VOLUME =
>1019086.6920
>> >> Density =
>1.0541
>> >> Ewald error estimate: 0.4495E-04
>> >>
>> >>
>> >>
>> >> PMEMD (amber7_compat=1)
>> >>
>> >> NSTEP = 100 TIME(PS) = 3150.100 TEMP(K) = 296.72 PRESS =
>32.8
>> >> Etot = -287292.8061 EKtot = 64195.1417 EPtot
> = -351487.9478
>> >> BOND = 3498.2009 ANGLE = 9330.0138 DIHED =
>11174.5889
>> >> 1-4 NB = 4342.5814 1-4 EEL = 44008.7981 VDWAALS =
>48276.0446
>> >> EELEC = -472118.1755 EHBOND = 0.0000 RESTRAINT =
>6.5910
>> >> EAMBER (non-restraint) = -351494.5388
>> >> EKCMT = 25330.5359 VIRIAL = 24608.2331 VOLUME =
>1019085.5924
>> >> Density =
>1.0541
>> >> Ewald error estimate: 0.4460E-04
>> >>
>> >>
>> >>
>> >> **In case you wonder about the difference in restraint energies, the
>> >dynamics
>> >> runs have NMR-style constraints on the torsion we're driving, as well
>as a
>> >> "virtual bond" to a heme group accomplished using an NMR-style distance
>> >> constraint. The minimizations just have the virtual bond.
>> >>
>> >>
>> >>
>> >>
>> >> Finally, for good measure, dynamics input files:
>> >>
>> >>
>> >> sander (amber7)
>> >>
>> >> Biased PMF: Cox dimer w/radical, explicit water, 9.0 cut
>> >> &cntrl
>> >>
>> >> nmropt = 1,
>> >> ntx = 5, irest = 1, ntrx = 1, ntxo = 1,
>> >> ntpr = 100, ntwx = 100, ntwv = 0, ntwe = 0,
>> >>
>> >> ntf = 2, ntb = 2,
>> >> cut = 9.0,
>> >>
>> >> ibelly = 0, ntr = 0,
>> >>
>> >> imin = 0,
>> >> nstlim = 10000,
>> >> nscm = 100,
>> >> t = 0.0, dt = 0.001,
>> >>
>> >> temp0 = 298.0, tempi = 298.0,
>> >> ig = 71277, heat = 0.0,
>> >> ntt = 1,
>> >> tautp = 1.0,
>> >> vlimit = 20.0,
>> >>
>> >> ntp = 1, pres0 = 1.0, comp = 44.6,
>> >> taup = 1.0,
>> >>
>> >> ntc = 2, tol = 0.00001, watnam = 'SPC ',
>> >>
>> >> &end
>> >> &ewald
>> >>
>> >> nbflag = 1, skinnb = 2.0, use_pme = 1,
>> >>
>> >> &end
>> >> &wt type='END' &end
>> >> DISANG=tor0.rst
>> >>
>> >>
>> >> pmemd input file just has the additional flag, amber7_compat=1,
>> >>
>> >> ****************************************************
>> >> Kristina E. Furse
>> >> Department of Chemistry
>> >> Center for Structural Biology
>> >> Vanderbilt University
>> >> email: kfurse_at_structbio.vanderbilt.edu
>> >>
>> >>
>> >> -----------------------------------------------------------------------
>> >> The AMBER Mail Reflector
>> >> To post, send mail to amber_at_scripps.edu
>> >> To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu
>> >>
>> >>
>> >
>> >
>> >
>> >-----------------------------------------------------------------------
>> >The AMBER Mail Reflector
>> >To post, send mail to amber_at_scripps.edu
>> >To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu
>>
>> ****************************************************
>> Kristina E. Furse
>> Department of Chemistry
>> Center for Structural Biology
>> Vanderbilt University
>> email: kfurse_at_structbio.vanderbilt.edu
>>
>>

****************************************************
Kristina E. Furse
Department of Chemistry
Center for Structural Biology
Vanderbilt University
email: kfurse_at_structbio.vanderbilt.edu

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu