Re: [AMBER] Optimal values for harm and output_freq in SMD

From: Feng Pan <fpan3.ncsu.edu>
Date: Mon, 28 Nov 2016 21:16:18 -0500

Hi, Ruth

Sorry for replying late because of Thanksgiving break

>As I haven’t set the random seed I guess it was different for each run..
On this matter, when doing smd for my protein and its ligand, should I run
the simulation with a random seed? If yes, how many replicas should I do in
order to obtain “reliable” results?

Here what "reliable" results do you mean? If you mean the free energy by
got from the work, you should use the Jarzynski relationship equation, the
section 22.7 in manual has a pretty good explanation of this.

>And, which would be the “ideal” value for the string constant in order not
to have such huge differences in the total work? I’ve tried the example
with harm=100 kcal/mol and the CV vs work plot is not as bad as for harm=10
kcal/mol, but is 100 kcal/mol a reasonable value (not too high)? Also, I’ve
just got the output files of a SMD of my protein performed with harm= 1000
kcal/mol and:
_During many steps the Self Consistency is not achieved
_The water molecule involved in a proton transfer I wish to observe
completely blows as one of its hydrogens is 18 A° away from the oxygen atom
and yet the proton transfer doesn’t occur.. But, I’should confirm this
results again as I think I’ve made a mistake on the steering path..

The steering strength really depends on the system you perform, here I also
think 1000 kcal/mol is too large. Generally if you run a longer SMD you
could use a smaller strength constant and the steering will be smoother.
Here I suggest if it is not too expensive, you can run for a longer
simulation time and using a smaller strength like 20 or 50 to see if it is
good.
I have not tried the example in the tutorial so I don't have a clue why
your water molecule blows up, probably the steering path problem. I may try
it when I have time this week to see what happens.

Best
Feng

On Fri, Nov 25, 2016 at 6:47 AM, Ruth Helena Tichauer <rhtichau.laas.fr>
wrote:

> Dear Feng,
>
> Indeed, I’m running the example and I was surprised to observe such huge
> differences only by changing the output_freq. Here are the input files and
> a part of the output files referring to the NFE section for:
>
> _harm=10 kcal/mol and output_freq=1000:
>
> Malonaldehyde: NCSU steered MD: LCOD rxn coord
> &cntrl
> imin = 0,
> irest = 1,
> ntx = 5,
> ntb = 0,
> cut = 999.0,
> tempi = 300.0,
> temp0 = 300.0,
> ntt = 3,
> gamma_ln = 1.0,
> ntf = 2, ntc = 2, tol = 0.00001,
> dt = 0.0005,
> ntpr = 50, ntwr=100, ntwx = 100,
> nstlim = 100000,
> ifqnt = 1, infe=1
> /
> &qmmm
> qmmask=':*',
> qmcharge=0,
> qm_theory='DFTB',
> qmshake=0,
> writepdb=1,
> /
> &smd
> output_file = 'smd_10.txt'
> output_freq = 1000
> cv_file='cv.in'
> /
>
> NFE : output_file = smd_10.txt
> NFE : output_freq = 1000 (0.5000 ps)
> NFE :
> NFE : CV #1
> NFE :
> NFE : <> path = (0.7000, -0.7000)
> NFE : <> path_mode = SPLINE
> NFE : <> harm = (10.0000)
> NFE : <> harm_mode = SPLINE
> NFE :
> NFE : type = 'LCOD' (Linear Combination Of Distances)
> NFE : 1.000 * (8 [O2] <=> 9 [H4])
> NFE : -1.000 * (7 [O1] <=> 9 [H4])
>
> _harm=10 kcal/mol and output_freq=500:
>
> Malonaldehyde: NCSU steered MD: LCOD rxn coord
> &cntrl
> imin = 0,
> irest = 1,
> ntx = 5,
> ntb = 0,
> cut = 999.0,
> tempi = 300.0,
> temp0 = 300.0,
> ntt = 3,
> gamma_ln = 1.0,
> ntf = 2, ntc = 2, tol = 0.00001,
> dt = 0.0005,
> ntpr = 50, ntwr=100, ntwx = 100,
> nstlim = 100000,
> ifqnt = 1, infe=1
> /
> &qmmm
> qmmask=':*',
> qmcharge=0,
> qm_theory='PM3',
> qmshake=0,
> writepdb=1,
> /
> &smd
> output_file = 'smd_10-500.txt'
> output_freq = 500
> cv_file='cv.in'
> /
>
> NFE : output_file = smd_10-500.txt
> NFE : output_freq = 500 (0.2500 ps)
> NFE :
> NFE : CV #1
> NFE :
> NFE : <> path = (0.7000, -0.7000)
> NFE : <> path_mode = SPLINE
> NFE : <> harm = (10.0000)
> NFE : <> harm_mode = SPLINE
> NFE :
> NFE : type = 'LCOD' (Linear Combination Of Distances)
> NFE : 1.000 * (8 [O2] <=> 9 [H4])
> NFE : -1.000 * (7 [O1] <=> 9 [H4])
>
> I’ve plotted the CV vs work and as you point out the CV doesn’t follow
> very well the steering path when harm=10 kcal/mol.. That is when I started
> changing other parameters.
>
> As I haven’t set the random seed I guess it was different for each run..
> On this matter, when doing smd for my protein and its ligand, should I run
> the simulation with a random seed? If yes, how many replicas should I do in
> order to obtain “reliable” results?
>
> And, which would be the “ideal” value for the string constant in order not
> to have such huge differences in the total work? I’ve tried the example
> with harm=100 kcal/mol and the CV vs work plot is not as bad as for harm=10
> kcal/mol, but is 100 kcal/mol a reasonable value (not too high)? Also, I’ve
> just got the output files of a SMD of my protein performed with harm= 1000
> kcal/mol and:
> _During many steps the Self Consistency is not achieved
> _The water molecule involved in a proton transfer I wish to observe
> completely blows as one of its hydrogens is 18 A° away from the oxygen atom
> and yet the proton transfer doesn’t occur.. But, I’should confirm this
> results again as I think I’ve made a mistake on the steering path..
>
> I really appreciate the time you spend in order to enlighten me on this
> subject.
> Thank you for any forward suggestion/enlightenment,
>
> Sincerely,
>
> Ruth
>
>
> > On 24 Nov 2016, at 19:55, Feng Pan <fpan3.ncsu.edu> wrote:
> >
> > Hi, Ruth
> >
> > The output_freq only affects the output so it should not affect the total
> > work. Are you sure you ran both cases with exact same conditions? The
> > random seed may also affect the
> > results of SMD.
> > Are you running the example? I think the harm strength 10kcal/mol is too
> > low because the CV path does not follow the Steering path very well. Try
> a
> > larger harm strength, could be better.
> > If the strength is too low, the steering does not take many effects, so
> the
> > total work may differ a lot from different runs.
> >
> > Feng
> >
> > On Thu, Nov 24, 2016 at 6:34 AM, Ruth Helena Tichauer <rhtichau.laas.fr>
> > wrote:
> >
> >> Dear Amber users,
> >>
> >> I’ve been following tutorial A10 in order to perform SMD for my protein
> >> and its ligand using the LCOD method. In the tutorial, the spring
> constant
> >> for the proton transfer to occur is set to 1000 kcal/mol whose value,
> >> according to a few related questions on the mailing list, is too high.
> >>
> >> In order to find a reasonable value, I’ve run the simulation several
> times
> >> and the proton transfer takes place for “low” values such as 10
> kcal/mol if
> >> the output_freq is set to 1000, like in the tutorial. But, if the
> >> output_freq is lowered (I’ve tried 500 and 50) the proton transfer
> doesn’t
> >> occur anymore (for harm=10kcal/mol).
> >>
> >> Moreover, in the output_file .txt the CV, handle_position and work are
> not
> >> the same for a given time during the MD, resulting in different values
> for
> >> the total work done.
> >>
> >> For harm=10 kcal/mol and output_freq=1000 I obtain:
> >>
> >> #
> >> # MD time (ps), CV, handle_position, spring_constant, work
> >> #
> >> 925.0000 0.50139387 0.70000000 10.00000000
> >> 0.00000000
> >> 925.5000 0.89660018 0.69958280 10.00000000
> >> 0.00030631
> >> 926.0000 0.83105576 0.69834240 10.00000000
> >> 0.00081418
> >> 926.5000 0.66819394 0.69629560 10.00000000
> >> 0.00204283
> >> 927.0000 0.60709575 0.69345920 10.00000000
> >> 0.00467988
> >> 927.5000 0.93546672 0.68985000 10.00000000
> >> 0.00719652
> >> 928.0000 0.56761201 0.68548480 10.00000000
> >> 0.01073920
> >> 928.5000 0.62737834 0.68038040 10.00000000
> >> 0.01357336
> >> 929.0000 0.92623955 0.67455360 10.00000000
> >> 0.01711340
> >> 929.5000 0.66861857 0.66802120 10.00000000
> >> 0.02039487
> >> 930.0000 0.91611501 0.66080000 10.00000000
> >> 0.02901486
> >> ...
> >> #
> >> # <> total work done: 0.7690942380
> >>
> >> For harm=10 kcal/mol and output_freq=500 I obtain:
> >>
> >> #
> >> # MD time (ps), CV, handle_position, spring_constant, work
> >> #
> >> 925.0000 0.50139387 0.70000000 10.00000000
> >> 0.00000000
> >> 925.2500 0.80711710 0.69989535 10.00000000
> >> 0.00022605
> >> 925.5000 1.00173702 0.69958280 10.00000000
> >> 0.00094337
> >> 925.7500 0.81672099 0.69906445 10.00000000
> >> 0.00198729
> >> 926.0000 0.81074280 0.69834240 10.00000000
> >> 0.00416796
> >> 926.2500 0.93641120 0.69741875 10.00000000
> >> 0.00658384
> >> 926.5000 0.79972556 0.69629560 10.00000000
> >> 0.00872512
> >> 926.7500 0.88288670 0.69497505 10.00000000
> >> 0.01112337
> >> 927.0000 0.99659450 0.69345920 10.00000000
> >> 0.01380751
> >> 927.2500 0.91748709 0.69175015 10.00000000
> >> 0.01664516
> >> 927.5000 0.98411987 0.68985000 10.00000000
> >> 0.01972750
> >> 927.7500 0.76115131 0.68776085 10.00000000
> >> 0.02316773
> >> 928.0000 0.82942938 0.68548480 10.00000000
> >> 0.02737620
> >> 928.2500 0.74668613 0.68302395 10.00000000
> >> 0.03189408
> >> 928.5000 0.95082231 0.68038040 10.00000000
> >> 0.03673629
> >> 928.7500 0.79738517 0.67755625 10.00000000
> >> 0.04177707
> >> 929.0000 1.03951245 0.67455360 10.00000000
> >> 0.04712890
> >> 929.2500 0.71139449 0.67137455 10.00000000
> >> 0.05324293
> >> 929.5000 0.85144838 0.66802120 10.00000000
> >> 0.05961067
> >> 929.7500 0.92964077 0.66449565 10.00000000
> >> 0.06619512
> >> 930.0000 0.84766718 0.66080000 10.00000000
> >> 0.07329916
> >> ...
> >> #
> >> # <> total work done: 11.2413859667
> >>
> >> So I wonder what exactly the output_freq achieves during the SMD? As I
> >> want to run the same kind of simulation for my protein and its ligand,
> I’d
> >> like to know how to choose this value.
> >>
> >> Thank you for any insight on this subject,
> >>
> >> Sincerely,
> >>
> >> Ruth
> >> _______________________________________________
> >> AMBER mailing list
> >> AMBER.ambermd.org
> >> http://lists.ambermd.org/mailman/listinfo/amber
> >>
> >
> >
> >
> > --
> > Feng Pan
> > Ph.D. Candidate
> > North Carolina State University
> > Department of Physics
> > Email: fpan3.ncsu.edu
> > _______________________________________________
> > 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
>



-- 
Feng Pan
Ph.D. Candidate
North Carolina State University
Department of Physics
Email:  fpan3.ncsu.edu
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Nov 28 2016 - 18:30:02 PST
Custom Search