AMBER Archive (2008)

Subject: AMBER: Force fields for REMD with GB solvent

From: guardiani_at_fi.infn.it
Date: Tue Jul 08 2008 - 11:24:27 CDT

Dear Amber experts,

I am currently performing REMD simulations in implicit solvent of
a peptide whose coordinates were taken from a larger protein.
My goal is to determine the equilibrium structure of the free
peptide in physiologic conditions. The peptide is first minimized
according to:

Minimization
  &cntrl
   imin = 1,
   maxcyc = 500,
   ncyc = 250,
   ntb = 0,
   igb = 1,
   gbsa = 1,
   saltcon = 0.2,
   cut = 300.0,
   ntpr = 10,
   ntx = 1,
  /

The minimized structure is then cloned into 12 replicas (12 corresponds
to the square root of the number of atoms in my peptide). The target
temperatures of the replicas are chosen in geometric progression and
are: 310, 332, 355, 380, 406, 435, 465, 498, 533, 570, 610, 652.
Each replica is heated from 0 K to the corresponding target temperature
in this way:

Equilibration
  &cntrl
   imin = 0,
   cut = 24.0,
   igb = 1,
   gbsa = 1,
   saltcon = 0.2,
   nstlim = 500000,
   dt = 0.001,
   ntt = 3, gamma_ln = 1.0,
   tempi = 0.0, temp0 = 310.0,
   ntx = 1, irest = 0, ntb = 0,
   nscm = 50,
   ntpr = 100, ntwr = 100, ntwx = 100
  /

I then perform 2.5 ns of REMD simulation that will be discarded in the
computation of energies and other properties of interest. The input
file for this part of the simulation is like this:

MD run
  &cntrl
   imin = 0,
   cut = 24.0,
   igb = 1,
   gbsa = 1,
   saltcon = 0.2,
   nstlim = 2500,
   numexchg = 800,
   dt = 0.001,
   ntt = 3, gamma_ln = 1.0,
   tempi = 310.0, temp0 = 310.0,
   ntx = 5, irest = 1, ntb = 0,
   nscm = 50,
   repcrd = 0,
   ntpr = 500, ntwr = 500, ntwx = 500
  /

After that I perform other 30 ns of simulation through 3 restarts with
input files of this kind:

MD run
  &cntrl
   imin = 0,
   cut = 24.0,
   igb = 1,
   gbsa = 1,
   saltcon = 0.2,
   nstlim = 2500,
   numexchg = 4000,
   dt = 0.001,
   ntt = 3, gamma_ln = 1.0,
   tempi = 310.0, temp0 = 310.0,
   ntx = 5, irest = 1, ntb = 0,
   nscm = 50,
   repcrd = 0,
   ntpr = 500, ntwr = 500, ntwx = 500
  /

My problem is that the equilibrium population is completely dominated by
helical structures, which is in disagreement with the simulation of some
colleague who is performin explicit-solvent simulations on the same peptide.

Basically I think there are two main sources of error:

1) I used the ff99 force field without noticing that it introduces a bias
    toward helical structures. I am now not sure which force field is more
    convenient for my simulations. Do you recommend I should use the ff03
    force field or rather the ff99SB force field, or maybe yet another one ?

2) I performed my simulations with the Hawkins, Cramer, Trular generalized
    Born model (igb=1). However I have read a paper
    (J. Chem. Theory Comput. Vol 3, No 1, pag 156-169, 2007)
    where the authors report on a comparison of REMD simulations on Ala10
    and they conclude that OBC-GB and HCT-GB introduce a strong bias for
    alpha-helical conformations while the GBn method (igb=7) produces a
    conformational population in better agreement with experimental data.
    Do you think the GBn method is the best choice of GB model for my
    simulations ? And which set of atomic radii (PBradii) do you recommend ?

I thank you very much for your help and I am looking forward to receiving
feedback from you.

Best regards,

               Carlo Guardiani

PS: If you notice any mistake in the REMD protocol that I am using, please
     tell me.

----------------------------------------------------------------
This message was sent using IMP, the Internet Messaging Program.

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
      to majordomo_at_scripps.edu