Dear Feng,
Not a problem at all, thank you for taking the time to answer me.
I’m a beginner on steered molecular dynamics so I guess I need to read more about it (if there is any literature you could suggest me, I’d be thankful). So by “reliable", according to my little understanding, I mean the total work achieved but also, yes, the free energy. I will carefully read the section of the manual you’ve pointed out. In fact, I would like to obtain the energy barrier of the reaction and I was thinking that I could get it with SMD.
The blowing water molecule didn’t occur in the example provided by the tutorial, it happened when I performed SMD for my protein and its ligand (178 residues in explicit solvent with Mg2+ and neutralising ions). But I had made a mistake on the steering path.. Now, I’ve re-run this calculation with the “appropriate” steering path, reducing the time step from 0.001 to 0.0005 and I haven’t seen any of the problems listed previously but the proton transfer doesn’t occur even if the spring constant is set to 1000 kcal/mol..
As the spring constant depends on the system, I guess I have to carry out multiple runs with different strengths. Is there a method to identify the most appropriate value for my system?
Thank you again for your time and any suggestion you can provide on this matter,
Sincerely,
Ruth
> On 29 Nov 2016, at 03:16, Feng Pan <fpan3.ncsu.edu> wrote:
>
> 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
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Nov 29 2016 - 08:00:02 PST
