AMBER Archive (2007)

Subject: AMBER: polarization and water density

From: Thérèse Malliavin (terez_at_pasteur.fr)
Date: Fri Apr 13 2007 - 04:36:44 CDT

Dear AMBER users,

I am trying since several days to run an MD simulation at constant
pressure (1atm), with periodic boundary conditions, PME and polarization
(force-field ff02) and a Carr-Parrinello scheme for the integration of the
dipole motion (INDMETH=3). A Langevin dynamics is used for the temperature
scaling (NTT=3 and GAMMA_LN=2). The SHAKE option is on, as required when
using polarization (NTC=2), and the timestep is 1fs and was even put to
0.5 fs (DT=.0005) without removing the problem.

My problem is that the density of water always increased up to around
1.1-1.15 and then the simulation stops, usually because of "vlimit
exceeded for step". I tried to change by hand the size of the box in the
initial restart file. This reduces the density, but after some simulation
timesteps, it starts to increase again and to become too big.
The timestep was reduced to 0.5 fs without improvig the situation.

The general setup of the equilibration is the following.

A) First, tleap to build the system:

solute = LoadPDB ./xaa
x = copy solute
addIons x Cl- 0
y = copy TIP3PBOX
solvatebox x y 10
saveamberparmpol x bmr4395.top bmr4395.crd

B) Then, 6 runs of minimisation where the solute is restrained with
decreasing the restraint energy constant from 100.0 to 2.5. A typical
input file is the following:

equilibrating water with restrained solute
 &cntrl
   imin=1, ntx=1, irest=0, ntrx=1, ntxo=1, ntpr=50, ntf=1, ntb=1,
   cut=10.0, nsnb=10, scee=1.2, ntr=0, maxcyc=500, ntmin=1,
   ipol=1, ntc = 2,
   restraint_wt=2.5, restraintmask=':1-101',
 &end
 &ewald
   indmeth = 3, irstdip = 0,
 &end

C) Thermalisation at constant volume to increase the temperature from
0 to 298 K, while keeping the solute restrained.

Thermalization from 0 to 298 K, at constant volume
 &cntrl
    imin = 0, ntx = 1, irest = 0,
    ntpr = 100, ntwx = 100,
    nscm = 100, ntf = 2, ntb = 1,
    nsnb = 10,
    scee = 1.2, ntr = 1, nstlim = 5000,
    cut = 10.0, t = 0.0, dt =.001,
    tempi = 0.0, temp0 = 298.0, ntt = 1,
    ig = 71277, restraint_wt=25.0, restraintmask=':1-101',
    ipol=1, ntc = 2,
 &end
 &ewald
   indmeth = 3, irstdip = 0,
 &end

D) One run of equilibration at constant volume:

Equilibration at 298 K at constant volume
 &cntrl
    imin = 0, ntx = 5, irest = 1,
    ntpr = 100, ntwx = 100,
    nscm = 100, ntf = 2, ntb = 1,
    nsnb = 10, scee = 1.2,
    ntr = 1, nstlim = 2500, ntp = 0,
    cut = 10.0, dt =.001,
    temp0 = 298.0, ntt = 3, gamma_ln = 2
    ipol = 1, ntc = 2,
    ig = 71277, restraint_wt=25.0, restraintmask=':1-101',
 &end
 &ewald
   indmeth = 3, irstdip = 0,
 &end

E) Five runs of equilibration at constant pressure, the restraint
on the solute decreasing from 25 to 2.5:

Equilibration at 298 K at constant pressure
 &cntrl
    imin = 0, ntx = 5, irest = 1,
    ntpr = 100, ntwx = 100,
    nscm = 100, ntf = 2, ntb = 2,
    nsnb = 10, scee = 1.2,
    ntr = 1, nstlim = 5000,
    ntp = 1, pres0 = 1.0, taup = 0.2,
    cut = 10.0, dt =.001,
    temp0 = 298.0, ntt = 3, gamma_ln = 2
    ipol = 1, ntc = 2,
    ig = 71277, restraint_wt=2.5, restraintmask=':1-101',
 &end
 &ewald
   indmeth = 3, irstdip = 0,
 &end

At the end of these runs, the water density is 0.9158. So
I run an additionnal equilibration without restraints on the
solute:

Equilibration at 298 K at constant pressure
 &cntrl
    imin = 0, ntx = 5, irest = 1,
    ntpr = 100, ntwx = 100,
    nscm = 100, ntf = 2, ntb = 2,
    nsnb = 10, scee = 1.2,
    ntr = 0, nstlim = 30000,
    ntp = 1, pres0 = 1.0, taup = 0.2,
    cut = 10.0, dt =.001,
    temp0 = 298.0, ntt = 3, gamma_ln = 2
    ipol = 1, ntc = 2,
    ig = 71277,
 &end
 &ewald
   indmeth = 3, irstdip = 0,
 &end

The density of water starts to increase as soon as the restraints on the
solute are removed: mean density: 1.0440 +/- 0.0431, and final
density: 1.0907. And, when the production run starts, after sometime the
density is too high and it stops.

The system I study is a protein of 94 amino-acids. It is one conformer of
an NMR structure, of medium quality (the PROCHECK report is below), but
200 ps of MD simulation without polarization was run without problem on
it. This protein has a beta-barrel structure, and the polarizable
forcefield may induce a repulsive unstability between the beta sheet
dipoles?

To conclude this (long) mail, I would like to know whether simulations
were run with AMBER using the polarization on proteins, and what are the
recommanded conditions.

 +----------<<< P R O C H E C K S U M M A R Y
 | xaa 1.5 94 residues |
|
*| Ramachandran plot: 67.5% core 20.5% allow 4.8% gener 7.2% disall |
*| All Ramachandrans: 15 labelled residues (out of 89) |
 | Chi1-chi2 plots: 0 labelled residues (out of 0) |
+| Main-chain params: 5 better 0 inside 1 worse|
 | Side-chain params: 5 better 0 inside 0 worse|
 |
|
*| Residue properties: Max.deviation: 4.0 Bad contacts: 1 |
*| Bond len/angle: 5.4 Morris et al class: 2 1 2 |
+| G-factors Dihedrals: -0.63 Covalent: -0.10 Overall: -0.38
 | M/c bond lengths: 89.7% within limits 10.3% highlighted
 | M/c bond angles: 98.0% within limits 2.0% highlighted

Thank you for your help,

Therese Malliavin
Unite de Bioinformatique Structurale
Institut Pasteur de Paris

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