Dear Bharat,

for a synchrounous reaction, you can define a single RC,

if the reaction is well described by your QM hamiltonian, all other processes will follow barrierless during the steered MD. You may use the   AMBER LCOD module to run the reaction using a single meaningfull RC.

The definition of RC in LCOD is not using the restraints command, so use the definition given in the LCOD SMD manual/tutorial.

http://ambermd.org/tutorials/advanced/tutorial10/

ncsu_smd

output_file = 'smd.txt'

output_freq = 1000

variable

type = LCOD ! Linear Combination of Distances

i = (8,9,7,9)

r = (1.0,-1.0)

path = (0.7,-0.7)

harm = (1000.0)

end variable

end ncsu_smd

Alternativly you can scan the reaction profile using the Nudged elastic band NEB in AMBER.

We recently published a reaction with one hydride transfer and three proton transfers for a Zinc-and NADH dependent enzyme, and it worked well with a single RC. 

ACS Catal., 2015, 5 (6), pp 3207–3215

Dear Bharat,

You can refer to the AMBER manual for more details. Here I'd like simply explain some of your confuse.

First, r1, r2, r3 and r4 together define the harmonic potential that is used to "fix" your molecule at correct configuration along the RC. But the potential is ONLY non-zero within the range r1 to r2, and r3 to r4. That means, if the bond length drops in between r2 and r3, AMBER will not exert restraint to the bond.

Second, the tag '$rst iat' defines where the restraint you want to impose. The two numbers on the r.h.s. are the atomic indices of the bond atoms.

Third, if you want to visualize how an entire reaction occurs, please refer to Marco's recommendation, using steered MD or targeted MD (both supported by AMBER). US is not able to do so because it is a piece-wise or segmented sampling method, you can never witness a complete reaction within any windows of US.

Hope this helps,

Justin

Dear Prof. Bocola,

Thank you for your response. The reaction which I am studying occurs in a single step with a single transition state. The entire reaction happens like this:

1. Proton transfer from acid catalyst (Glu) to the leaving group O4 (O4 of the glycosidic bond to be broken)

2. breakage of C-O4 bond

3. Nucleophilic attack of water on anomeric C of C-O4 bond assisted by base catalyst Asp

Here's the link for the schematic view of the reaction: http://journals.plos.org/plosone/article/figure/image?size=large&id=info:doi/10.1371/journal.pone.0012947.g003

I am now confused which RC shall I choose for the reaction to proceed. Right now , I am trying SMD with the RC for step 1 i.e. proton transfer step. I will post you the results once simulation is finished to get some more advice on this topic. 

I read you paper which you referred to me and it seems three proton transfers were studied individually, whereas in my case reaction seems to be quite complex, as it involves simultaneous bond breakage and proton transfers (Ref. http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0012947). 

Could you please tell me how did you choose the value for reaction coordinate from +1.25 to  -1.25  for proton or hydride transfer? 

Thanks

BHARAT

Dear Dr. Jun Zhang, 

Thank you for the explanation. It seems it is not possible to define difference in distance restraints using the restraint file. I searched in amber manual and I found that this can be done using LCOD and ncsu_pmd approach.

I am able to break the bonds using restraints using usual restraint file and umbrella sampling, but I am not able to model the entire reaction using this approach. I tried to contact the authors of the paper who already published this work, but they haven't replied yet. Here is an excerpt from the their paper where they mention about using 4 different RC and umbrella sampling

"

The following simulation protocol was employed to calculate the reaction mechanism. The reaction coordinate which quantifies the glycosidic bond length, RC2, was used as a progress variable. An efficient iterative optimization procedure was used to obtain the optimized structures and minimum energy pathways along the reaction coordinates [40]. The minimum-energy pathway along RC2 subject to constraints on the other coordinates (RC1, RC3, RC4) traces out a pseudo-one dimensional stripe through a four-dimensional space. In the first step, the structure of the QM/MM system in each RC2 window was optimized with constrained RC1, RC3, and RC4 coordinates. An energetically reasonable path was found to be that in which the four RCs were incremented simultaneously.

Subsequently, umbrella sampling calculations [41] were performed along the reaction pathway in windows starting from the corresponding optimized structures. The Weighted Histogram Analysis Method (WHAM) [42] was used to unbias the umbrella-sampling results along the RC2 reaction coordinate, i.e. along the glycosidic bond cleavage. Harmonic restraints with force constants of 100 kcal/mol/Å for proton transfer (RC1, RC4) and 500 kcal/mol/Å for covalent bond breaking/formation (RC2, RC3) were used to guide the inverting reaction mechanism. 42 windows were computed, each with 10 ps of equilibration and 10 ps of production run."

I can't figure out how did they do this step? "In the first step, the structure of the QM/MM system in each RC2 window was optimized with constrained RC1, RC3, and RC4 coordinates. An energetically reasonable path was found to be that in which the four RCs were incremented simultaneously."

Regards

BHARAT