AMBER Archive (2002)

Subject: Perturbation with sander : trouble in keeping perturbed part bonded to unperturbed (fwd)

From: amber_at_heimdal.compchem.ucsf.edu
Date: Mon Aug 26 2002 - 16:50:02 CDT

---------- Forwarded message ----------
Date: Thu, 22 Aug 2002 14:22:28 +0200
From: hagop demirdjian <hagop.demirdjian_at_ens.fr>
To: amber list <amber_at_heimdal.compchem.ucsf.edu>
Subject: Perturbation with sander : trouble in keeping perturbed part
    bonded to unperturbed

Hi AMBER community,

I am a phd student currently learning thermodynamic integration
techniques to compute free energies. The system I chose is acetic acid
perturbed to formic acid in a 511 water box.

I built acetic acid using bcc charges from antechamber, and equilibrated
the system for 30 ps. It went fine (energy decreased to a plateau,
density reached 1). Bond lenght remained unchanged but angles varied a
little : HC-CT-HC went from 109° to 106-108 ; OH-C -O went from 129.5°
to 125 and dihedral HC-CT-C -O from 0° to 11°.

I then used the rstrt file to run equilibration for a perturbation run
at lambda=0.
Large jumps in energy (100 kCal rise in a few steps when fluctuations
are only about 5 kCal) occur after a few ps. The energy then relaxes for
a few more ps before the next jump.
I visualized the acid with VMD and gOpenmol at several points of the
simulation and both softwares describe it with a distorted
geometry. Just after an energy jump, it is very distorted (methyl C at
0.3 A from carbonyl C instead of 1.5 ; linked to both oxygens ; methyl
hydrogens not linked to C and between 3 to 6 A). When energy relaxes,
distortion is less important : one H may be far from the methyl...

I checked the .top file created with saveamberparmpert : the perturbed
part seems fine (normal number of perturbed atoms, bonds....) but the
unperturbed part seems a bit odd to me : the part corresponding to the
methyl fragment does not appear anymore.
e.g

%FLAG BONDS_WITHOUT_HYDROGEN

%FORMAT(10I8)

        3 6 2 3 9 3 0 3 4

from unperturbed.top

becomes

%FLAG BONDS_WITHOUT_HYDROGEN

%FORMAT(10I8)

        3 6 2 3 9 3

in perturbed.top

(here, that means that CT-C bond is not taken in account anymore)

Maybe this is the normal way the perturbed.top files are build, but it
may be the sign that something is wrong with my frcmod file.

I used the following one :

       acetic to formic acid
MASS
Z2 1.008 Z2 IS DUMMY PROTON (HYDROGEN)
Z1 12.01 Z1 IS aliphatic carbon linked hydrogen

NONB
Z2 0.0 0.0
Z1 0.0 0.0

BOND
Z1-Z2 331. 1.080
Z1-C 331. 1.080

ANGLE
Z2-Z1-Z2 35. 109.5
C-Z1-Z2 35. 109.5
Z1-C-O 50. 120.
Z1-C-OH 54. 113.7

DIHE
OH-C-Z1-Z2 1 0.0 0. 2.
OH-C-Z1-Z2 1 0.0 0. 1.
O-C-Z1-Z2 1 0.0 0. 3.
O-C-Z1-Z2 1 0.0 0. 1.
OH-C-CT-HC 1 0.0 0. 1.
OH-C-CT-HC 1 0.0 0. 2.

I took the parameters from PARM99.dat except for DIHE which were
suggested by xleap when I tried to saveamberparmpert.

I used a formatted file too (using a tutorial template) because I was
wondering if something was not read in. There were no differences in
output files.

Do you have any idea of how to keep the methyl fragment bonded to the
unperturbed one or do you see any obvious errors in the way I build the
system ?

I also ran T_only equilibration and tried several tolerances for SHAKE
(0.0005, 0.00001 and 0.000001). Did runs with lamda =0.5 too. Same problems.

Here is my input file in case that might give a clue :

  acetic acid in water
  &cntrl
        ntx = 5,
        irest = 1,
        ntwe = 1,

        ntf = 3,
        ntb =2,

        nstlim = 20000,

        tempi = 300.0,
        ntt = 1,

        ntp=1,

        ntc =3,
        jfastw = 0,

        icfe = 1,
        clambda = 0.00000,
        klambda = 4,
  &end

Any help would be greatly appreciated,
thank you,
Hagop.