Brian,
An Amber 12 test run in a 30 A box seems working fine. I will do more
runs to make sure it works fine. Since Amber 11 runs fail while an Amber
12 run succeed, may I know if the Amber team has modified the related codes
in Amber 12 ?
Thanks,
Henry
On Thu, Aug 9, 2012 at 9:30 AM, Brian Radak <radak004.umn.edu> wrote:
> Henry,
>
> During QM/MM simulations with periodic boundaries, all atoms within qmcut
> (default is to match cut) of ANY qm atom are included in the direct space
> (and for QM/MM Lennard-Jones interactions). If that selection reaches
> outside the primary image, you will get the error you observed. However,
> the check is based on orthorhombic conditions and may fail unnecessarily in
> other boxes (Ross or Andreas, please correct me if this is wrong), but this
> is probably not a big deal in 90% of cases. You probably do just need a
> bigger box, which is recommendable for TI calculations that change charge
> too. Another alternative is to use a smaller qmcut. In terms of energy
> conservation 8 works just about as well as 10 when using qmmm_switch = 1,
> especially when no d-orbitals are present. I've also observed some strange
> simulations where an SCF instability causes the qm region to "explode" and
> trigger that error message.
>
> As for restraints, that could be a functional solution, although I'm not
> sure what you would restrain or if it would change the result of the TI. My
> general inclination is to avoid extra non-physically motivated energy terms
> whenever possible, but they can be pretty justifiable when the alternative
> is catastrophic failure of the simulation.
>
> Regards,
> Brian
>
> On Wed, Aug 8, 2012 at 6:02 PM, Pin-Chih Su (Henry Su) <psu4.uic.edu>
> wrote:
>
> > Brian,
> >
> > Thanks for your response. After I tried tempi= 100 or 150, the result
> is
> > identical. As you said, the energy is not conserved when enter the QMMM
> > production run and at the end of each stage, the energy will drop as
> > before. Please see the previous pictures as followed.
> >
> > *TI_energy<
> > http://i1076.photobucket.com/albums/w454/happypsu4/TI_energy.png>
> > (**http://i1076.photobucket.com/albums/w454/happypsu4/TI_energy.png**)*
> > *TI_pressure<
> > http://i1076.photobucket.com/albums/w454/happypsu4/TI_pressure.png>(
> > http://i1076.photobucket.com/albums/w454/happypsu4/TI_pressure.png)*
> > *TI_temperature<
> > http://i1076.photobucket.com/albums/w454/happypsu4/TI_temperature.png>(
> > http://i1076.photobucket.com/albums/w454/happypsu4/TI_temperature.png)*
> >
> > So now I am trying to finish the whole TI cycles and do the integral
> to
> > calculate the relative free energy between two ligands. However, when I
> am
> > running the ligand-water (no protein) run, the following error messages
> > pops out in each production run.
> >
> > ****************************************************
> > ERROR: QM region + cutoff larger than box dimension:
> > QM-MM Cutoff = 10.0000
> > Coord Lower Upper Size Radius of largest sphere inside
> unit
> > cell
> > X -12.960 12.051 25.011 14.723
> > Y -12.337 12.147 24.484 14.723
> > Z -14.754 14.694 29.448 14.723
> > ****************************************************
> > SANDER BOMB in subroutine QM_CHECK_PERIODIC<qm_mm.f>
> > QM region + cutoff larger than box
> > cannot continue, need larger box.
> >
> > I am aware there was the same issue in the other thread:
> > http://archive.ambermd.org/201110/0139.html
> >
> > I can see from the MD movies that the ligand deviate from the water box
> > center significantly and thus the run fails. Thus, I have tried a 30A
> water
> > box (the ligand itself is about 5A and the cutoff is 10A. So the box size
> > should be at least 10+10+5=25A) in Amber 11 but the same error message
> pops
> > out. I am now trying in Amber 12 again. If it doesn't work, I might try
> a
> > bigger box. But would you recommend if I restrain the ligand in a
> > "ligand-water only" QMMM production run with a weak force so that the
> > ligand would not flow away? Will the restraint sabotage the accuracy in
> a
> > TI run?
> >
> > The following is the proposed restraint script for convenience.
> >
> > ##################################################################
> > mask0=':MAA.F14,CL17'
> > mask1=
> >
> > for X in 1
> >
> > do
> >
> > ###################### minimization step ################################
> >
> > imin = 1,
> > maxcyc = 2000,
> > ntmin = 2,
> > ntpr = 100,
> > ntf = 2,
> > ntc = 1, (I also tried ntc=2 here)
> > ntb = 1,
> > cut = 10.0,
> > icfe = 1,
> > clambda = 0.${X},
> > ifsc=0,
> > crgmask='${mask0}',
> > &end
> > ###########################Constant volume equilibrium to equilibrate
> > temperature ######################################
> >
> > imin = 0,
> > irest = 0,
> > ntx = 1,
> > ntb = 1,
> > cut = 10.0,
> > ntr = 1,
> > ntc = 2,
> > ntf = 2,
> > tempi = 0.0,
> > temp0 = 300.0,
> > ntt = 3,
> > gamma_ln = 2.0,
> > nstlim = 50000,
> > dt = 0.001,
> > ntpr = 250,
> > ntwx = 250,
> > ntwr = 2500,
> > icfe = 1,
> > clambda = 0.${X},
> > ifsc=0,
> > crgmask='${mask0}',
> > */*
> > *Keep the ligand fixed with weak restraints*
> > *10.0*
> > *.1-23
> > *
> > */*
> >
> > END
> >
> > ##################################### Constant Pressure equilibrium to
> get
> > a proper density #########################
> > imin = 0,
> > irest = 1,
> > ntx = 7,
> > ntb = 2,
> > pres0 = 1.0,
> > ntp = 1,
> > taup = 2.0,
> > cut = 10,
> > ntr = 0,
> > ntc = 2,
> > ntf = 2,
> > tempi = 300.0,
> > temp0 = 300.0,
> > ntt = 3,
> > gamma_ln = 2.0,
> > nstlim = 50000,
> > dt = 0.001
> > ntpr = 250,
> > ntwx = 250,
> > ntwr = 2500,
> > icfe = 1,
> > clambda = 0.${X},
> > ifsc=0,
> > crgmask='${mask0}',
> > &end
> > ############################QMMM production
> > run###############################################
> > imin = 0,
> > ntx = 5,
> > irest = 1,
> > ntb = 2,
> > ntp = 1,
> > pres0 = 1.0,
> > taup = 2.0,
> > ntf = 1,
> > ntc = 1,
> > cut = 10.0,
> > temp0 = 300.0,
> > ntt = 3,
> > gamma_ln = 2.0,
> > nstlim = 50000,
> > dt = 0.001,
> > ntpr = 250,
> > ntwr = 2500,
> > ntwx = 250,
> > icfe=1,
> > clambda = 0.${X},
> > ifsc=0,
> > crgmask='${mask0}',
> > ifqnt = 1
> > /
> > &qmmm
> > iqmatoms= 1-23,
> > qmchage=-1,
> > qm_theory='PM3',
> > qmshake=0,
> > qm_ewald=1,
> > qm_pme=1,
> > * 10.0
> > .1-23*
> > /
> > &end
> >
> >
> >
> ############################################################################
> >
> >
> > Cheers,
> > Henry
> >
> >
> > On Tue, Aug 7, 2012 at 9:09 AM, Brian Radak <radak004.umn.edu> wrote:
> >
> > > Henry,
> > >
> > > I'm not an expert on ligand binding, but I believe the cycle in the
> > > tutorial will be compatible with QM/MM. The complication would be that
> > > there are different heats of formation for the ligands if they differ
> in
> > > the number and kinds of atoms. This will appear as a (probably unknown)
> > > constant shift in the *relative* free energy of transformation between
> > > ligands because the potential energies are on different scales.
> > >
> > > However, if you calculate two relative free energies (one in solution
> and
> > > one in the protein), then the shift will be the same in both legs and
> > > formally cancel out when you calculate a relative binding free energy.
> > You
> > > can do the math yourself by assuming a constant in each QM Hamiltonian
> > that
> > > is only a function of the kinds of QM atoms and not their coordinates.
> > The
> > > constants will come out of the configuration integrals and thus be
> > additive
> > > to the free energy.
> > >
> > > As for changing tempi, that may in fact be advantageous, but will not
> > > change the fact that energy is not conserved. That behavior should only
> > be
> > > desired/expected when ntt=0. Of course you will probably not get
> > physically
> > > meaningful results unless the system is first equilibrated near the
> > > temperature of interest (look up ensemble equivalence in the book by
> > Allen
> > > and Tildesley for example).
> > >
> > > Regards,
> > > Brian
> > >
> > > On Mon, Aug 6, 2012 at 6:51 PM, Pin-Chih Su (Henry Su) <psu4.uic.edu>
> > > wrote:
> > >
> > > > Hi Brian,
> > > >
> > > > Thanks for the response. I am trying tempi = 100 & 150 and will
> let
> > > you
> > > > know the result.
> > > >
> > > > Could you explain in more details about "how to carefully
> calculated
> > or
> > > > canceled out byan appropriate thermodynamic cycle" to overcome the
> > subtle
> > > > QMMM TI difficulties? If I follow and modify the protocol of amber TI
> > > > tutorial, <http://ambermd.org/tutorials/advanced/tutorial9/> will I
> > > > overcome the subtle QMMM TI difficulty? Or there are something
> > important
> > > > tips not mentioned in the tutorial?
> > > >
> > > > When visualized the QMMM TI trajectory, I don't see weird points.
> > > The
> > > > ligand still binds to the protein active site with normal
> fluctuation;
> > no
> > > > blow up. The following is the RMSD v.s. time plot of the 100 ps
> > > > production run in step 3, V0 group.
> > > >
> > > > http://i1076.photobucket.com/albums/w454/happypsu4/RMSD.jpg
> > > >
> > > > Here is my RMSD script
> > > >
> > > > ######################################
> > > > cat << EOF > AAA.calc_rms
> > > >
> > > > trajin BBB.mdcrd
> > > >
> > > > reference C.inpcrd
> > > >
> > > > rms reference mass out AAA.rms .N,CA,C time 0.5
> > > >
> > > > EOF
> > > >
> > > > $AMBERHOME/bin/ptraj C.prmtop < AAA.calc_rms
> > > > ###############################################
> > > >
> > > > Please let me know if any comment.
> > > >
> > > > Many thanks,
> > > > Henry
> > > >
> > > > On Fri, Aug 3, 2012 at 1:44 PM, Brian Radak <radak004.umn.edu>
> wrote:
> > > >
> > > > > Hi Henry,
> > > > >
> > > > > You appear to be running thermostatted MD (ntt != 1 and gamma_ln >
> > 0),
> > > in
> > > > > which case energy shouldn't be conserved, especially when you set
> > > tempi =
> > > > > 0., as the system quickly converts potential energy into kinetic
> > > energy.
> > > > > You might consider speeding that up by setting tempi > 0 or using
> > > various
> > > > > heating protocols.
> > > > >
> > > > > As for QM/MM, I believe it is quite normal for the energy to jump
> as
> > > the
> > > > QM
> > > > > code includes atomic heats of formation and therefore has a
> different
> > > > zero
> > > > > of potential energy. This is a subtle difficulty in performing TI
> > with
> > > > > QM/MM because the zero of energy will be an offset in any free
> > energies
> > > > > that you compute and must either be carefully calculated or
> canceled
> > > out
> > > > by
> > > > > an appropriate thermodynamic cycle.
> > > > >
> > > > > Are your dynamics unstable? Relatively large changes in energy
> during
> > > > > equilibration are not necessarily cause for concern unless you are
> > > > getting
> > > > > vlimit or SHAKE errors.
> > > > >
> > > > > Regards,
> > > > > Brian
> > > > >
> > > > > On Fri, Aug 3, 2012 at 1:44 PM, Pin-Chih Su (Henry Su) <
> psu4.uic.edu
> > >
> > > > > wrote:
> > > > >
> > > > > > Dear Amber,
> > > > > >
> > > > > > When I run a TI job, the energy doesn't seem to conserve in
> each
> > > > > clambda
> > > > > > step in no matter step 1, step 2 (soft core potential step) and
> > step
> > > 3.
> > > > > > The total energy always drops significantly at the beginning of
> > each
> > > > > stage
> > > > > > (minimization, constant volume equilibrium, constant
> > > > > > pressure equilibrium and production run). Moreover, when MD
> enters
> > > > into a
> > > > > > QMMM production run, the total energy will drop again and then
> > reach
> > > a
> > > > > > higher value than the total energy of minimization and
> equilibrium.
> > > > > Please
> > > > > > see the following Xmgrace figures for more details.
> > > > > >
> > > > > > *TI_energy<
> > > > > > http://i1076.photobucket.com/albums/w454/happypsu4/TI_energy.png
> >
> > > > > > (**
> > > > http://i1076.photobucket.com/albums/w454/happypsu4/TI_energy.png**)*
> > > > > > *TI_pressure<
> > > > > >
> http://i1076.photobucket.com/albums/w454/happypsu4/TI_pressure.png
> > >(
> > > > > >
> > http://i1076.photobucket.com/albums/w454/happypsu4/TI_pressure.png)*
> > > > > > *TI_temperature<
> > > > > >
> > > http://i1076.photobucket.com/albums/w454/happypsu4/TI_temperature.png
> > > > >(
> > > > > >
> > > >
> > http://i1076.photobucket.com/albums/w454/happypsu4/TI_temperature.png)*
> > > > > >
> > > > > >
> > > > > > The following is my Amber TI script in step 1 (to simplify, I
> > just
> > > > show
> > > > > > half of the script, as in a group ):
> > > > > >
> > > > > > mask0=':MAA.F14,CL17'
> > > > > > mask1=
> > > > > >
> > > > > > for X in 1,2,3
> > > > > >
> > > > > > do
> > > > > >
> > > > > > ###################### minimization step
> > > > ################################
> > > > > >
> > > > > > imin = 1,
> > > > > > maxcyc = 2000,
> > > > > > ntmin = 2,
> > > > > > ntpr = 100,
> > > > > > ntf = 2,
> > > > > > ntc = 1, (I also tried ntc=2 here)
> > > > > > ntb = 1,
> > > > > > cut = 10.0,
> > > > > > icfe = 1,
> > > > > > clambda = 0.${X},
> > > > > > ifsc=0,
> > > > > > crgmask='${mask0}',
> > > > > > &end
> > > > > > ###########################Constant volume equilibrium to
> > equilibrate
> > > > > > temperature ######################################
> > > > > >
> > > > > > imin = 0,
> > > > > > irest = 0,
> > > > > > ntx = 1,
> > > > > > ntb = 1,
> > > > > > cut = 10.0,
> > > > > > ntr = 1,
> > > > > > ntc = 2,
> > > > > > ntf = 2,
> > > > > > tempi = 0.0,
> > > > > > temp0 = 300.0,
> > > > > > ntt = 3,
> > > > > > gamma_ln = 2.0,
> > > > > > nstlim = 50000,
> > > > > > dt = 0.001,
> > > > > > ntpr = 250,
> > > > > > ntwx = 250,
> > > > > > ntwr = 2500,
> > > > > > icfe = 1,
> > > > > > clambda = 0.${X},
> > > > > > ifsc=0,
> > > > > > crgmask='${mask0}',
> > > > > > /
> > > > > > Keep the Protein fixed with weak restraints
> > > > > > 10.0
> > > > > > RES 1 258
> > > > > > /
> > > > > >
> > > > > > END
> > > > > >
> > > > > > ##################################### Constant Pressure
> equilibrium
> > > to
> > > > > get
> > > > > > a proper density #########################
> > > > > > imin = 0,
> > > > > > irest = 1,
> > > > > > ntx = 7,
> > > > > > ntb = 2,
> > > > > > pres0 = 1.0,
> > > > > > ntp = 1,
> > > > > > taup = 2.0,
> > > > > > cut = 10,
> > > > > > ntr = 0,
> > > > > > ntc = 2,
> > > > > > ntf = 2,
> > > > > > tempi = 300.0,
> > > > > > temp0 = 300.0,
> > > > > > ntt = 3,
> > > > > > gamma_ln = 2.0,
> > > > > > nstlim = 50000,
> > > > > > dt = 0.001
> > > > > > ntpr = 250,
> > > > > > ntwx = 250,
> > > > > > ntwr = 2500,
> > > > > > icfe = 1,
> > > > > > clambda = 0.${X},
> > > > > > ifsc=0,
> > > > > > crgmask='${mask0}',
> > > > > > &end
> > > > > > ############################QMMM production
> > > > > > run###############################################
> > > > > > imin = 0,
> > > > > > ntx = 5,
> > > > > > irest = 1,
> > > > > > ntb = 2,
> > > > > > ntp = 1,
> > > > > > pres0 = 1.0,
> > > > > > taup = 2.0,
> > > > > > ntf = 1,
> > > > > > ntc = 1,
> > > > > > cut = 10.0,
> > > > > > temp0 = 300.0,
> > > > > > ntt = 3,
> > > > > > gamma_ln = 2.0,
> > > > > > nstlim = 50000,
> > > > > > dt = 0.001,
> > > > > > ntpr = 250,
> > > > > > ntwr = 2500,
> > > > > > ntwx = 250,
> > > > > > icfe=1,
> > > > > > clambda = 0.${X},
> > > > > > ifsc=0,
> > > > > > crgmask='${mask0}',
> > > > > > ifqnt = 1
> > > > > > /
> > > > > > &qmmm
> > > > > > iqmatoms= 3868,3869,3870,3871,3872,
> > > > > > 3873,3874,3875,3876,3877,
> > > > > > 3878,3879,3880,3881,3882,
> > > > > > 3883,3884,3885,3886,3887,
> > > > > > 3888,3889,3890,1395,1394,
> > > > > > 1396,1398,1399,2328,2327,
> > > > > > 3833,3834,3835,3836,3837,
> > > > > > 3838,3839,3840,3841,3862,
> > > > > > 3863,3864,3865,3866,3867,
> > > > > > qmchage=-1,
> > > > > > qm_theory='PM3',
> > > > > > qmshake=0,
> > > > > > qm_ewald=1,
> > > > > > qm_pme=1,
> > > > > > /
> > > > > > &end
> > > > > >
> > > > > >
> > > > > >
> > > > >
> > > >
> > >
> >
> ############################################################################
> > > > > >
> > > > > > I have tried
> > > > > >
> > > > > > a. turn off QMMM. Or keep QMMM but with qmshake=1,. The same
> issue
> > > > > > b. Run in both Amber 11 & 12, The same issue
> > > > > > c. Turn on / off SHAKE. The same issue.
> > > > > > d. Turn off Langevin Dynamic. The same issue
> > > > > > e. Switch to Andersen temperature control. The same issue.
> > > > > > f. Prepare new .prmtop and .inpcrd files and run again. The same
> > > issue
> > > > > > g. I have set ntpr=1 and check all the .out files by "gvim -d"
> and
> > > > total
> > > > > > energies, temperature, pressure and all the MM components are the
> > > same
> > > > in
> > > > > > the same group file.
> > > > > > h. The restart files from V0 to V1 are identical.
> > > > > >
> > > > > > I only tried 9 windows so far.
> > > > > > clambda = 0.${X}, for X in 1 2 3 4 5 6 7 8 9
> > > > > > But I guess number of windows is not the issue here.
> > > > > >
> > > > > > Please let me know if any comment. Many comments.
> > > > > >
> > > > > > Henry
> > > > > > _______________________________________________
> > > > > > AMBER mailing list
> > > > > > AMBER.ambermd.org
> > > > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > > > >
> > > > >
> > > > >
> > > > >
> > > > > --
> > > > > ================================ Current Address
> > > =======================
> > > > > Brian Radak :
> > BioMaPS
> > > > > Institute for Quantitative Biology
> > > > > PhD candidate - York Research Group : Rutgers, The State
> > > > > University of New Jersey
> > > > > University of Minnesota - Twin Cities : Center for
> > > > Integrative
> > > > > Proteomics Room 308
> > > > > Graduate Program in Chemical Physics : 174 Frelinghuysen
> > Road,
> > > > > Department of Chemistry : Piscataway,
> > NJ
> > > > > 08854-8066
> > > > > radak004.umn.edu :
> > > > > radakb.biomaps.rutgers.edu
> > > > >
> ====================================================================
> > > > > Sorry for the multiple e-mail addresses, just use the institute
> > > > appropriate
> > > > > address.
> > > > > _______________________________________________
> > > > > AMBER mailing list
> > > > > AMBER.ambermd.org
> > > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > > >
> > > >
> > > >
> > > >
> > > > --
> > > >
> > > > Pin-Chih Su (Henry Su)
> > > >
> > > > Graduate Student
> > > >
> > > > Center for Pharmaceutical Biotechnology (MC 870)
> > > >
> > > > College of Pharmacy, University of Illinois at Chicago
> > > >
> > > > 900 South Ashland Avenue, Room 1052
> > > >
> > > > Chicago, IL 60607-7173
> > > >
> > > > office 312-996-5388
> > > >
> > > > fax 312-413-9303
> > > > _______________________________________________
> > > > AMBER mailing list
> > > > AMBER.ambermd.org
> > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > >
> > >
> > >
> > >
> > > --
> > > ================================ Current Address
> =======================
> > > Brian Radak : BioMaPS
> > > Institute for Quantitative Biology
> > > PhD candidate - York Research Group : Rutgers, The State
> > > University of New Jersey
> > > University of Minnesota - Twin Cities : Center for
> > Integrative
> > > Proteomics Room 308
> > > Graduate Program in Chemical Physics : 174 Frelinghuysen Road,
> > > Department of Chemistry : Piscataway, NJ
> > > 08854-8066
> > > radak004.umn.edu :
> > > radakb.biomaps.rutgers.edu
> > > ====================================================================
> > > Sorry for the multiple e-mail addresses, just use the institute
> > appropriate
> > > address.
> > > _______________________________________________
> > > AMBER mailing list
> > > AMBER.ambermd.org
> > > http://lists.ambermd.org/mailman/listinfo/amber
> > >
> >
> >
> >
> > --
> >
> > Pin-Chih Su (Henry Su)
> >
> > Graduate Student
> >
> > Center for Pharmaceutical Biotechnology (MC 870)
> >
> > College of Pharmacy, University of Illinois at Chicago
> >
> > 900 South Ashland Avenue, Room 1052
> >
> > Chicago, IL 60607-7173
> >
> > office 312-996-5388
> >
> > fax 312-413-9303
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
>
>
>
> --
> ================================ Current Address =======================
> Brian Radak : BioMaPS
> Institute for Quantitative Biology
> PhD candidate - York Research Group : Rutgers, The State
> University of New Jersey
> University of Minnesota - Twin Cities : Center for Integrative
> Proteomics Room 308
> Graduate Program in Chemical Physics : 174 Frelinghuysen Road,
> Department of Chemistry : Piscataway, NJ
> 08854-8066
> radak004.umn.edu :
> radakb.biomaps.rutgers.edu
> ====================================================================
> Sorry for the multiple e-mail addresses, just use the institute appropriate
> address.
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Aug 10 2012 - 00:30:03 PDT
