Version：Gromacs2022.5 cp2k2022.1 Plumed2.8.2

##############################################################################

actin.gro and actin.top files are from actin.inpcrd and actin.prmtop, and the charges are normal amber charges and LJ parameters are 0 for HW (hydrogen of water) and OH (hydrogen of Ser, Thr and Tyr).

Note: do not change the LJ parameters when you use gromacs combined with cp2k, but one shuold adjust the charges for the partial residues caused by truncation.

One must change the charges and LJ parameters if you use only cp2k to run qm/mm/md.

##############################################################################

see https://www.cp2k.org/howto:biochem_qmmm

Note the FORCE_EVAL/QMMM#LINKsubsection. This tells cp2k where the two subsystems are linked, in this case through a bond between MM atom 1411 (Cα) and QM atom 1413 (Cβ), as well as how to treat the link. We are using the IMOMM method for the link. The arginine sidechain has now been moved to the QM region. This means we need to adjust the QM charge to -1, but this also affects the MM region. The partial charges on the atoms of the arginine residue sum up to +1. When we move the sidechain atoms to the QM region, their charges are no longer counted in the MM region, leading to a residual partial charge. In our case, the total charge amounts to -0.0362. We can neutralize the system by distributing an equal but opposite charge across all 6 mainchain atoms of the residue. For example, the charge on the N is -0.3479. Subtracting -0.0362/6 from that gives us the new charge, -0.3419. Do the same calculation for the 5 remaining mainchain atoms and check that the modified charges for the six atoms add up to zero. Then, apply the changes using parmed’s change command:

change charge @1409 -0.3419

Before running the QM/MM simulation we need to amend our prmtop file, because in the classical forcefield, several hydrogen atom types do not have separate Lennard-Jones parameters or they are set up to 0.0. The AMBER FF atom types are HO (from the TIP3P water model), HG (from serine, threonine and tyrosine residues). This will cause unphysical interactions with the QM region and the simulation will fail due for stability reasons. The hydrogen atoms with these atom types will strongly interact with the QM electron clouds, eventually crashing the system. To prevent this from happening, we have to add the correct LJ parameters to the mentioned atom types using the ambertools utility parmed.

Parmed reads the prmtop file and allows you to modify the parameters defined in it. The parmed command we are going to use is: changeLJSingleType, which changes the LJ type for all of the atoms that share the LJ type, or addLJType, for more selective modifications. We will change their LJ parameters to those of a GAFF2 type alcohol that has 0.3019 0.047 values for the LJ type.

Either if we are using changeLJSingleType or addLJType, we need to ensure that we are changing the right atom index. To do so, we need to use printDetails or printLJTypes that print the information on the selected atom index.

$ parmed complex.prmtop
changeLJSingleType :WAT@H1 0.3019 0.047
changeLJSingleType :*@HO 0.3019 0.047
outparm complex_LJ_mod.prmtop
quit

##############################################################################

1. gmx_mpi make_ndx -f actin.gro

adding QMatomics group in the index.ndx file

[ QMatoms ]

5828 5829 5830 5831 5832 5833 5834 5835 5836 5837 5838 5839 5840 5841 5842

5843 5844 5845 5846 5847 5848 5849 5850 5851 5852 5853 5854 5855 5856 5857

5858 5859 5860 5861 5862 5863 5864 5865 5866 5867 5868 5869 5870

5871 6022 6023 6024 6034 6035 6036 6058 6059 6060 6118 6119 6120

2441 2442 2443 2444 2445 2446 2447 2448 2449 2450 2451

2362 2363 2364 2365 2366 2367

2108 2109 2110 2111 2112 2113 2114 2115 2116 2117 2118

6256 6257 6258 5992 5993 5994

2. gmx_mpi grompp -f em-qmmm.mdp -p actin.top -c actin.gro -n index.ndx -o em-qmmm.tpr

em-qmmm.mdp

title = QMMM energy minimization

integrator = steep

tinit = 0

dt = 0.002

nsteps = 1000

nstcomm = 1

comm_grps = system

emtol = 100.0

emstep = 0.01

nstxout = 1

nstvout = 1

nstfout = 1

nstlog = 1

nstenergy = 1

nstxout-compressed = 1

compressed-x-precision = 1000

compressed-x-grps = system

energygrps =

nstlist = 10

ns_type = grid

pbc = xyz

rlist = 1.0

coulombtype = PME

rcoulomb-switch = 1.0

rcoulomb = 1.6

vdwtype = Cut-off

rvdw-switch = 0

rvdw = 1.6

DispCorr = No

constraints = none

qmmm-cp2k-active = true

qmmm-cp2k-qmgroup = QMatoms

qmmm-cp2k-qmmethod = PBE

qmmm-cp2k-qmcharge = -3

qmmm-cp2k-qmmultiplicity = 1

3. mpirun gmx_mpi mdrun -deffnm em-qmmm

For Mg²⁺ ion,

error The requested basis set <DZVP-MOLOPT-GTH> for element <MG> was not found in the basis set files <BASIS_MOLOPT>

For cp2k input file, one may need to change some parameters.

&FORCE_EVAL

METHOD QMMM

&PRINT

&FORCES ON

&END FORCES

&END PRINT

to print the force of each atom

and change DZVP-MOLOPT-GTH to DZVP-MOLOPT-SR-GTH, and

qmmm-cp2k-active = true

qmmm-cp2k-qmgroup = QMatoms

qmmm-cp2k-qmmethod = PBE

qmmm-cp2k-qmcharge = -3

qmmm-cp2k-qmmultiplicity = 1

qmmm-cp2k-active = true

qmmm-cp2k-qmgroup = QMatoms

qmmm-cp2k-qmmethod = INPUT

4. gmx_mpi grompp -f em-qmmm.mdp -qmi em-qmmm.inp -p actin.top -c actin.gro -n index.ndx -o em-qmmm.tpr

5. sbatch sub.sh

#!/bin/bash

#SBATCH -p v6_384

#SBATCH -N 1

#SBATCH -c 4

#SBATCH -n 24

source /public5/home/sch7642/soft/cp2k-gromacs2022/cp2k-2022.1/tools/toolchain/install/setup

export PATH=/public5/home/sch7642/soft/cp2k-gromacs2022/gromacs-2022.5/install/bin:$PATH

export PATH=/public5/home/sch7642/soft/cp2k-gromacs2022/cp2k-2022.1/exe/local:$PATH

gmx_mpi mdrun -deffnm em-qmmm

6. run metadynamics. (srun gmx_mpi mdrun -deffnm Qmmm-metad -plumed plumed.dat)

md.mdp

; md-qmmm-nvt.mdp

integrator = md ; MD using leap-frog integrator

dt = 0.001 ; 1fs time-step

nsteps = 1000 ; 1 ps simulation

; Set output frequency to each step

nstxout = 1 ; Coordinates to trr

nstvout = 1 ; velocities to trr

nstfout = 1 ; forces to trr

nstlog = 1 ; Energies to md.log

nstcalcenergy = 1 ; Energies to ener.edr

nstenergy = 1 ; Energies to ener.edr

nstxout-compressed = 1 ;Coordinates to xtc

compressed-x-precision = 1000

compressed-x-grps = system

energygrps =

; Set cut-offs

rlist = 1.2

coulombtype = PME

coulomb-modifier = Potential-shift-Verlet

rcoulomb-switch = 1.0

rcoulomb = 1.2

vdwtype = Cut-off

vdw-modifier = Force-switch

rvdw-switch = 1.0

rvdw = 1.2

cutoff-scheme = Verlet

nstlist = 1 ; Freq, to update neighbour list

pbc = xyz ; With periodic boundary conditions

;Temperature coupling options

tcoupl = v-rescale

nsttcouple = 10

tc-grps = System

tau-t = 0.1

ref-t = 300

qmmm-cp2k-active = true

qmmm-cp2k-qmgroup = QMatoms

qmmm-cp2k-qmmethod = INPUT ;PBE

;qmmm-cp2k-qmcharge = -3

;qmmm-cp2k-qmmultiplicity = 1

1) generate a new tpr

gmx_mpi grompp -f md-metad.mdp -qmi Qmmm-metad_cp2k.inp -p actin.top -c em-qmmm.gro -n index.ndx -o Qmmm-metad.tpr

2) write a plumed.dat file

#######################################################################################################

# define atoms for dist

d1: DISTANCE ATOMS=5829,6256

d2: DISTANCE ATOMS=5829,5832

m: MATHEVAL ARG=d1,d2 FUNC=x-y PERIODIC=NO

# metadynamics argument

metad: METAD ARG=m PACE=20 HEIGHT=1 SIGMA=0.05 FILE=HILLS BIASFACTOR=15 TEMP=300(for well tempered metadynamics)

# in HILLS the following variables will be written: time, angle, dist, sigma, height, biasfactor

# monitor angle

PRINT STRIDE=1 ARG=d1,d2,m,metad.bias FILE=COLVAR

# print the calculated angle and bias in COLVAR every 500 steps, COLVAR will contain three columns time, angl, dist and bias

#####################################################################################################

# define atoms for dist

d1: DISTANCE ATOMS=5829,6256

d2: DISTANCE ATOMS=5829,5832

m: MATHEVAL ARG=d1,d2 FUNC=x-y PERIODIC=NO

uwall: UPPER_WALLS ARG=m AT=0.18 KAPPA=500000.0 EXP=2 EPS=1 OFFSET=0

lwall: LOWER_WALLS ARG=m AT=-0.18 KAPPA=500000.0 EXP=2 EPS=1 OFFSET=0

# metadynamics argument

metad: METAD ARG=m PACE=20 HEIGHT=1 SIGMA=0.05 FILE=HILLS GRID_MIN=-0.2 GRID_MAX=0.2 GRID_BIN=500

# GRID_MIN should be smaller than lwall; GRID_MAX should be larger than uwall

# in HILLS the following variables will be written: time, angle, dist, sigma, height, biasfactor

PRINT STRIDE=1 ARG=d1,d2,m,uwall.bias,lwall.bias,metad.bias FILE=COLVAR

# print the calculated angle and bias in COLVAR every 500 steps, COLVAR will contain three columns time, angl, dist and bias

#####################################################################################################

#define atoms for dist

d1: DISTANCE ATOMS=5829,6256

d2: DISTANCE ATOMS=5829,5832

uwall: UPPER_WALLS ARG=d1,d2 AT=0.32,0.32 KAPPA=500000,500000 EXP=2,2 EPS=1,1 OFFSET=0,0

lwall: LOWER_WALLS ARG=d1,d2 AT=0.15,0.15 KAPPA=500000,500000 EXP=2,2 EPS=1,1 OFFSET=0,0

# metadynamics argument

metad: METAD ARG=d1,d2 PACE=20 HEIGHT=1 SIGMA=0.1,0.1 FILE=HILLS GRID_MIN=0.10,0.10 GRID_MAX=0.37,0.37 GRID_BIN=100,100

# in HILLS the following variables will be written: time, angle, dist, sigma, height, biasfactor

PRINT STRIDE=1 ARG=d1,d2,uwall.bias,lwall.bias,metad.bias FILE=COLVAR

# print the calculated angle and bias in COLVAR every 500 steps, COLVAR will contain three columns time, angl, dist and bias

#####################################################################################################

GRID_MIN: the minimum value your CV can adopt
GRID_MAX: the maximum value your CV can adopt
GRID_SPACING: the bin size，Remember that the bins in the grid should be small enough to describe the changes in the bias created by the Gaussian kernels. Normally a bin of size σ/5(with σ the gaussian width) is small enough.

https://www.plumed.org/doc-v2.8/user-doc/html/_d_i_s_t_a_n_c_e__f_r_o_m__c_o_n_t_o_u_r.html

In the simplest possible implementation of a metadynamics calculation the expense of a metadynamics calculation increases with the length of the simulation as one has to, at every step, evaluate the values of a larger and larger number of Gaussian kernels. To avoid this issue you can store the bias on a grid. This approach is similar to that proposed in [8] but has the advantage that the grid spacing is independent on the Gaussian width. Notice that you should provide the grid boundaries (GRID_MIN and GRID_MAX) and either the number of bins for every collective variable (GRID_BIN) or the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use the most conservative choice (highest number of bins) for each dimension. In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING) PLUMED will use 1/5 of the Gaussian width (SIGMA) as grid spacing if the width is fixed or 1/5 of the minimum Gaussian width (SIGMA_MIN) if the width is variable. This default choice should be reasonable for most applications.

##################################

Height 1 kj/Mol (even less),

Sigma 0.35 rad for dihedrals and 0.05 nm for distances,

the biasfactor you can choose 15, or 20, or higher,

Plumed and gromacs can do all these things, for defining "good" CVS, with physically relevant barriers, that's is more complicated, and very system dependent.

########################################

Plumed units

By default the PLUMED inputs and outputs quantities in the following units:

Energy - kJ/mol
Length - nanometers
Time - picoseconds

Restart the calculation

in plumed.dat file, adding

RESTART The RESTART action with label

For gromacs, run gmx_mpi mdrun -deffnm Qmmm-metad -cpi Qmmm-metad.cpt -plumed plumed.dat