Hi
Mahmood, you mention that there is a big jump in the value of your reaction
coordinate. Have you checked what happened in your system when the jump
happened. I have some experience with systems in which a "second reaction
happened", for example a rupture of a Hydroxyl group that is protonated to
give a water molecule. In those cases to improve things I do a 2D umbrella
sampling with both reaction coordinates, but only with the values of the
original reaction coordinate that are near the jump. But all of these is
valid if you see something in your system that can justify the jump.
Hope it helps.
2015-10-27 15:00 GMT+01:00 Brian Radak <brian.radak.accts.gmail.com>:
> In my experience, large force constants are not as necessary for bond
> difference coordinates as for simple bond coordinates. The final result
> of course depends on the shape of the underlying free energy surface,
> but I think 50-100 kcal/mol-A^2 is usually enough (that's an entirely
> opinion based estimate).
>
> I am sometimes skeptical of results using extremely large force
> constants (I've seen 500+), because those windows tend to have the
> /appearance/ of fast convergence when in fact there are several
> orthogonal degrees of freedom that are now slower and thus not converged.
>
> My reservations are often mitigated by using HREMD, in which case
> exchange amongst windows should either 1) even out the lack of
> convergence or 2) make any substantial errors very obvious.
> Unfortunately, only 1D-HREMD is usually affordable and feasible for
> QM/MM systems. The 0.1 A "rule of thumb" espoused by Adrian is also
> quite good for replica placement too in my experience.
>
> Regards,
> Brian
>
> On 10/27/2015 07:30 AM, Adrian Roitberg wrote:
> > Hi
> >
> > When we do this, for a simple bond breaking/forming, we tend to use
> > windows roughly 0.1 apart, and force constants of the order of 200.
> >
> > See like the overlap is better than before, but I would have more
> > windows between the ones you have now. Try also much larger force
> > constant around your transition state, which is probably around a value
> > of zero for your generalized coordinate.
> >
> > Adrian
> >
> >
> > On 10/27/15 8:05 AM, Mahmood Jasim wrote:
> >> Hi Adrian,
> >>
> >> Many thanks for your comments; they were very useful. Regarding ig==-1,
> I
> >> am using it except in the single input file I copied in the email. I am
> now
> >> running the simulation using the generalized distance coordinate as the
> >> difference between the lenghths of the bond to be formed and the one to
> be
> >> broken with 0.2 A windows. This is my restraint file:
> >>
> >> iat=1937,6803,6803,6798, r1=-10.0, r2=1.5, r3=1.5, r4=10.0, rstwt=1,-1,
> >> rk2=50, rk3=50,
> >>
> >> and the overlap has greatly improved. The problem starts to happen when
> the
> >> coordinate gets closer to the formation of the bond where a big gap
> appears
> >> in the graph and the coordinate jumps from 0.15 to -0.3 A. I tried
> >> increasing the force near the region but it did not work. This is how
> the
> >> reaction coordinate versus time looks like:
> >>
> >> [image: Inline images 2]
> >>
> >> Any advice please
> >>
> >> Best Regards,
> >> Mahmood
> >>
> >> On 23 October 2015 at 14:39, Adrian Roitberg <roitberg.ufl.edu> wrote:
> >>
> >>> Hi
> >>>
> >>> Several comments.
> >>>
> >>> Never run the WHOLE set of simulations before you are satisfied with
> the
> >>> overlap. Just run 3 or 4 windows and then see if you need to change
> >>> something before going further.
> >>>
> >>> As for your qeustions:
> >>>
> >>> Is that sufficient overlap? If it is not, do I need to reduce the
> spacing
> >>> or increase the force constant to get better overlap?
> >>>
> >>> You do not show the overlap, so it is hard to tell you, but my guess
> >>> from your first figure is that no, it is not enough. By the way, if you
> >>> want to increase overlap, you need either a closer spacing and/or a
> >>> LOWER force constant(not increased!).
> >>>
> >>> PLEASE add ig=-1 to your mdin !
> >>>
> >>> As for your coordinate, it is NEVER a good idea to use just one
> distance
> >>> for a reaction. Look at published work on this, and you will see that
> >>> peopel tend to use a reaction coordinate that has RC=distance(bond
> being
> >>> broken) - distance(bond being made)
> >>>
> >>> That produced much more stable surfaces.
> >>>
> >>> adrian
> >>>
> >>>
> >>>
> >>> On 10/23/15 7:56 AM, Mahmood Jasim wrote:
> >>>> Hi AMBER users,
> >>>>
> >>>> I am trying to simulate a reaction involving the formation of a
> covalent
> >>>> bond using umbrella sampling combine with QM treatment of the reaction
> >>>> atoms. I have the SG atom of a cysteine residue attacking an
> electrophilc
> >>>> carbon on a ligand. The distance between these 2 atoms is the reaction
> >>>> coordinate. The windows of umbrella sampling are spaced 0.3 Angstrom
> >>> apart.
> >>>> The restraint file looks like this:
> >>>>
> >>>> #
> >>>> # 123 CYM SG 432 UNK C1 3.43
> >>>> &rst
> >>>> ixpk= 0, nxpk= 0, iat=1937,6803, r1= 2.93, r2= 3.43, r3= 3.43,
> r4=
> >>> 3.93,
> >>>> rk2=0.0, rk3=100.0, ir6=1, ialtd=0,
> >>>> &end
> >>>>
> >>>> and the md.in file looks like this:
> >>>>
> >>>> &cntrl
> >>>> imin = 0,
> >>>> irest = 1,
> >>>> ntx = 7,
> >>>> ntb = 2,
> >>>> cut = 12,
> >>>> ntr = 0,
> >>>> ntc = 1,
> >>>> ntf = 1,
> >>>> igb = 0
> >>>> ntp = 1
> >>>> tempi = 300.0,
> >>>> temp0 = 300.0,
> >>>> ntt = 3,
> >>>> gamma_ln = 1.0,
> >>>> nstlim =10000, dt = 0.001
> >>>> ntpr = 1000, ntwx = 1000, ntwr = 1000, nmropt = 1,
> >>>> ifqnt=1,
> >>>> /
> >>>>
> >>>>
> >>>> &qmmm
> >>>>
> >>>> qmmask =
> >>>> '.1934-1937,6794-6848',
> >>>> qmcharge =
> >>>> 0,
> >>>> qm_theory =
> >>>> 'PM3',
> >>>> qmcut =
> >>>> 12.0,
> >>>> /
> >>>> &wt type='DUMPFREQ', istep1=100,/
> >>>> &wt type='END' /
> >>>> DISANG=RST-1.dist
> >>>> LISTOUT=POUT
> >>>> DUMPAVE=DIST_1.dat
> >>>>
> >>>> The run was able to simulate the formation of the bond at the end and
> the
> >>>> next graph represents the values of the coordinates across the
> windows.
> >>>> [image: Inline images 1]
> >>>>
> >>>> Is that sufficient overlap? If it is not, do I need to reduce the
> spacing
> >>>> or increase the force constant to get better overlap?
> >>>>
> >>>> I used the WHAM code by Alan Grossfieldto calculate PMF with this
> >>> command:
> >>>> wham 1.59 3.73 70 0.00001 300 0 meta.dat result.dat
> >>>>
> >>>> and this was the meta.dat file:
> >>>>
> >>>> DIST_1.dat 3.43 200
> >>>> DIST_2.dat 3.13 200
> >>>> DIST_3.dat 2.83 200
> >>>> DIST_4.dat 2.53 200
> >>>> DIST_5.dat 2.23 200
> >>>> DIST_5.dat 1.93 200
> >>>> DIST_7.dat 1.63 200
> >>>>
> >>>> and the histogram I am getting looks like:
> >>>>
> >>>> [image: Inline images 2]
> >>>>
> >>>> I was expecting the histogram to look different, the PMF to go higher
> and
> >>>> then decreases as it approcahes the products. If I managed to improve
> the
> >>>> overlapping, will this change the histogrm?
> >>>>
> >>>> Many thanks,
> >>>> Mahmood Jasim
> >>>> Aston University
> >>>> UK
> >>>>
> >>>>
> >>>>
> >>>> _______________________________________________
> >>>> AMBER mailing list
> >>>> AMBER.ambermd.org
> >>>> http://lists.ambermd.org/mailman/listinfo/amber
> >>> --
> >>> Dr. Adrian E. Roitberg
> >>> Professor.
> >>> Department of Chemistry
> >>> University of Florida
> >>> roitberg.ufl.edu
> >>> 352-392-6972
> >>>
> >>> _______________________________________________
> >>> 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
>
> --
> Brian Radak
> Postdoctoral Scholar
> Gordon Center for Integrative Science, W323A
> Department of Biochemistry & Molecular Biology
> University of Chicago
> 929 E. 57th St.
> Chicago, IL 60637-1454
> Tel: 773/834-2812
> email: radak.uchicago.edu
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Dr. Emilio Lence Quintana
Organic Chemistry Department-CIQUS
University of Santiago de Compostela
15782 Santiago de Compostela, Spain
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Oct 27 2015 - 09:30:05 PDT
