Dear Feng,
Thank you a lot for your helpful answer.
I’ve plotted the CV vs work and the CV doesn’t follow the steering path, so I’ll first try a longer simulation time.
Thank you again for your time,
Sincerely,
Ruth
> On 30 Nov 2016, at 02:22, Feng Pan <fpan3.ncsu.edu> wrote:
>
> Hi, Ruth
>
> Based on the situation that you didn't observe the proton transfer in your
> own system. What I suggest is:
>
> - Check the SMD output file to see if the CV is well steered along you
> path. If not, maybe the strength is not enough or your simulation time is
> too short.
> - If the steering is OK, but no transfer happens, you analyze your
> trajectories using cpptraj or visualize it with VMD to check why the CV is
> steered but proton is not transferred.
> In this case, always you CV (here LCOD?) can not describe the transition
> of states well, so maybe new CV needs to be chosen.
>
> I am not familiar with proton transfer in protein but I hope this may be
> helpful to you
>
> Best
> Feng
>
> On Tue, Nov 29, 2016 at 10:59 AM, Ruth Helena Tichauer <rhtichau.laas.fr>
> wrote:
>
>> 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
>>
>
>
>
> --
> 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 Wed Nov 30 2016 - 02:00:03 PST
