ISDB: setting up a Metadynamics Metainference simulation

Aims

The aim of this tutorial is to introduce the users to the ISDB module and in particular to Metadynamics Metainference [23] [24] ensemble determination. We will reproduce the setup of the simulation for a simple system [86] . For a general overview of the problem of ensembles determination please read [25] .

Objectives

Once this tutorial is completed students will be able to:

  • Setup their own PLUMED-ISDB simulation.

Resources

The TARBALL for this project contains the following files:

  • charmm36-eef1sb.ff: the force-field files for gromacs (not needed)
  • system: a folder with reference files for gromacs (not needed)
  • reference-impl: a folder to perform a simple implicit solvent simulation
  • reference-impl-pbmetad: a folder to perform a pbmetad implicit solvent simulation
  • m_and_m: a folder to perform a metadynamics metainference simulation

Introduction

Molecular dynamics simulations are the ideal tool to determine at atomistic resolution the behavior of complex molecules. This great resolution power comes at the cost of approximations that affects the agreement with actual experimental observables. At the same time experimental data alone are generally speaking not enough to determine a structural ensemble due the inverse nature of the problem, that is to go from few observables to many atoms in many different configurations. Furthermore, experimental data are affected by errors of multiple nature, from noise, systematic errors and errors in their atomistic interpretation. Most important experimental data are the result of the averaging over the ensemble of structure so it is not trivial to deconvolve this signal. One possibility is that of employing MD simulations together with experimental data to generate simulations already corrected for the data themselves. With METAINFERENCE this is done on-the-fly by adding an additional energy to the system that takes into account the agreement with the experimental data considering the multiple sources of errors.

Run a reference simulation

The system we use is the EGAAWAASS peptide used in ref. [86] . First of all we will run a simulation in implicit solvent using the EEF1-SB CHARMM36 force field. EEF1-SB includes a correction to the standard backbone torsion potential of CHARMM36, an electrostatic interaction with a distance dependent dielectric constant and a simple gaussian form for the solvation energy. The first two terms are implemented in the force field and using table potentials while the latter is implemented as a collective variable in PLUMED, EEFSOLV .

Click on the labels of the actions for more information on what each action computes
tested on v2.9
# this is optional and tell to VIM that this is a PLUMED file
# vim: ft=plumed
# see comments just below this input file
#SETTINGS MOLFILE=user-doc/tutorials/others/isdb-1/reference-impl/egaawaass.pdb
MOLINFO MOLTYPEcompulsory keyword ( default=protein )
what kind of molecule is contained in the pdb file - usually not needed since protein/RNA/DNA
are compatible =protein STRUCTUREcompulsory keyword 
a file in pdb format containing a reference structure. =egaawaass.pdb 
WHOLEMOLECULES ENTITY0the atoms that make up a molecule that you wish to align. =1-111 
# EEF1SB Implicit solvation
#SETTINGS AUXFILE=user-doc/tutorials/others/isdb-1/reference-impl/index.ndx
protein-h: GROUP NDX_FILEthe name of index file (gromacs syntax) =index.ndx NDX_GROUPthe name of the group to be imported (gromacs syntax) - first group found is used
by default =Protein-H 
solv: EEFSOLV ATOMSThe atoms to be included in the calculation, e.g. =protein-h NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 NL_STRIDEcompulsory keyword ( default=40 )
The frequency with which the neighbor list is updated. =20 NL_BUFFERcompulsory keyword ( default=0.1 )
The buffer to the intrinsic cutoff used when calculating pairwise interactions. =0.1
bias: BIASVALUE ARG=solv

This can be run using gromacs (unfortunately recent versions of gromacs do not support Verlet groups with table potentials, so performances are currently sub-optimal on the gromacs side)

gmx_mpi mdrun -s run.tpr -table table.xvg -tablep table.xvg -plumed plumed-eef1.dat -v

In order to have a converged sampling for this reference ensemble calculation it is useful to setup a Metadynamics calculation. In particular we will use PBMETAD because it is then a natural choice for Metadynamics Metainference later.

Click on the labels of the actions for more information on what each action computes
tested on v2.9
# this is optional and tell to VIM that this is a PLUMED file
# vim: ft=plumed
# see comments just below this input file
#SETTINGS MOLFILE=user-doc/tutorials/others/isdb-1/reference-impl/egaawaass.pdb
MOLINFO MOLTYPEcompulsory keyword ( default=protein )
what kind of molecule is contained in the pdb file - usually not needed since protein/RNA/DNA
are compatible =protein STRUCTUREcompulsory keyword 
a file in pdb format containing a reference structure. =egaawaass.pdb 
WHOLEMOLECULES ENTITY0the atoms that make up a molecule that you wish to align. =1-111 
# EEF1SB Implicit solvation
#SETTINGS AUXFILE=user-doc/tutorials/others/isdb-1/reference-impl/index.ndx
protein-h: GROUP NDX_FILEthe name of index file (gromacs syntax) =index.ndx NDX_GROUPthe name of the group to be imported (gromacs syntax) - first group found is used
by default =Protein-H 
solv: EEFSOLV ATOMSThe atoms to be included in the calculation, e.g. =protein-h NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 NL_STRIDEcompulsory keyword ( default=40 )
The frequency with which the neighbor list is updated. =20 NL_BUFFERcompulsory keyword ( default=0.1 )
The buffer to the intrinsic cutoff used when calculating pairwise interactions. =0.1 
bias: BIASVALUE ARGthe input for this action is the scalar output from one or more other actions. =solv 
# CVs, Psi9, Phi1 are not defined
psi1: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-1 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi2: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-2 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi3: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-3 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi4: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-4 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi5: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-5 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi6: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-6 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi7: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-7 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi8: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-8 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi2: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi3: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-3 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi4: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-4 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi5: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-5 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi6: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-6 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi7: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-7 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi8: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-8 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi9: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-9 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
ahc: ALPHARMSD RESIDUESthis command is used to specify the set of residues that could conceivably form part
of the secondary structure. =all TYPEcompulsory keyword ( default=DRMSD )
the manner in which RMSD alignment is performed. =OPTIMAL LESS_THANcalculate the number of variables less than a certain target value. ={RATIONAL R_0=0.12}  
# Bulky Trp residue dihedral
dihtrp_cacb: TORSION ATOMSthe four atoms involved in the torsional angle =67,47,49,52 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
dihtrp_cbcg: TORSION ATOMSthe four atoms involved in the torsional angle =47,49,52,53 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
protein-ca: GROUP NDX_FILEthe name of index file (gromacs syntax) =index.ndx NDX_GROUPthe name of the group to be imported (gromacs syntax) - first group found is used
by default =C-alpha 
gyr: GYRATION TYPEcompulsory keyword ( default=RADIUS )
The type of calculation relative to the Gyration Tensor you want to perform =RADIUS ATOMSthe group of atoms that you are calculating the Gyration Tensor for. =protein-ca NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
# PBMetaD
pb: PBMETAD ...
   
   ARGthe input for this action is the scalar output from one or more other actions. =phi2,phi3,phi4,phi5,phi6,phi7,phi8,phi9,psi1,psi2,psi3,psi4,psi5,psi6,psi7,psi8,dihtrp_cacb,dihtrp_cbcg,ahc.lessthan 
   SIGMAcompulsory keyword 
the widths of the Gaussian hills =0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.01 
   HEIGHTthe height of the Gaussian hills, one for all biases. =0.5 
   PACEcompulsory keyword 
the frequency for hill addition, one for all biases =400 
   BIASFACTORuse well tempered metadynamics with this bias factor, one for all biases. =20 
   GRID_MINthe lower bounds for the grid =-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,0 
   GRID_MAXthe upper bounds for the grid =pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,5 
   GRID_WSTRIDEfrequency for dumping the grid =5000 
   WALKERS_MPI( default=off ) Switch on MPI version of multiple walkers - not compatible with WALKERS_*
options other than WALKERS_DIR  
...
PRINT FILEthe name of the file on which to output these quantities =COLVAR ARGthe input for this action is the scalar output from one or more other actions. =phi2,phi3,phi4,phi5,phi6,phi7,phi8,phi9,psi1,psi2,psi3,psi4,psi5,psi6,psi7,psi8,dihtrp_cacb,dihtrp_cbcg,ahc.lessthan STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200 
PRINT FILEthe name of the file on which to output these quantities =ENERGY ARGthe input for this action is the scalar output from one or more other actions. =bias.bias,pb.bias STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200

In this case we are running a multiple-replica simulation where the sampling is used to parallelize the Metadynamics time-dependent potential through the use of multiple walkers.

mpiexec -np 14 gmx_mpi mdrun -s topolnew -multi 14 -plumed plumed-eef1-pbmetad.dat -table table.xvg -tablep table.xvg >& log.out &

Metadynamics Metainference

The former simulations should provide a converged (check for this) ensemble for the peptide. As shown in [86] the agreement with the multiple available NMR experimental data is not perfect. In order to generate an ensemble compatible with most of the available experimental data it is possible to include them in the simulation using METAINFERENCE . To do so the forward models for the data sets should be defined in the input file. In this case we have backbone chemical shifts, CS2BACKBONE ; residual dipolar couplings for two bonds, RDC ; and J-couplings for multiple atoms, JCOUPLING. Once the forward models are defined for the data sets, the calculated data together with the corresponding experimental values can be used to calculate the metainference score. The metainference score is additive so it can be split into multiple METAINFERENCE entries. In this case we are using two metainference entries for the two sets of RDCs because these are compared with the experimental data modulo a constant that should be unique each data set. Then we use one metainference for all the jcouplings and another one for the chemical shifts.

Click on the labels of the actions for more information on what each action computes
tested on v2.9
# this is optional and tell to VIM that this is a PLUMED file
# vim: ft=plumed
# see comments just below this input file
#SETTINGS NREPLICAS=2
#SETTINGS MOLFILE=user-doc/tutorials/others/isdb-1/reference-impl/egaawaass.pdb
MOLINFO MOLTYPEcompulsory keyword ( default=protein )
what kind of molecule is contained in the pdb file - usually not needed since protein/RNA/DNA
are compatible =protein STRUCTUREcompulsory keyword 
a file in pdb format containing a reference structure. =egaawaass.pdb 
WHOLEMOLECULES ENTITY0the atoms that make up a molecule that you wish to align. =1-111 
# EEF1SB Implicit solvation
#SETTINGS AUXFILE=user-doc/tutorials/others/isdb-1/reference-impl/index.ndx
protein-h: GROUP NDX_FILEthe name of index file (gromacs syntax) =index.ndx NDX_GROUPthe name of the group to be imported (gromacs syntax) - first group found is used
by default =Protein-H 
solv: EEFSOLV ATOMSThe atoms to be included in the calculation, e.g. =protein-h NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 NL_STRIDEcompulsory keyword ( default=40 )
The frequency with which the neighbor list is updated. =20 NL_BUFFERcompulsory keyword ( default=0.1 )
The buffer to the intrinsic cutoff used when calculating pairwise interactions. =0.1 
bias: BIASVALUE ARGthe input for this action is the scalar output from one or more other actions. =solv 
# CVs, Psi9, Phi1 are not defined
psi1: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-1 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi2: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-2 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi3: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-3 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi4: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-4 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi5: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-5 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi6: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-6 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi7: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-7 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
psi8: TORSION ATOMSthe four atoms involved in the torsional angle =@psi-8 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi2: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-2 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi3: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-3 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi4: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-4 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi5: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-5 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi6: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-6 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi7: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-7 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi8: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-8 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
phi9: TORSION ATOMSthe four atoms involved in the torsional angle =@phi-9 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
ahc: ALPHARMSD RESIDUESthis command is used to specify the set of residues that could conceivably form part
of the secondary structure. =all TYPEcompulsory keyword ( default=DRMSD )
the manner in which RMSD alignment is performed. =OPTIMAL LESS_THANcalculate the number of variables less than a certain target value. ={RATIONAL R_0=0.12}  
# Bulky Trp residue dihedral
dihtrp_cacb: TORSION ATOMSthe four atoms involved in the torsional angle =67,47,49,52 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
dihtrp_cbcg: TORSION ATOMSthe four atoms involved in the torsional angle =47,49,52,53 NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
protein-ca: GROUP NDX_FILEthe name of index file (gromacs syntax) =index.ndx NDX_GROUPthe name of the group to be imported (gromacs syntax) - first group found is used
by default =C-alpha 
gyr: GYRATION TYPEcompulsory keyword ( default=RADIUS )
The type of calculation relative to the Gyration Tensor you want to perform =RADIUS ATOMSthe group of atoms that you are calculating the Gyration Tensor for. =protein-ca NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
# PBMetaD
pb: PBMETAD ...
   
   ARGthe input for this action is the scalar output from one or more other actions. =phi2,phi3,phi4,phi5,phi6,phi7,phi8,phi9,psi1,psi2,psi3,psi4,psi5,psi6,psi7,psi8,dihtrp_cacb,dihtrp_cbcg,ahc.lessthan 
   SIGMAcompulsory keyword 
the widths of the Gaussian hills =0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.01 
   HEIGHTthe height of the Gaussian hills, one for all biases. =0.5 
   PACEcompulsory keyword 
the frequency for hill addition, one for all biases =400 
   BIASFACTORuse well tempered metadynamics with this bias factor, one for all biases. =20 
   GRID_MINthe lower bounds for the grid =-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,-pi,0 
   GRID_MAXthe upper bounds for the grid =pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,pi,5 
   GRID_WSTRIDEfrequency for dumping the grid =5000 
   WALKERS_MPI( default=off ) Switch on MPI version of multiple walkers - not compatible with WALKERS_*
options other than WALKERS_DIR  
...
PRINT FILEthe name of the file on which to output these quantities =COLVAR ARGthe input for this action is the scalar output from one or more other actions. =phi2,phi3,phi4,phi5,phi6,phi7,phi8,phi9,psi1,psi2,psi3,psi4,psi5,psi6,psi7,psi8,dihtrp_cacb,dihtrp_cbcg,ahc.lessthan STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200 
PRINT FILEthe name of the file on which to output these quantities =ENERGY ARGthe input for this action is the scalar output from one or more other actions. =bias.bias,pb.bias STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200 
# EXPERIMENTAL DATA SECTION
# RDCs (Grzesiek et al.)
# xGAAWAASS
nh: RDC ...
   GYROMcompulsory keyword ( default=1. )
Add the product of the gyromagnetic constants for the bond. =-72.5388 
   SCALEcompulsory keyword ( default=1. )
Add the scaling factor to take into account concentration and other effects. =0.0001 
   NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
   ATOMS1the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =18,19 COUPLING1Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-5.4 
   ATOMS2the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =25,26 COUPLING2Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-1.26 
   ATOMS3the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =35,36 COUPLING3Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-5.22 
   ATOMS4the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =45,46 COUPLING4Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-0.91 
   ATOMS5the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =69,70 COUPLING5Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=2.33 
   ATOMS6the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =79,80 COUPLING6Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-2.88 
   ATOMS7the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =89,90 COUPLING7Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-8.37 
   ATOMS8the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =100,101 COUPLING8Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-3.78 
   
...
# ExAAWAASx
caha: RDC ...
   GYROMcompulsory keyword ( default=1. )
Add the product of the gyromagnetic constants for the bond. =179.9319 
   SCALEcompulsory keyword ( default=1. )
Add the scaling factor to take into account concentration and other effects. =0.0001 
   NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
   ATOMS1the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =5,6 COUPLING1Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=12.95 
   ATOMS2the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =27,28 COUPLING2Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=11.5 
   ATOMS3the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =37,38 COUPLING3Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=21.42 
   ATOMS4the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =47,48 COUPLING4Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=-9.37 
   ATOMS5the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =71,72 COUPLING5Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=10.01 
   ATOMS6the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =81,82 COUPLING6Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=15.01 
   ATOMS7the couple of atoms involved in each of the bonds for which you wish to calculate
the RDC. =91,92 COUPLING7Add an experimental value for each coupling (needed by SVD and useful for \ref STATS)..
=15.73 
   
...
#RDCS
byrdcnh: METAINFERENCE ...
   ARGthe input for this action is the scalar output from one or more other actions. =(nh\.rdc-.*),pb.bias 
   PARARGreference values for the experimental data, these can be provided as arguments without
derivatives =(nh\.exp-.*) 
   REWEIGHT( default=off ) simple REWEIGHT using the latest ARG as energy  
   NOISETYPEcompulsory keyword ( default=MGAUSS )
functional form of the noise (GAUSS,MGAUSS,OUTLIERS,MOUTLIERS,GENERIC) =MGAUSS 
   OPTSIGMAMEANcompulsory keyword ( default=NONE )
Set to NONE/SEM to manually set sigma mean, or to estimate it on the fly =SEM_MAX AVERAGINGStride for calculation of averaged weights and sigma_mean =400 
   SCALEDATA( default=off ) Set to TRUE if you want to sample a scaling factor common to all
values and replicas  SCALE_PRIORcompulsory keyword ( default=FLAT )
either FLAT or GAUSSIAN =GAUSSIAN SCALE0 could not find this keyword =8.0 DSCALEmaximum MC move of the scaling factor =0.5 
   SIGMA0 could not find this keyword =5.0 SIGMA_MINcompulsory keyword ( default=0.0 )
minimum value of the uncertainty parameter =0.0001 SIGMA_MAXcompulsory keyword ( default=10. )
maximum value of the uncertainty parameter =15.0 DSIGMAmaximum MC move of the uncertainty parameter =0.1 
   WRITE_STRIDEcompulsory keyword ( default=10000 )
write the status to a file every N steps, this can be used for restart/continuation
=10000 
   
...
#RDCS
byrdccaha: METAINFERENCE ...
   ARGthe input for this action is the scalar output from one or more other actions. =(caha\.rdc-.*),pb.bias 
   PARARGreference values for the experimental data, these can be provided as arguments without
derivatives =(caha\.exp-.*) 
   REWEIGHT( default=off ) simple REWEIGHT using the latest ARG as energy  
   NOISETYPEcompulsory keyword ( default=MGAUSS )
functional form of the noise (GAUSS,MGAUSS,OUTLIERS,MOUTLIERS,GENERIC) =MGAUSS 
   OPTSIGMAMEANcompulsory keyword ( default=NONE )
Set to NONE/SEM to manually set sigma mean, or to estimate it on the fly =SEM_MAX AVERAGINGStride for calculation of averaged weights and sigma_mean =400 
   SCALEDATA( default=off ) Set to TRUE if you want to sample a scaling factor common to all
values and replicas  SCALE_PRIORcompulsory keyword ( default=FLAT )
either FLAT or GAUSSIAN =GAUSSIAN SCALE0 could not find this keyword =9.0 DSCALEmaximum MC move of the scaling factor =0.5 
   SIGMA0 could not find this keyword =5.0 SIGMA_MINcompulsory keyword ( default=0.0 )
minimum value of the uncertainty parameter =0.0001 SIGMA_MAXcompulsory keyword ( default=10. )
maximum value of the uncertainty parameter =15.0 DSIGMAmaximum MC move of the uncertainty parameter =0.1 
   WRITE_STRIDEcompulsory keyword ( default=10000 )
write the status to a file every N steps, this can be used for restart/continuation
=10000 
   
...
# xGxAWxASx
jhan: JCOUPLING ...
   TYPEcompulsory keyword 
Type of J-coupling to compute (HAN,HAHN,CCG,NCG,CUSTOM) =HAN 
   NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
   ATOMS1the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@psi-2 COUPLING1Add an experimental value for each coupling. =-0.49 
   ATOMS2the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@psi-4 COUPLING2Add an experimental value for each coupling. =-0.54 
   ATOMS3the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@psi-5 COUPLING3Add an experimental value for each coupling. =-0.53 
   ATOMS4the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@psi-7 COUPLING4Add an experimental value for each coupling. =-0.39 
   ATOMS5the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@psi-8 COUPLING5Add an experimental value for each coupling. =-0.39 
   
...
# xxAAWAASS
jhahn: JCOUPLING ...
   TYPEcompulsory keyword 
Type of J-coupling to compute (HAN,HAHN,CCG,NCG,CUSTOM) =HAHN 
   NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
   ATOMS1the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@phi-2 COUPLING1Add an experimental value for each coupling. =6.05 
   ATOMS2the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@phi-3 COUPLING2Add an experimental value for each coupling. =5.95 
   ATOMS3the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@phi-4 COUPLING3Add an experimental value for each coupling. =6.44 
   ATOMS4the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@phi-5 COUPLING4Add an experimental value for each coupling. =6.53 
   ATOMS5the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@phi-6 COUPLING5Add an experimental value for each coupling. =5.93 
   ATOMS6the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@phi-7 COUPLING6Add an experimental value for each coupling. =6.98 
   ATOMS7the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@phi-8 COUPLING7Add an experimental value for each coupling. =7.16 
   
...
# xxxxWxxxx
jccg: JCOUPLING ...
   TYPEcompulsory keyword 
Type of J-coupling to compute (HAN,HAHN,CCG,NCG,CUSTOM) =CCG 
   NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
   ATOMS1the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@chi1-5 COUPLING1Add an experimental value for each coupling. =1.59 
   
...
# xxxxWxxxx
jncg: JCOUPLING ...
   TYPEcompulsory keyword 
Type of J-coupling to compute (HAN,HAHN,CCG,NCG,CUSTOM) =NCG 
   NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 
   ATOMS1the 4 atoms involved in each of the bonds for which you wish to calculate the J-coupling.
=@chi1-5 COUPLING1Add an experimental value for each coupling. =1.21 
   
...
#JC METAINFERENCE
byj: METAINFERENCE ...
   ARGthe input for this action is the scalar output from one or more other actions. =(jhan\.j-.*),(jhahn\.j-.*),(jccg\.j.*),(jncg\.j.*),pb.bias 
   PARARGreference values for the experimental data, these can be provided as arguments without
derivatives =(jhan\.exp-.*),(jhahn\.exp-.*),(jccg\.exp.*),(jncg\.exp.*) 
   REWEIGHT( default=off ) simple REWEIGHT using the latest ARG as energy  
   NOISETYPEcompulsory keyword ( default=MGAUSS )
functional form of the noise (GAUSS,MGAUSS,OUTLIERS,MOUTLIERS,GENERIC) =MGAUSS 
   OPTSIGMAMEANcompulsory keyword ( default=NONE )
Set to NONE/SEM to manually set sigma mean, or to estimate it on the fly =SEM_MAX AVERAGINGStride for calculation of averaged weights and sigma_mean =400 
   SIGMA0 could not find this keyword =5.0 SIGMA_MINcompulsory keyword ( default=0.0 )
minimum value of the uncertainty parameter =0.0001 SIGMA_MAXcompulsory keyword ( default=10. )
maximum value of the uncertainty parameter =15.0 DSIGMAmaximum MC move of the uncertainty parameter =0.1 
   WRITE_STRIDEcompulsory keyword ( default=10000 )
write the status to a file every N steps, this can be used for restart/continuation
=10000 
   
...

# Chemical shifts
#SETTINGS AUXFOLDER=user-doc/tutorials/others/isdb-1/m_and_m/data
# Chemical shifts
cs: CS2BACKBONE ...
   NOPBC( default=off ) ignore the periodic boundary conditions when calculating distances
 ATOMSThe atoms to be included in the calculation, e.g. =1-111 DATADIRcompulsory keyword ( default=data/ )
The folder with the experimental chemical shifts. =data TEMPLATEcompulsory keyword ( default=template.pdb )
A PDB file of the protein system. =egaawaass.pdb DOSCORE( default=off ) activate metainference  ARGthe input for this action is the scalar output from one or more other actions. =pb.bias 
   NOISETYPEcompulsory keyword ( default=MGAUSS )
functional form of the noise (GAUSS,MGAUSS,OUTLIERS,MOUTLIERS,GENERIC) =MGAUSS REWEIGHT( default=off ) simple REWEIGHT using the ARG as energy  OPTSIGMAMEANcompulsory keyword ( default=NONE )
Set to NONE/SEM to manually set sigma mean, or to estimate it on the fly =SEM_MAX AVERAGINGStride for calculation of averaged weights and sigma_mean =400 
   SIGMA0 could not find this keyword =4.0 
   SIGMA_MINcompulsory keyword ( default=0.0 )
minimum value of the uncertainty parameter =0.0001 
   SIGMA_MAXcompulsory keyword ( default=10. )
maximum value of the uncertainty parameter =5.00 
   WRITE_STRIDEcompulsory keyword ( default=10000 )
write the status to a file every N steps, this can be used for restart/continuation
=10000 
...
mcs: BIASVALUE ARGthe input for this action is the scalar output from one or more other actions. =cs.score 
# output from METAINFERENCE
PRINT ARGthe input for this action is the scalar output from one or more other actions. =byrdcnh.* STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200 FILEthe name of the file on which to output these quantities =BAYES.RDC.NH 
PRINT ARGthe input for this action is the scalar output from one or more other actions. =byrdccaha.* STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200 FILEthe name of the file on which to output these quantities =BAYES.RDC.CAHA 
PRINT ARGthe input for this action is the scalar output from one or more other actions. =byj.* STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200 FILEthe name of the file on which to output these quantities =BAYES.J 
PRINT ARGthe input for this action is the scalar output from one or more other actions. =cs.* STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =200 FILEthe name of the file on which to output these quantities =BAYES.CS 
# the following are useful for the analysis on-the-fly of the quality of the agreement with the experimentl data
ens: ENSEMBLE ...
   ARGthe input for this action is the scalar output from one or more other actions. =(nh\.rdc-.*),(caha\.rdc-.*),(jhan\.j-.*),(jhahn\.j-.*),(jccg\.j-.*),(jncg\.j-.*),(cs\...-.*),pb.bias REWEIGHT( default=off ) simple REWEIGHT using the latest ARG as energy  
   
...
nhst: STATS ...
   ARGthe input for this action is the scalar output from one or more other actions. =(ens\.nh\.rdc-.*) PARARGthe input for this action is the scalar output from one or more other actions without
derivatives. =(nh\.exp-.*) 
   
...
cahast: STATS ...
   ARGthe input for this action is the scalar output from one or more other actions. =(ens\.caha\.rdc-.*) PARARGthe input for this action is the scalar output from one or more other actions without
derivatives. =(caha\.exp-.*) 
   
...
csst: STATS ...
   ARGthe input for this action is the scalar output from one or more other actions. =(ens\.cs\...-.*) PARARGthe input for this action is the scalar output from one or more other actions without
derivatives. =(cs\.exp.*) 
   
...
jhanst: STATS ...
   ARGthe input for this action is the scalar output from one or more other actions. =(ens\.jhan\.j-.*) PARARGthe input for this action is the scalar output from one or more other actions without
derivatives. =(jhan\.exp-.*) 
   
...
jhahnst: STATS ...
   ARGthe input for this action is the scalar output from one or more other actions. =(ens\.jhahn\.j-.*) PARARGthe input for this action is the scalar output from one or more other actions without
derivatives. =(jhahn\.exp-.*) 
   
...
jw5ccyst: STATS ...
   ARGthe input for this action is the scalar output from one or more other actions. =(ens\.jccg\.j.*),(ens\.jccg\.j.*) PARARGthe input for this action is the scalar output from one or more other actions without
derivatives. =(jccg\.exp-.*),(jccg\.exp-.*) 
   SQDEVSUM( default=off ) calculates only SQDEVSUM  
   
...
jw5ncyst: STATS ...
   ARGthe input for this action is the scalar output from one or more other actions. =(ens\.jncg\.j.*),(ens\.jncg\.j.*) PARARGthe input for this action is the scalar output from one or more other actions without
derivatives. =(jncg\.exp-.*),(jncg\.exp-.*) 
   SQDEVSUM( default=off ) calculates only SQDEVSUM  
   
...
#output from STATS
PRINT ARGthe input for this action is the scalar output from one or more other actions. =nhst.* STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =2000 FILEthe name of the file on which to output these quantities =ST.RDC.NH 
PRINT ARGthe input for this action is the scalar output from one or more other actions. =cahast.* STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =2000 FILEthe name of the file on which to output these quantities =ST.RDC.CAHA 
PRINT ARGthe input for this action is the scalar output from one or more other actions. =csst.* STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output =2000 FILEthe name of the file on which to output these quantities =ST.CS
PRINT ARG=jhanst.*,jhahnst.*,jw5ccyst.*,jw5ncyst.* STRIDE=2000 FILE=ST.J

As for the former case we are running a multiple-replica simulation where in addition to multiple-walker metadynamics we are also coupling the replicas through Metainference. The use of multiple-walkers metadynamics is here key in order to have the same bias defined for all the replicas. This allows us to calculate a weighted average of the experimental observables where the weights are defined univocally from the bias [24] .

mpiexec -np 14 gmx_mpi mdrun -s topolnew -multi 14 -plumed plumed-eef1-pbmetad-m_m.dat -table table.xvg -tablep table.xvg >& log.out &