Targeted molecular dynamics (tMD) using GROMACS and PLUMED

Prepare your initial structure, topology file (e.g., .top), and parameter files for GROMACS as you would for a regular molecular dynamics simulation. I have used GROMACS 2018.2.

You can plumed from their official website https://www.plumed.org/doc-v2.8/user-doc/html/_installation.html . PLUMED patched with GROMACS provides enhanced sampling capabilities, including targeted molecular dynamics.

In tMD, you have to specify a set of collective variables that describe the progress of your targeted transition. These could be distances, angles, torsions, etc., that are relevant to your system's behaviour. In this example, I have used RMSD based CV.

Create a PLUMED input file (usually named "plumed.dat") where you define the collective variables and any restraints or biases you want to apply during the simulation. The PLUMED input file will have directives to specify the CVs, targets, and forces acting on the system.

My plumed.dat looks like this

UNITS LENGTH=A TIME=ps ENERGY=kcal/mol

rmsd: RMSD REFERENCE=target-BB.pdb TYPE=SIMPLE

restraint: ...

MOVINGRESTRAINT

ARG=rmsd

AT0=43.1 STEP0=0 KAPPA0=10

AT1=0.0 STEP1=90000000 KAPPA1=10

...

PRINT STRIDE=1000 ARG=* FILE=colar.dat

FIT_TO_TEMPLATE STRIDE=10 REFERENCE=reference-BB.pdb TYPE=SIMPLE

One thing to note is that the atoms IDs should exactly match between target and reference structure. For example, if working CA atoms, all the CA atoms should have the same atoms ID. Above file means that we are using simple moving restraints which are added to the reference structure. AT0 is the RMSD difference between the reference and target structure. We are asking it to be 0 (AT1) using a force constant of 10 kcal/mol (KAPPA0/KAPPA1) while completing the transition in 90000000 timesteps (total time = 90000000 * dt used in the mdp file). The results are printed every 1000 steps (PRINT STRIDE). FIT_TO_TEMPLATE STRIDE parameter affects how often code talks to PLUMED for calculating and applying the restraints. A stride of 1 will affect the performance a lot, in practice I have found that a value of 10 to 100 works best with a maximum 20-30% loss in the GROMACS performance.

Run the MD using: gmx_mpi mdrun -s topol.tpr -plumed plumed.dat

You can plot the RMSD change with respect to time using the generated covar.dat file. Analyze the trajectory for the progress of the transition and behaviour of the system. If any artefacts are found you can try changing the force constants and the number of transition steps.