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
Received on Wed Aug 08 2012 - 15:30:03 PDT
