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
Received on Fri Aug 03 2012 - 12:00:04 PDT
