Dear Prof. Case and AMBER list,
My question is regarding the Thermodynamic Integration implementation
in AMBER 9. In the amber mailing list, Petr Kulháne kwrote that in a TI
calculation, mass changes are not reflected in perturbation calculations.
http://amber.ch.ic.ac.uk/archive/200504/0432.html
I am doing perturbation calculations on some systems where O -> N or N ->
O perturbations are done. I was testing one of my system, to see if there
is any anomalities in it. What I have done is the following: I have two
prmtop files, prmtop1 and prmtop2, which represent the initial and final
states of the system. I have also an inpcrd file which is the coordinate
file for the simulation.
I have done two 1-step simulation;
i) from prmtop1 to prmtop2 (lambda=0, klambda=6, icfe=2),
ii) from prmtop2 to prmtop1 (lambda=1, klambda=6, icfe=2).
The energy results should be the same for these two test simulations, but
I am getting different kinetic energy results. Here is the 1-step energy
results:
--- prmtop1 to prmtop2 with lambda=0 ------------------
NSTEP = 0 TIME(PS) = 120.000 TEMP(K) = 36.05 PRESS = 0.0
Etot = -27314.0287 EKtot = 9137.4836 EPtot = -36451.5123
BOND = 143.3686 ANGLE = 252.1235 DIHED = 413.6544
1-4 NB = 148.3946 1-4 EEL = -2256.9910 VDWAALS = 4269.4955
EELEC = -39421.5580 EHBOND = 0.0000 RESTRAINT = 0.0000
Ewald error estimate: 0.4419E-04
--------------------------------------------------------------------
NSTEP = 1 TIME(PS) = 120.001 TEMP(K) = 330.42 PRESS = 0.0
Etot = -29527.6465 EKtot = 6923.8658 EPtot = -36451.5123
BOND = 143.3686 ANGLE = 252.1235 DIHED = 413.6544
1-4 NB = 148.3946 1-4 EEL = -2256.9910 VDWAALS = 4269.4955
EELEC = -39421.5580 EHBOND = 0.0000 RESTRAINT = 0.0000
Ewald error estimate: 0.4419E-04
--------------------------------------------------------------------
--- prmtop2 to prmtop1 with lambda=1 ------------------
NSTEP = 0 TIME(PS) = 120.000 TEMP(K) = 436.05 PRESS = 0.0
Etot = -27314.0287 EKtot = 9137.4836 EPtot = -36451.5123
BOND = 143.3686 ANGLE = 252.1235 DIHED = 413.6544
1-4 NB = 148.3946 1-4 EEL = -2256.9910 VDWAALS = 4269.4955
EELEC = -39421.5580 EHBOND = 0.0000 RESTRAINT = 0.0000
Ewald error estimate: 0.4419E-04
--------------------------------------------------------------------
NSTEP = 1 TIME(PS) = 120.001 TEMP(K) = 330.35 PRESS = 0.0
Etot = -29529.0292 EKtot = 6922.4831 EPtot = -36451.5123
BOND = 143.3686 ANGLE = 252.1235 DIHED = 413.6544
1-4 NB = 148.3946 1-4 EEL = -2256.9910 VDWAALS = 4269.4955
EELEC = -39421.5580 EHBOND = 0.0000 RESTRAINT = 0.0000
Ewald error estimate: 0.4419E-04
--------------------------------------------------------------------
The 0th step has the same potential/kinetic energy but at step=1, the
kinetic energy changes.
When the perturbation is from H->DH (where DH is a dummy atom with an
H-mass), everything is fine; energies are the same for both systems.
Assuming Hamiltonian mechanics is used in calculations of the velocities
and coordinates in AMBEM 9, I can see what might be causing this. As Petr
wrote in one of his emails (above link), this might be due to the mass
difference of the atoms that are perturbed. Does that mean that I cannot
do any mass perturbation (like a perturbation from oxygen to nitrogen) in
TI Approach?
If both states (initial and final states) have the same conformations, I
would assume that this should not create any problem, but I would like to
know your ideas about this. I will greatly appreciate your comments on
this.
Best regards,
--
Ilyas Yildirim
---------------------------------------------------------------
= Department of Chemistry - =
= University of Rochester - =
= Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) =
= http://www.pas.rochester.edu/~yildirim/ =
---------------------------------------------------------------
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Dec 12 2007 - 06:07:20 PST
