Re: inconsistence in NAMD and AMBER energy: AMBER truncted octahedron

From: Vlad Cojocaru (Vlad.Cojocaru_at_eml-r.villa-bosch.de)
Date: Wed Jul 22 2009 - 04:44:45 CDT

Tom,

Yes, this differs but its just a rotation .... These vectors were
written sometime ago on the AMBER mailing list and it took me a while to
find them ...I guess it would be good to include them in the NAMD manual
(also in the paragraph describing simulations with the AMBER force field).

You'll not see it as a nice truncated octahedron in VMD .. you'd rather
see it as a general triclinic box .. probably similar to what ptraj
would give you if you don't specify "familiar" on the image command line
(this, if you are familiar with AMBER). This is because NAMD is
converting all non-cubic boxes to the more general triclinic
representation ...

But this is just a visualization issue ... The simulations are fine ...
And you can get back the truncated octahedron if you image the
trajectory in ptraj using the "familiar" keyword in the image command
line ..

Cheers
Vlad

Thomas C. Bishop wrote:
> Thanks Vlad,
> This differs a bit from info in Chap 3.2 Gromacs manual but appears to be
> consistent (up to a rotation or something) and Gromacs uses only 3 sig digits
> IN THE MANUAL s.t. 4/9 sqrt(3 ) = 0.770
>
> It's been a long time since I tried this.. does the trajectory
> appear "normal" in in VMD, i.e. truncated oct or w/ cells angles of
> 72/109/72?
>
> TOm
>
> On Tuesday 21 July 2009, Vlad Cojocaru wrote:
>
>> Tom, Yi,
>>
>> The correct vectors to run periodic simulations with amber-generated
>> truncated octahedrons in NAMD are:
>> d, 0, 0
>> (-1/3)*d, (2/3)sqrt2*d, 0
>> (-1/3)*d, (-1/3)sqrt2*d, (-1/3)sqrt6*d
>>
>> where d is the cell dimension. The volume of the amber truncated
>> octahedron is 0.7698*d^3 (from amber code)
>> If you use these vectors it should work.
>>
>> Cheers
>> Vlad
>>
>> Yi Shang wrote:
>>
>>> Hi Thomas,
>>> I thought about this before.
>>> It's true you could get very precise cell dimensions to fit your
>>> truncated octahedron. But there is no discriminance in NAMD that could
>>> tell if that's a octahedron or a box (no angle parameters, for
>>> example). It only takes three dimension parameters and assumes it's a
>>> box.
>>>
>>> So, personally I don't think any imaging methods other than box would
>>> work with NAMD. Basically we could not restart amber simulation in
>>> truncated octahedron, rerun it with periodic box instead...
>>>
>>> On Tue, Jul 21, 2009 at 11:14 AM, Thomas C. Bishop <bishop_at_tulane.edu
>>> <mailto:bishop_at_tulane.edu>> wrote:
>>>
>>> The important take home message is as you say.
>>>
>>> "And I suggest all people that are or want to do NAMD using AMBER
>>> parameters
>>> need to keep in mind that VDW energy is NOT going to be the same
>>> if you are
>>> using periodic boundary and NAMD 2.7b1, or older version."
>>>
>>> w/r/t periodic BC's ( i'm doing this from memory) but I thought
>>> that even w/
>>> PBC on I could get the same answer up to tyhe VDW issue. NOTE I
>>> never figured out how to get truncated octahedron to work ( the gromacs
>>> manual has
>>> a nice list of how to convert cell dimensions for various
>>> geometries) but I
>>> thought a cubic cell worked.
>>>
>>>
>>> Tom
>>>
>>> On Tuesday 21 July 2009, Yi Shang wrote:
>>> > Hi all,
>>> > Thank you for replying.
>>> > Thomas, are you using periodic boundary condition or not? Because
>>> > if you look at my email. The energy gap is there only if I use
>>> > periodic boundary condition. I was not using any SHAKE, temperature
>>> > or
>>>
>>> pressure
>>>
>>> > control.
>>> > Goutham, I am using same amber parameter file for both
>>>
>>> simulations when I
>>>
>>> > did comparison, so I believe force parameters are the same.
>>> >
>>> > Floris, you are right, by turning of vdw correction in amber
>>>
>>> simulation
>>>
>>> > (vdwmeth = 0), I get following energies, notice that vdw energy
>>>
>>> is exactly
>>>
>>> > the same as NAMD!! now all energies are consistent...
>>> > NSTEP = 0 TIME(PS) = 0.000 TEMP(K) = 0.00
>>>
>>> PRESS =
>>>
>>> > 0.0
>>> > Etot = -57081.5138 EKtot = 0.0000 EPtot =
>>> > -57081.5138
>>> > BOND = 341.5077 ANGLE = 716.0017 DIHED =
>>> > 2372.7669
>>> > 1-4 NB = 1076.3715 1-4 EEL = 9100.8801 VDWAALS =
>>> > 6132.6047
>>> > EELEC = -76821.6464 EHBOND = 0.0000 RESTRAINT =
>>> > 0.0000
>>> > Ewald error estimate: 0.2976E-03
>>> > Now I need to think about how to get away with it. And I suggest
>>>
>>> all people
>>>
>>> > that are or want to do NAMD using AMBER parameters need to keep
>>>
>>> in mind
>>>
>>> > that VDW energy is NOT going to be the same if you are using
>>>
>>> periodic
>>>
>>> > boundary and NAMD 2.7b1, or older version.
>>> >
>>> > On Tue, Jul 21, 2009 at 10:33 AM, Thomas C. Bishop
>>>
>>> <bishop_at_tulane.edu <mailto:bishop_at_tulane.edu>>wrote:
>>> > > I have been able to obtain the exact same energies from NAMD
>>>
>>> and amber
>>>
>>> > > for several time steps using suggested parameters in the NAMD
>>>
>>> manual.
>>>
>>> > > For this it is easiest to start with a temperature of 0 but
>>>
>>> you can also
>>>
>>> > > read
>>> > > in a velocity file s.t. they have same coords AND momentum
>>> > >
>>> > > Obviously you have to turn off anything that relies on random
>>>
>>> number
>>>
>>> > > generator
>>> > > so NVE it primary target here.
>>> > >
>>> > > bonds, angle, dihedrals, improper agree w/out any problem
>>> > >
>>> > > vdw and and elec are tricky.
>>> > > first start with a simple cut-off to make sure your parameters
>>> > > are consistent.
>>> > > then turn on long range interactions. PAY CAREFUL ATTENTION TO
>>>
>>> NAMD
>>>
>>> > > SWITCHING
>>> > > OPTIONS! see manual on swithcing functions too. This
>>>
>>> introduces a shift
>>>
>>> > > in the elec and vdw terms that will make it look like NAMD and
>>>
>>> AMBER are
>>>
>>> > > radically different but in terms of forces they are nearly
>>>
>>> same. make
>>>
>>> > > switching agree and energies agree also.
>>> > >
>>> > > finally amber includes a long range vdw correction that will
>>>
>>> affect the
>>>
>>> > > trajectory over time but in first few steps ~10 this is
>>>
>>> minimal. W/
>>>
>>> > > cutoff they should give almost same result for much longer
>>> > >
>>> > > Please post what you finally come up w/
>>> > >
>>> > > Tom
>>> > >
>>> > > On Tuesday 21 July 2009, Floris Buelens wrote:
>>> > > > Hi,
>>> > > >
>>> > > > I'm not sure about exactly how AMBER handles this issue, but
>>>
>>> I suspect
>>>
>>> > > the
>>> > >
>>> > > > difference might arise from the neglect of van der Waals
>>>
>>> interactions
>>>
>>> > > > beyond the cutoff. I think AMBER applies an analytical
>>>
>>> correction by
>>>
>>> > > > default for these interactions (search the manual for
>>>
>>> 'vdwmeth'), while
>>>
>>> > > > NAMD as of 2.7b1 doesn't yet. I would suggest re-running
>>>
>>> your analysis
>>>
>>> > > > using an increasing series of cutoff values, I think you
>>>
>>> should see the
>>>
>>> > > > NAMD value converging towards that of AMBER. Alternatively,
>>>
>>> the option
>>>
>>> > > for
>>> > >
>>> > > > analytical corrections for long-range van der Waals has now
>>>
>>> been added
>>>
>>> > > > to the latest development code, which you can access through
>>>
>>> CVS as
>>>
>>> > > described
>>> > >
>>> > > > on the NAMD website, if you use that code with the option
>>>
>>> "LJCorrection
>>>
>>> > > on"
>>> > >
>>> > > > I hope you shouldn't see such a big discrepancy. Best wishes,
>>> > > >
>>> > > > Floris
>>> > > >
>>> > > >
>>> > > > ________________________________
>>> > > > From: Yi Shang <mirandaisbest_at_gmail.com
>>>
>>> <mailto:mirandaisbest_at_gmail.com>>
>>>
>>> > > > To: namd-l_at_ks.uiuc.edu <mailto:namd-l_at_ks.uiuc.edu>
>>> > > > Sent: Tuesday, 21 July, 2009 1:31:47
>>> > > > Subject: namd-l: inconsistence in NAMD and AMBER energy
>>> > > >
>>> > > >
>>> > > > Hi all,
>>> > > > I did some comparison of energy calculated by NAMD 2.7b1 and
>>>
>>> AMBER 10.
>>>
>>> > > > I found something intriguing. During simulation, VDW energy
>>>
>>> from AMBER
>>>
>>> > > > simulation is constantly much lower from NAMD. So I tried to
>>>
>>> compare
>>>
>>> > > point
>>> > >
>>> > > > energy first. I turned off SHAKE, temperature, and pressure
>>>
>>> control and
>>>
>>> > > > only looked at basic energy calculation of one single
>>>
>>> structure:
>>> > > > bonds, angles, dihedrals, electrostatic, and vdw energies. For
>>> > > > non-periodic boundary, I had no problem, energies are
>>>
>>> exactly the same.
>>>
>>> > > > For periodic boundary condition here comes the problem:
>>> > > > -electro energy off by ~14kcal/mol, which is ok, parameters
>>>
>>> might not
>>>
>>> > > > be identical in two methods. -vdw energy off by > 600
>>>
>>> kcal/mol. That's
>>>
>>> > > > terrible.
>>> > > >
>>> > > > I tested these inputs on two different periodic systems, and
>>>
>>> I get the
>>>
>>> > > same
>>> > >
>>> > > > result (same energy gap). So I don't think it's system specific
>>> > > > problem. Here is my inputs and outputs from one system's
>>>
>>> comparison (I
>>>
>>> > > > combined AMBER VDWAALS and 1-4 NB energies to compare to
>>>
>>> NAMD vdw
>>>
>>> > > > energy, likewise for electrostatic interaction):
>>> > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>>> > > > NAMD config.in <http://config.in/>:
>>> > > > amber on
>>> > > > parmfile box.top
>>> > > > ambercoor box.crd
>>> > > > binaryoutput no
>>> > > > outputname md6
>>> > > > outputpressure 1
>>> > > > outputenergies 1
>>> > > > exclude scaled1-4
>>> > > > 1-4scaling 0.83333
>>> > > > timestep 2
>>> > > > temperature 0
>>> > > > cutoff 8
>>> > > > switching off
>>> > > > PME yes
>>> > > > PMEGridSpacing 1.0
>>> > > > cellbasisvector1 79.6614 0 0
>>> > > > cellbasisvector2 0 60.2065 0
>>> > > > cellbasisvector3 0 0 63.7403
>>> > > > cellorigin 0 0 0
>>> > > > run 0
>>> > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>>> > > > amber md.in <http://md.in/>:
>>> > > > &cntrl
>>> > > > imin = 0, nstlim = 1, dt = 0.002,
>>> > > > irest = 0, ntx = 1,
>>> > > > tempi = 0.0, temp0 = 0.0,
>>> > > > ntwx = 0, ntwe = 0, ntwr = 0, ntpr = 1,
>>> > > > cut = 8.0,
>>> > > > /
>>> > > >
>>> > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>>> > > > namd energies:
>>> > > > ETITLE: TS BOND ANGLE DIHED
>>> > >
>>> > > IMPRP
>>> > >
>>> > > > ELECT VDW BOUNDARY
>>> > > > MISC KINETIC TOTAL TEMP POTENTIAL TOTAL3
>>> > > > TEMPAVG PRESSURE GPRESSURE VOL UME
>>> > > > PRESSAVG GPRESSAVG
>>> > > > ENERGY: 0 341.5078 716.0018 2372.7673
>>> > >
>>> > > 0.0000
>>> > >
>>> > > > -67739.0356 7208.9795 0.0000
>>>
>>> 0.0000 0.0000
>>>
>>> > > > -57099.7791 0.0000 -57099.7791 -56529.2524
>>> > >
>>> > > 0.0000
>>> > >
>>> > > > 11740.3005 -276.1608 305707.0 250
>>>
>>> 11740.3005
>>>
>>> > > > -276.1608
>>> > > >
>>> > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>>> > > > amber energies:
>>> > > > NSTEP = 0 TIME(PS) = 0.000 TEMP(K) =
>>>
>>> 0.00 PRESS =
>>>
>>> > > > 0.0 Etot = -57686.8091 EKtot = 0.0000 EPtot
>>>
>>> =
>>>
>>> > > > -57686.8091 BOND = 341.5077 ANGLE =
>>>
>>> 716.0017 DIHED
>>>
>>> > > =
>>> > >
>>> > > > 2372.7669 1-4 NB = 1076.3715 1-4 EEL =
>>> > > > 9100.8801 VDWAALS = 5527.3094 EELEC = -76821.6464
>>> > > > EHBOND = 0.0000 RESTRAINT = 0.0000 Ewald error
>>> > > > estimate:
>>>
>>> 0.2976E-03
>>>
>>> > > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>>> > > >
>>> > > > It's strange because in NAMD seems most important parameter
>>>
>>> I would
>>>
>>> > > > play around with VDW energy is switching, and it's turned
>>>
>>> off. Same
>>>
>>> > > > cutoff is set in two inputs, but still, I get very different
>>>
>>> results.
>>>
>>> > > > any
>>> > >
>>> > > suggestions
>>> > >
>>> > > > would be appreciated! Thanks!
>>> > >
>>> > > --
>>> > > **********************
>>> > > Thomas C. Bishop *
>>> > > Office: 504-862-3370 *
>>> > > Fax: 504-862-8392 *
>>> > > **********************
>>>
>>> --
>>> **********************
>>> Thomas C. Bishop *
>>> Office: 504-862-3370 *
>>> Fax: 504-862-8392 *
>>> **********************
>>>
>>>
>>>
>>>
>>> --
>>> Miranda
>>>
>
>
>
>

-- 
----------------------------------------------------------------------------
Dr. Vlad Cojocaru
EML Research gGmbH
Schloss-Wolfsbrunnenweg 33
69118 Heidelberg
Tel: ++49-6221-533202
Fax: ++49-6221-533298
e-mail:Vlad.Cojocaru[at]eml-r.villa-bosch.de
http://projects.villa-bosch.de/mcm/people/cojocaru/
----------------------------------------------------------------------------
EML Research gGmbH
Amtgericht Mannheim / HRB 337446
Managing Partner: Dr. h.c. Klaus Tschira
Scientific and Managing Director: Prof. Dr.-Ing. Andreas Reuter
http://www.eml-r.org
----------------------------------------------------------------------------

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:53:04 CST