QM/MM Umbrella Sampling: How to generate initial conformations?

QM/MM for Biomolecular Simulation
1

Hi all,

I was fortunate enough to attend a recent CP2K/GROMACS course and am now trying to test out the interface on some other toy systems. The reaction I’m trying to model is abstraction of the hydrogen from a carbon by an oxygen attached to a ferryl atom (I’ve attached a picture for the sake of clarity with the distance (reaction co-ordinate) between the two highlighted). Where I’m getting stuck in the generation of the initial conformations to pass onto the umbrella sampling protocol

image .

Would I need to manually pull the hydrogen (i.e. change the co-ordinate) between the two or is there a more systematic way to do it? I appreciate for a dihedral-type problem I could use a strong pulling force in classical MD but this obviously won’t work for a bond-breaking type problem. Finally, for these types of calculations (where a bond involving hydrogen is being pulled), should I disable LINCS in GROMACS?

Many thanks for any help on this issue. Greatly appreciated.

Best,
Matthew

2

Hi Matthew,

I could suggest you the following workflow for performing that simulation:

  1. Take initial structure from classical MD (similar to the one you shown on the picture). Do QMMM equilibration for like ~10ps.

  2. After equilibration use pulling code to drive system from reagents to products, do so you could pull along that coordinate with some speed. I would suggest pulling over ~18 ps from 2.8 to 1 A:
    pull_coord1_rate = -0.01 ; pulling speed in nm/ps

  3. The you could extract frames from that pulling trajectory and use them as starting points for umbrella sampling simulations.

Some additional things which could be important:

  1. You have here Fe Ion in the QM system, check for suitable DFT functional, Basis set and pseudopotential, PW cut-off, which are suggested for iron in CP2K (they will probably be different from defaults). Also beware of multiplicity in QM system, it could be non 1.

  2. From my experience, for such kind of reactions better coordinate of the reaction would be not one distance, like here H - O distance (2.8 A), but rather difference in distances between O - H and C - H d(O-H) - d(C-H). Such simulation is also possible, but then you will need to know how to use Plumed in combination with GROMACS.

should I disable LINCS in GROMACS?

Yes.

Hope that ideas would help in your simulations

Best,
Dmitry