Dear AMBER users,
I am performing solvation free energy calculation of ACE-ALA-ALA-NME or something like that in water by using the softcore potential TI and FEP by using AMBER 11 and 12.
The tolal solvation free energy calculations consist of three steps.
step1w/ elec part, V0 is no charge and V1 is full charge in TIP3P water.
step1w2/ elec part, the same as step1w but in vacuum.
step2w/ vdW part, softcore potential.
I used 21 lambda points for step1w and step1w2 (0.05 interval), and 49 points for step2w (0.02 interval). all 5 ns simulations.
But the output file from step2w/ included ELEC term for DV/DL.
DV/DL, AVERAGES OVER 2500000 STEPS
NSTEP = 2500000 TIME(PS) = 6000.000 TEMP(K) = 7.64 PRESS = 0.0
Etot = 23.4943 EKtot = 23.4943 EPtot = 7.8078
BOND = 0.0000 ANGLE = 0.0000 DIHED = 0.0000
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS = -6.9264
EELEC = 0.5888 EHBOND = 0.0000 RESTRAINT = 14.1454
EAMBER (non-restraint) = -6.3376
DV/DL = 7.8078
EKCMT = 0.0000 VIRIAL = 0.0000 VOLUME = 1.4092
Density = 0.0229
Ewald error estimate: 0.6058E-07
------------------------------------------------------------------------------
Could you please explain the origin of "EELEC=0.5888" above?
In my understanding, the solute has no charge by crgmask in this case.
Why can elec term have the values for DV/DL?
To clarify the condition of the simulations, I write the input files for step2w for V1 (ACE-ala-ala-NME in water system)
&cntrl
imin = 0, ntx = 7, irest = 1,
ntpr = 1000, ntwr = 100000, ntwx = 1000,
ntf = 1, ntc = 2,
ntb = 2, iwrap=1,
nstlim = 2500000, dt = 0.002,
temp0 = 300.0, ntt = 3, gamma_ln = 2,
ntp = 1, pres0 = 1.0, taup = 5.0,
icfe=1, clambda = 0.50,
ifmbar=1, bar_intervall=10,
bar_l_min=0.48, bar_l_max=0.52,
bar_l_incr=0.02, logdvdl=1,
ntr=1,
restraint_wt=10.0,
restraintmask=':1-4',
ifsc=1,
crgmask=':1-4',
scmask=':1-4',
&end
The input file for V0 is similar but crgmask and scmask is just '' and topology file of V0 is a water system with the same water coordinates as V1.
I used the restraints on all the solute atoms here to constrain the conformation to the extended one, but it is not the problem here.
I confirmed that free simulations also had nonzero values for EELEC DV/DB term from other simulations.
This is not specific for this system and toluene, methane, ethane, ala1 without constraints also had nonzero EELEC terms.
(It seems that we should use "EPtot" - " RESTRAINT" (not DV/DL) for the constrained solvation free energy calculations above.)
The input file of step1w was
&cntrl
imin = 0, ntx = 7, irest = 1,
ntpr = 1000, ntwr = 100000, ntwx = 1000,
ntf = 2, ntc = 2,
ntb = 2, iwrap=1,
nstlim = 2500000, dt = 0.002,
temp0 = 300.0, ntt = 3, gamma_ln = 2,
ntp = 1, pres0 = 1.0, taup = 5.0,
icfe=1, clambda = 0.50,
ifmbar=1, bar_intervall=10,
bar_l_min=0.45, bar_l_max=0.55,
bar_l_incr=0.05,
ntr=1,
restraint_wt=10.0,
restraintmask=':1-4',
ifsc=0,
crgmask='',
&end
and the output was
DV/DL, AVERAGES OVER 2500000 STEPS
NSTEP = 2500000 TIME(PS) = 6000.000 TEMP(K) = 0.00 PRESS = 0.0
Etot = 0.0000 EKtot = 0.0000 EPtot = -45.5407
BOND = 0.0000 ANGLE = 0.0000 DIHED = 0.0000
1-4 NB = 0.0000 1-4 EEL = 118.6249 VDWAALS = 0.0000
EELEC = -164.1657 EHBOND = 0.0000 RESTRAINT = 0.0000
DV/DL = -45.5407
Ewald error estimate: -0.5779E-05
------------------------------------------------------------------------------
In this elec case, there is no vdW term as expected.
and that of step1w2 was
&cntrl
imin = 0, ntx = 7, irest = 1,
ntpr = 1000, ntwr = 100000, ntwx = 1000,
ntf = 2, ntc = 2, nscm=1000,
ntb = 0, cut=999.0,
nstlim = 2500000, dt = 0.002,
temp0 = 300.0, ntt = 3, gamma_ln = 2,
icfe=1, clambda = 0.50,
ifmbar=1, bar_intervall=10,
bar_l_min=0.45, bar_l_max=0.55,
bar_l_incr=0.05,
ntr=1,
restraint_wt=10.0,
restraintmask=':1-4',
ifsc=0,
crgmask='',
&end
DV/DL, AVERAGES OVER 2500000 STEPS
NSTEP = 2500000 TIME(PS) = 6000.000 TEMP(K) = 0.00 PRESS = 0.0
Etot = 0.0000 EKtot = 0.0000 EPtot = -28.4696
BOND = 0.0000 ANGLE = 0.0000 DIHED = 0.0000
1-4 NB = 0.0000 1-4 EEL = 119.0847 VDWAALS = 0.0000
EELEC = -147.5543 EHBOND = 0.0000 RESTRAINT = 0.0000
DV/DL = -28.4696
------------------------------------------------------------------------------
As far as I understand, I needed to subtract step1w2 values from step1w to obtain the elec solvation free energy term correctly.
All these behaviors may be due to the architecture of amber solvation free energy codes.
Best regards,
Hironori Kokubo
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Sep 28 2012 - 04:00:02 PDT
