AMBER Archive (2003)

Subject: AMBER: PMEMD NMR Restraints fix

From: Robert Duke (rduke_at_email.unc.edu)
Date: Fri Aug 08 2003 - 17:13:08 CDT

Dear Amber-ites,
I have fixed the PMEMD bug reported by Kristina Furst, and have appended the
mail I sent to her below. I originally cc'd the fix to the amber reflector,
but it probably got bounced due to the presence of the attached fixed source
file, pmemd.src/pmemd/force_cit.f90. There is a 1 line diff indicated in
the mail below which may be used to do the fix, or if you prefer, you can
send me mail requesting a copy of the whole file.
Regards - Bob Duke

Original mail thread with Kristina, sans attachment:

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
>
>

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