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
Received on Tue Nov 29 2016 - 17:30:02 PST
