PLUMED基础知识-01

发表于 | 分类于

基础知识

输入文件格式

第一个关键字指定要进行的 action，然后是对应的参数，可以是单个参数(eg. NOPBC 告诉plumed在计算cv时不使用周期性边界条件)或者参数后面接等号和内容(eg. ATOMS=1,2,3,4 告诉plumed在计算cv时使用1-4号原子)。

逗号分隔的列表可以使用空格加上大括号来替代，eg. ATOMS={1 2 3 4}

注释符号： #

连续行输入

1
2
DISTANCE ATOMS=1,2 LABEL=d1
PRINT ARG=d1 FILE=colvar STRIDE=10

等价的

1
2
d1: DISTANCE ATOMS=1,2
PRINT ARG=d1 FILE=colvar STRIDE=10

单位

默认( 输出? )单位:

1
2
3
Energy  - kJ/mol
Length  - nanometers
Time    - picoseconds

使用UNITS关键字来改变单位:

1
UNITS LENGTH=A

貌似从 .xyz 文件读取轨迹的话，长度会扩大十倍。

Groups and Virtual Atoms:

指定原子:

1
2
3
4
5
GROUP ATOMS=10,11,15,20 LABEL=g1                 逗号分隔，包含两端
GROUP ATOMS=10-20 LABEL=g2                       包含两端
GROUP ATOMS=10-100:10 LABEL=g3                   包含两端
GROUP ATOMS=100-10:-10 LABEL=g4                  包含两端
GROUP ATOMS=1,2,10-20,40-60:5,100-70:-2 LABEL=g5 混合

MOLINFO

使用 MOLINFO 和 TORSION 计算第一个和第四个残基的二面角 phi and psi:

1
2
3
4
5
MOLINFO MOLTYPE=protein STRUCTURE=myprotein.pdb
# t: TORSION ATOMS=1,2,3,4
t1: TORSION ATOMS=@phi-3
t2: TORSION ATOMS=@psi-4
PRINT ARG=t1,t2 FILE=colvar STRIDE=10

WHOLEMOLECULES

用来重建被周期性边界条件分割开的分子

1
2
3
4
5
计算包含100个原子的聚合物的 end-to-end distance

WHOLEMOLECULES ENTITY0=1-100
e2e: DISTANCE ATOMS=1,100 NOPBC
RESTRAINT ARG=e2e KAPPA=1 AT=5

Note:
DISTANCE的 NOPBC 参数，在距离大于模拟盒子的一半长度的时候使用。
由于某些MD软件中分子会穿过盒子的边界，所以可能还需要使用WHOLEMOLECULES 关键字。当然，需要写在 DISTANCE 之前。

1
2
3
4
5
6
7
在分子重建之前输出轨迹以查看差异 (100个原子)
结果发现复合物dump之后的坐标不对啊?
做的时候针对具体情况仔细检查

# DUMPATOMS FILE=dump-broken.xyz ATOMS=1-100 STRIDE=1000
WHOLEMOLECULES STRIDE=1 ENTITY0=1-100
DUMPATOMS FILE=dump.xyz ATOMS=1-100 STRIDE=1000

其它方式:

1
2
3
Using the FIT_TO_TEMPLATE they can be aligned to a template structure.
Using WRAPAROUND you can bring a set of atom as close as possible to another set of atoms.
Using RESET_CELL you can rotate the periodic cell.

Virtual Atoms

Keyword Description
CENTER_OF_MULTICOLVAR Calculate a a weighted average position based on the value of some multicolvar.
CENTER Calculate the center for a group of atoms, with arbitrary weights.
COM Calculate the center of mass for a group of atoms.
FIXEDATOM Add a virtual atom in a fixed position.
GHOST Calculate the absolute position of a ghost atom with fixed coordinatesin the local reference frame formed by three atoms. The computed ghost atom is stored as a virtual atom that can be accessed inan atom list through the the label for the GHOST action that creates it.

CENTER:

1
2
3
4
5
6
7
8
9
10
11
12
# 虚原子位于atom 1和10之间三分之一距离的位置

c1: CENTER ATOMS=1,1,10
# 等价于
c1bis: CENTER ATOMS=1,10 WEIGHTS=2,1

# 质心
c2: CENTER ATOMS=2,3,4,5 MASS

d1: DISTANCE ATOMS=c1,c2

PRINT ARG=d1

COM

1
2
3
4
COM ATOMS=1-7         LABEL=c1
COM ATOMS=15,20       LABEL=c2
DISTANCE ATOMS=c1,c2  LABEL=d1
PRINT ARG=d1

FIXEDATOM

在固定位置定义虚原子

1
2
3
4
5
6
计算两个原子与z轴之间的夹角

a: FIXEDATOM AT=0,0,0
b: FIXEDATOM AT=0,0,1
an: ANGLE ATOMS=a,b,15,20
RESTRAINT ARG=an AT=0.0 KAPPA=100.0
1
2
3
4
5
6
先align蛋白，然后计算其与固定点之间的距离

FIT_TO_TEMPLATE STRIDE=1 REFERENCE=ref.pdb TYPE=SIMPLE
a: FIXEDATOM AT=10,20,30
d: DISTANCE ARG=a,20
PRINT ARG=d FILE=colvar

CV Documentation

Keyword Description
ALPHABETA Measures a distance including pbc between the instantaneous values of a set of torsional angles and set of reference values.
ALPHARMSD Probe the alpha helical content of a protein structure.
ANGLE Calculate an angle.
ANTIBETARMSD Probe the antiparallel beta sheet content of your protein structure.
CELL Calculate the components of the simulation cell
CONSTANT Return one or more constant quantitieswith or without derivatives.
CONTACTMAP Calculate the distances between a number of pairs of atoms and transform each distance by a switching function.
COORDINATION Calculate coordination numbers.
CS2BACKBONE This collective variable calculates the backbone chemical shifts for a protein.
DHENERGY Calculate Debye-Huckel interaction energy among GROUPA and GROUPB.
DIHCOR Measures the degree of similarity between dihedral angles.
DIPOLE Calculate the dipole moment for a group of atoms.
DISTANCE_FROM_CONTOUR Calculate the perpendicular distance from a Willard-Chandler dividing surface.
DISTANCE Calculate the distance between a pair of atoms.
ENERGY Calculate the total energy of the simulation box.
ERMSD Calculate eRMSD with respect to a reference structure.
FAKE This is a fake colvar container used by cltools or various other actionsand just support input and period definition
FRET Calculate the FRET efficiency between a pair of atoms.The efficiency is calculated using the Forster relation:
GPROPERTYMAP Property maps but with a more flexible framework for the distance metric being used.
GYRATION Calculate the radius of gyration, or other properties related to it.
JCOUPLING Calculates 3J coupling constants for a dihedral angle.
NOE Calculates NOE intensities as sums of 1/r^6, also averaging over multiple equivalent atomsor ambiguous NOE.
PARABETARMSD Probe the parallel beta sheet content of your protein structure.
PATHMSD This Colvar calculates path collective variables.
PATH Path collective variables with a more flexible framework for the distance metric being used.
PCAVARS Projection on principal component eigenvectors or other high dimensional linear subspace
POSITION Calculate the components of the position of an atom.
PRE Calculates the Paramegnetic Resonance Enhancement intensity ratio between two atoms.
PROPERTYMAP Calculate generic property maps.
PUCKERING Calculate sugar pseudorotation coordinates.
RDC Calculates the (Residual) Dipolar Coupling between two atoms.
TEMPLATE This file provides a template for if you want to introduce a new CV.
TORSION Calculate a torsional angle.
VOLUME Calculate the volume of the simulation box.

Distances from reference configurations

Keyword Description
DRMSD Calculate the distance RMSD with respect to a reference structure.
MULTI-RMSD Calculate the RMSD distance moved by a number of separated domains from their positions in a reference structure.
PCARMSD Calculate the PCA components ( see [56] and [54] ) for a number of provided eigenvectors and an average structure. Performs optimal alignment at every step and reports the rmsd so you know if you are far or close from the average structure.It takes the average structure and eigenvectors in form of a pdb.Note that beta and occupancy values in the pdb are neglected and all the weights are placed to 1 (differently from the RMSD colvar for example)
RMSD Calculate the RMSD with respect to a reference structure.
TARGET This function measures the pythagorean distance from a particular structure measured in the space defined by some set of collective variables.

Functions

分析bias / CV的函数而不是CV本身

Keyword Description
COMBINE Calculate a polynomial combination of a set of other variables.
ENSEMBLE Calculates the replica averaging of a collective variable over multiple replicas.
FUNCPATHMSD This function calculates path collective variables.
FUNCSUMHILLS This function is intended to be called by the command line tool sum_hillsand it is meant to integrate a HILLS file or an HILLS file interpreted as a histogram i a variety of ways. Therefore it is not expected that you use this during your dynamics (it will crash!)
LOCALENSEMBLE Calculates the average over multiple arguments.
MATHEVAL Calculate a combination of variables using a matheval expression.
PIECEWISE Compute a piecewise straight line through its arguments that passes througha set of ordered control points.
SORT This function can be used to sort colvars according to their magnitudes.
STATS Calculates statistical properties of a set of collective variables with respect to a set of reference values.In particular it calculates and store as components the sum of the squared deviations, the correlation, theslope and the intercept of a linear fit.

MultiColvar

计算某一类CV的分布情况
eg. DISTANCES 能够计算一对或多对原子之间的距离，最小距离，小于特定值的距离的个数，处于特定范围的距离的个数..
离散到连续，还得查看具体的定义。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# 计算atom 3和5之间，以及1和2之间的距离，并输出两个距离的MIN(具体查看min的定义)。

DISTANCES ATOMS1=3,5 ATOMS2=1,2 MIN={BETA=0.1} LABEL=d1
PRINT ARG=d1.min


计算atom 3和5之间，以及1和2之间的距离，并输出小于2A (默认单位为nm，前面使用 UNIT LENGTH=A 关键字改变长度单位为A)的距离的个数。(具体查看LESS_THAN的定义)

DISTANCES ATOMS1=3,5 ATOMS2=1,2 LABEL=d1 LESS_THAN={RATIONAL R_0=2}
PRINT ARG=d1.lessthan


计算Group内所有原子之间的距离(1 and 2, 1 and 3, 2 and 3), 并输出平均值。

DISTANCES GROUP=1-3 MEAN LABEL=d1
PRINT ARG=d1.mean


计算GroupA与GroupB之间所有的距离(1 and 2, 1 and 3)，并输出大于1A的个数

DISTANCES GROUPA=1 GROUPB=2,3 MORE_THAN={RATIONAL R_0=1}
PRINT ARG=d1.morethan 


d1: DISTANCES GROUPA=1-10 GROUPB=11-20 MIN={BETA=500.} 
PRINT ARG=d1.min FILE=colvar STRIDE=10


DISTANCES ...
 GROUPA=1-10 GROUPB=11-20 
 LESS_THAN={RATIONAL R_0=1.0} 
 MORE_THAN={RATIONAL R_0=2.0} 
 BETWEEN={GAUSSIAN LOWER=1.0 UPPER=2.0} 
 MIN={BETA=500.}
... DISTANCES
PRINT ARG=d1.lessthan,d1.morethan,d1.between,d1.min FILE=colvar STRIDE=10

Exploiting contact matrices

Contact matrix 是N×N的矩阵，其第i行j列的元素表面第i/j个原子(分子)是否相邻。
有多种尺度定义是否相邻:

Keyword Description
ALIGNED_MATRIX Adjacency matrix in which two molecule are adjacent if they are within a certain cutoff and if they have the same orientation.
CONTACT_MATRIX Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff.
HBOND_MATRIX Adjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them.
SMAC_MATRIX Adjacency matrix in which two molecules are adjacent if they are within a certain cutoff and if the angle between them is within certain ranges.

计算 Contact matrix 之后，可进行一系列的分析:

Keyword Description
CLUSTER_WITHSURFACE Find the various connected components in an adjacency matrix and then output averageproperties of the atoms in those connected components.
COLUMNSUMS Sum the columns of a contact matrix
DFSCLUSTERING Find the connected components of the matrix using the DFS clustering algorithm.
ROWSUMS Sum the rows of a contact matrix
SPRINT Calculate SPRINT topological variables from an adjacency matrix.

继续分析:

Keyword Description
CLUSTER_DIAMETER Print out the diameter of one of the connected components.
CLUSTER_DISTRIBUTION Calculate functions of the distribution of properties in your connected components..
CLUSTER_NATOMS Gives the number of atoms in the connected component.
CLUSTER_PROPERTIES Calculate properties of the distribution of some quantities that are part of a connected component.
OUTPUT_CLUSTER Output the indices of the atoms in one of the clusters identified by a clustering object.

Analysis

PLUMED 既能在MD过程中进行分析，也能通过后处理进行分析。

Keyword Description
COMMITTOR Does a committor analysis.
DUMPATOMS Dump selected atoms on a file.
DUMPDERIVATIVES Dump the derivatives with respect to the input parameters for one or more objects (generally CVs, functions or biases).
DUMPFORCES Dump the force acting on one of a values in a file.
DUMPMASSCHARGE Dump masses and charges on a selected file.
DUMPMULTICOLVAR Dump atom positions and multicolvar on a file.
DUMPPROJECTIONS Dump the derivatives with respect to the input parameters for one or more objects (generally CVs, functions or biases).
PRINT Print quantities to a file.
UPDATE_IF Conditional update of other actions.

COMMITTOR

告诉PLUMED当某些条件满足之后停止计算，只能在MD时使用。

1
2
3
4
5
6
7
8
9
10
11
12
# 监控两个二面角，并定义两个basins，当模拟进入其中一个basin的时候停止

TORSION ATOMS=1,2,3,4 LABEL=r1
TORSION ATOMS=2,3,4,5 LABEL=r2
COMMITTOR ...
  ARG=r1,r2 
  STRIDE=10
  BASIN_LL1=0.15,0.20 
  BASIN_UL1=0.25,0.40 
  BASIN_LL2=-0.15,-0.20 
  BASIN_UL2=-0.25,-0.40 
... COMMITTOR

DUMPATOMS

在文件(.xyz or .gro)中输出指定原子的位置信息
Note:
如果编译的时候使用了xdrfile library，也能输出 xtc 和 trr 轨迹。

1
2
3
4
5
6
7
8
9
COM ATOMS=11-20 LABEL=c1
DUMPATOMS STRIDE=10 FILE=file.xyz ATOMS=1-10,c1

# this is required to have proper atom names:
MOLINFO STRUCTURE=reference.pdb
# if omitted, atoms will have "X" name...

COM ATOMS=11-20 LABEL=c1
DUMPATOMS STRIDE=10 FILE=file.gro ATOMS=1-10,c1

DUMPDERIVATIVES

输出导数

1
2
3
4
The following input instructs plumed to write a file called deriv that contains both the analytical and numerical derivatives of the distance between atoms 1 and 2.
DISTANCE ATOM=1,2 LABEL=distance
DISTANCE ATOM=1,2 LABEL=distanceN NUMERICAL_DERIVATIVES
DUMPDERIVATIVES ARG=distance,distanceN STRIDE=1 FILE=deriv

DUMPFORCES

Dump the force acting on one of a values in a file.

For a CV this command will dump the force on the CV itself. Be aware that in order to have the forces on the atoms you should multiply the output from this argument by the output from DUMPDERIVATIVES. Furthermore, also note that you can output the forces on multiple quantities simultaneously by specifying more than one argument.

1
2
3
4
Write a file called forces that contains the force acting on the distance between atoms 1 and 2.

DISTANCE ATOM=1,2 LABEL=distance
DUMPFORCES ARG=distance STRIDE=1 FILE=forces

DUMPMASSCHARGE

输出质量和电荷信息
只在MD中的第一步运行，输出文件能在driver中再次使用。如果不提供原子列表，则输出所有原子的信息。
一般放在输入文件的末尾。

1
2
3
4
5
c1: COM ATOMS=1-10
c2: COM ATOMS=11-20
PRINT ARG=c1,c2 FILE=colvar STRIDE=100

DUMPMASSCHARGE FILE=mcfile

这样就能通过-mc选项在driver命令中使用相同的质量来分析轨迹。

1
plumed driver --mc mcfile --plumed plumed.dat --ixyz traj.xyz

1
2
3
solute_ions: GROUP ATOMS=1-121,200-2012
DUMPATOMS FILE=traj.gro ATOMS=solute_ions STRIDE=100
DUMPMASSCHARGE FILE=mcfile ATOMS=solute_ions

Note: 如果希望通过driver来处理电荷(e.g. reading traj.gro)，则必须使用以下脚本来修改原子序号:

1
2
3
4
awk 'BEGIN{c=0}{
  if(match($0,"#")) print ; else {print c,$2,$3; c++}
}' < mc > newmc
}'

然后

1
plumed driver --mc newmc --plumed plumed.dat --ixyz traj.gro

DUMPMULTICOLVAR

输出原子位置和multicolvar

1
2
3
4
5
6
7
计算两个group之间的距离，并写入到文件。对于每一对距离，同时写入两个原子中点的坐标和距离。

pos:   GROUP ATOMS=220,221,235,236,247,248,438,439,450,451,534,535
neg:   GROUP ATOMS=65,68,138,182,185,267,270,291,313,316,489,583,621,711
DISTANCES GROUPA=pos GROUPB=neg LABEL=slt

DUMPMULTICOLVAR DATA=slt FILE=MULTICOLVAR.xyz

DUMPPROJECTIONS

Dump the derivatives with respect to the input parameters for one or more objects (generally CVs, functions or biases).

1
2
3
4
5
6
Compute the distance between two groups and write on a file the derivatives of this distance with respect to all the atoms of the two groups

x1: CENTER ATOMS=1-10
x2: CENTER ATOMS=11-20
d: DISTANCE ATOMS=x1,x2
DUMPPROJECTIONS ARG=d FILE=proj STRIDE=20

PRINT

1
2
3
4
DISTANCE ATOMS=2,5 LABEL=distance
ENERGY             LABEL=energy
PRINT ARG=distance          STRIDE=10   FILE=COLVAR
PRINT ARG=distance,energy   STRIDE=1000 FILE=COLVAR_ALL

FLUSH

This command instructs plumed to flush all the open files with a user specified frequency. Notice that all files are flushed anyway every 10000 steps.

1
2
3
4
5
6
7
8
9
10
# A command like this in the input will instruct plumed to flush all the output files every 100 steps

d1: DISTANCE ATOMS=1,10
PRINT ARG=d1 STRIDE=5 FILE=colvar1

FLUSH STRIDE=100

d2: DISTANCE ATOMS=2,11
# also this print is flushed every 100 steps:
PRINT ARG=d2 STRIDE=10 FILE=colvar2

UPDATE_IF

Conditional update of other actions.

1
2
3
4
5
6
7
8
根据特定的值来决定是否进行更新。比如，通过判断groupA 和groupB之间相互作用的数目来决定是否输出坐标。

solute: GROUP ATOMS=1-124
coord: COORDINATION GROUPA=solute GROUPB=500 R_0=0.5

UPDATE_IF ARG=coord LESS_THAN=0.5
DUMPATOMS ATOMS=solute,500 FILE=output.xyz
UPDATE_IF ARG=coord END

Reweight

以上的分析中，输出结果都是一种enseble average，如果MD是有偏的，则需要进行reweight来移除所施加的bias带来的影响。以下的方法可以用来计算权重:

Keyword Description
REWEIGHT_BIAS Calculate weights for ensemble averages that negate the effect the bias has on the region of phase space explored
REWEIGHT_METAD Calculate the weights configurations should contribute to the histogram in a simulation in which a metadynamics bias acts upon the system.
REWEIGHT_TEMP Calculate weights for ensemble averages allow for the computing of ensemble averages at temperatures lower/higher than that used in your original simulation.

然后通过以下方法计算系综平均:

Keyword Description
AVERAGE Calculate the ensemble average of a collective variable
HISTOGRAM Accumulate the average probability density along a few CVs from a trajectory.
MULTICOLVARDENS Evaluate the average value of a multicolvar on a grid.

For many of the above commands data is accumulated on the grids. These grids can be further analysed using one of the actions:

Keyword Description
CONVERT_TO_FES Convert a histogram, H(x), to a free energy surface using F(x)=−kBTlnH(x).
DUMPCUBE Output a three dimensional grid using the Gaussian cube file format.
DUMPGRID Output the function on the grid to a file with the PLUMED grid format.
FIND_CONTOUR_SURFACE Find an isocontour by searching along either the x, y or direction.
FIND_CONTOUR Find an isocontour in a smooth function.
FIND_SPHERICAL_CONTOUR Find an isocontour in a three dimensional grid by searching over a Fibonacci sphere.
FOURIER_TRANSFORM Compute the Discrete Fourier Transform (DFT) by means of FFTW of data stored on a 2D grid.
INTERPOLATE_GRID Interpolate a smooth function stored on a grid onto a grid with a smaller grid spacing.

example:

1
2
3
x: DISTANCE ATOMS=1,2
h: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1 STRIDE=5
DUMPGRID GRID=h FILE=histo STRIDE=100

Note: 进行上述分析的时候，第一帧总是被当做MD的初始构象而被忽略

1
2
3
x: DISTANCE ATOMS=1,2
h: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1 STRIDE=5
DUMPGRID GRID=h FILE=histo

计算平均值，PLUMED默认的方式是计算前 step N, step 2N… 的累积平均值；可以通过 CLEAR关键字来指定分别计算 first step N, second step N… 的平均值。

Dimensionality Reduction

降维

Keyword Description
CLASSICAL_MDS Create a low-dimensional projection of a trajectory using the classical multidimensional scaling algorithm.
PCA Perform principal component analysis (PCA) using either the positions of the atoms a large number of collective variables as input.

Bias

增强采样的方法:

Keyword Description
ABMD Adds a ratchet-and-pawl like restraint on one or more variables.
BIASVALUE Takes the value of one variable and use it as a bias
EXTENDED_LAGRANGIAN Add extended Lagrangian.
EXTERNAL Calculate a restraint that is defined on a grid that is read during start up
LOWER_WALLS Defines a wall for the value of one or more collective variables, which limits the region of the phase space accessible during the simulation.
METAD Used to performed MetaDynamics on one or more collective variables.
METAINFERENCE Calculate the Metainference Score for a set of back calculated experimental data.
MOVINGRESTRAINT Add a time-dependent, harmonic restraint on one or more variables.
PBMETAD Used to performed Parallel Bias MetaDynamics.
RESTRAINT Adds harmonic and/or linear restraints on one or more variables.
UPPER_WALLS Defines a wall for the value of one or more collective variables, which limits the region of the phase space accessible during the simulation.

METAD or PBMETAD，这类依赖历史信息的bias能通过RESTART关键字进行restart

ABMD

Adds a ratchet-and-pawl like restraint on one or more variables.
This action can be used to evolve a system towards a target value in CV space using an harmonic potential moving with the thermal fluctuations of the CV.

1
2
3
4
5
6
# The two target values are defined using TO and the two strength using KAPPA. The total energy of the bias is printed.

DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
ABMD ARG=d1,d2 TO=1.0,1.5 KAPPA=5.0,5.0 LABEL=abmd
PRINT ARG=abmd.bias,abmd.d1_min,abmd.d2_min

BIASVALUE

把CV的值作为bias
最简单的bias: bias势等于CV值。
eg: 把关于CV的函数作为bias势

1
2
3
4
5
6
The following input tells plumed to use the value of the distance between atoms 3 and 5 and the value of the distance between atoms 2 and 4 as biases. It then tells plumed to print the energy of the restraint.

DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=3,6 LABEL=d2
BIASVALUE ARG=d1,d2 LABEL=b
PRINT ARG=d1,d2,b.d1,b.d2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
t: TIME
# this just print cos and sin of time
cos: MATHEVAL ARG=t VAR=t FUNC=cos(t) PERIODIC=NO 
sin: MATHEVAL ARG=t VAR=t FUNC=sin(t) PERIODIC=NO
c1: COM ATOMS=1,2
c2: COM ATOMS=3,4
d: DISTANCE COMPONENTS ATOMS=c1,c2
PRINT ARG=t,cos,sin,d.x,d.y,d.z STRIDE=1 FILE=colvar FMT=%8.4f
# this calculates sine and cosine of a projected component of distance
mycos:  MATHEVAL ARG=d.x,d.y  VAR=x,y   FUNC=x/sqrt(x*x+y*y) PERIODIC=NO
mysin:  MATHEVAL ARG=d.x,d.y  VAR=x,y   FUNC=y/sqrt(x*x+y*y) PERIODIC=NO
# this creates a moving spring so that the system follows a circle-like dynamics 
# but it is not a bias, it is a simple value now
vv1:  MATHEVAL ARG=mycos,mysin,cos,sin VAR=mc,ms,c,s  FUNC=100*((mc-c)^2+(ms-s)^2) PERIODIC=NO
# this takes the value calculated with matheval and uses as a bias 
cc: BIASVALUE ARG=vv1 
# some printout
PRINT ARG=t,cos,sin,d.x,d.y,d.z,mycos,mysin,cc.bias.vv1 STRIDE=1 FILE=colvar FMT=%8.4f

LOWER_WALLS UPPER_WALLS

对一个或多个CV定义墙，来限定模拟能达到的相空间

1
2
3
4
5
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=uwall
LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=lwall
PRINT ARG=uwall.bias,lwall.bias

METAD

MetaDynamics

1
2
3
4
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

MOVINGRESTRAINT

Add a time-dependent, harmonic restraint on one or more variables.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
The following input is dragging the distance between atoms 2 and 4 from 1 to 2 in the first 1000 steps, then back in the next 1000 steps. In the following 500 steps the restraint is progressively switched off.

DISTANCE ATOMS=2,4 LABEL=d
MOVINGRESTRAINT ...
  ARG=d
  STEP0=0    AT0=1.0 KAPPA0=100.0
  STEP1=1000 AT1=2.0
  STEP2=2000 AT2=1.0
  STEP3=2500         KAPPA3=0.0
... MOVINGRESTRAINT


The following input is progressively building restraints distances between atoms 1 and 5 and between atoms 2 and 4 in the first 1000 steps. Afterwards, the restraint is kept static.

DISTANCE ATOMS=1,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
MOVINGRESTRAINT ...
  ARG=d1,d2 
  STEP0=0    AT0=1.0,1.5 KAPPA0=0.0,0.0
  STEP1=1000 AT1=1.0,1.5 KAPPA1=1.0,1.0
... MOVINGRESTRAINT

RESTRAINT

Adds harmonic and/or linear restraints on one or more variables.

Either or both of SLOPE and KAPPA must be present to specify the linear and harmonic force constants respectively.

1
2
3
4
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
RESTRAINT ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 LABEL=restraint
PRINT ARG=restraint.bias

RESTART

Activate restart.
此选项应当放在输入文件的最前面。
还可以仅对某个action进行restart，RESTART=AUTO表示针对整体，RESTART=YES or RESTART=NO表示对单个action的开或者关。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
d: DISTANCE ATOMS=1,2
PRINT ARG=d FILE=out
# the original file 'out' will be backed up. 


# On the other hand, using the following input:

RESTART
d: DISTANCE ATOMS=1,2
PRINT ARG=d FILE=out
# the file 'out' will be appended. 


# In the following case, file out1 will be backed up and file out2 will be concatenated

RESTART
d1: DISTANCE ATOMS=1,2
d2: DISTANCE ATOMS=1,2
PRINT ARG=d1 FILE=out1 RESTART=NO
PRINT ARG=d2 FILE=out2

Command Line Tools

PLUMED包含一系列的简单命令行工具，如下格式:

1
plumed <toolname> <list of input flags for that tool>

Keyword Description
driver-float Equivalent to driver, but using single precision reals.
driver driver is a tool that allows one to to use plumed to post-process an existing trajectory.
gentemplate gentemplate is a tool that you can use to construct template inputs for the variousactions
info This tool allows you to obtain information about your plumed version
kt Print out the value of kBT at a particular temperature
manual manual is a tool that you can use to construct the manual page for a particular action
simplemd simplemd allows one to do molecular dynamics on systems of Lennard-Jones atoms.
sum_hills sum_hills is a tool that allows one to to use plumed to post-process an existing hills/colvar file

driver

对已经存在的轨迹进行分析
driver并不事先知道质量和电荷，所以如果分析与质量和电荷有关的量，比如质心，就需要 -pdb or -mc 来读入质量或者电荷信息

1
plumed driver --plumed plumed.dat --ixyz trajectory.xyz

可以通过 --ixyz --igro --ixtc -itrr --mf_dcd 读入轨迹，其中--ixtc --itrr 需要提前安装xdrfile library，PLUMED编译的时候能自动检测到。

还可以通过--mf_dcd --mf_crd --mf_crdbox --mf_gro --mf_g96 --mf_trr --mf_xtc --mf_pdb读入轨迹。(需要VMD的支持)

Check the available molfile plugins and limitations.
To have support of all of VMD’s plugins you need to recompile PLUMED. You need to download the SOURCE of VMD, which contains a plugins directory. Adapt build.sh and compile it. At the end, you should get the molfile plugins compiled as a static library libmolfile_plugin.a. Locate said file and libmolfile_plugin.h, and customize the configure command with something along the lines of:

1
2
configure [...] LDFLAGS="-ltcl8.5 -L/mypathtomolfilelibrary/ -L/mypathtotcl" CPPFLAGS="-I/mypathtolibmolfile_plugin.h/"
and rebuild.

xdrfile的xtc或trr，比molfile的更稳健一些。

--trajectory-stride 设置为MD时保存轨迹的频率，-–timestep 设置为MD时的时间步长。

1
2
3
plumed driver --plumed plumed.dat --ixyz trajectory.xyz --trajectory-stride 1000 --timestep 0.002

plumed driver --plumed plumed.dat --pdb diala.pdb --mf_xtc traj.xtc --trajectory-stride 1000 --timestep 0.002

还可以在输入文件中使用READ关键字来读取MD过程中产生的CV文件.

1
2
3
rphi1:       READ FILE=input_colvar.data  VALUES=phi1
rphi2:       READ FILE=input_colvar.data  VALUES=phi2
PRINT ARG=rphi1,rphi2 STRIDE=500  FILE=output_colvar.data

gentemplate

为特定的 action 生成模板

1
plumed gentemplate --action DISTANCE

sum_hills

sum_hills 对已经存在的 hills或者cv文件进行后处理

1
plumed sum_hills  --hills PATHTOMYHILLSFILE

Tutorials

Keyword Description
Belfast tutorial Analyzing CVs This tutorial explains how to use plumed to analyze CVs
Belfast tutorial Adaptive variables I How to use path CVs
Belfast tutorial Adaptive variables II Dimensionality reduction and sketch maps
Belfast tutorial Umbrella sampling Umbrella sampling, reweighting, and weighted histogram
Belfast tutorial Out of equilibrium dynamics How to run a steered MD simulations and how to estimate the free energy
Belfast tutorial Metadynamics How to run a metadynamics simulation
Belfast tutorial Replica exchange I Parallel tempering and Metadynamics, Well-Tempered Ensemble
Belfast tutorial Replica exchange II and Multiple walkers Bias exchange and multiple walkers
Belfast tutorial NMR restraints NMR restraints
Belfast tutorial Steinhardt Parameters Steinhardt Parameters
Cambridge tutorial A short 2 hours tutorial that introduces Well-Tempered Metadynamics, Bias-Exchange Metadynamics and Replica-Average Metadynamics
Cineca tutorial A short 2 hours tutorial that introduces analysis, well-tempered metadynamics, and multiple-restraints umbrella sampling.
Using Hamiltonian replica exchange with GROMACS This tutorial explains how to use Hamiltonian replica exchange in GROMACS
Julich tutorial Developing CVs in plumed Implementing new collective variables in plumed
Moving from PLUMED 1 to PLUMED 2 This tutorial explains how plumed 1 input files can be translated into the new plumed 2 syntax.
Munster tutorial A short 3 hours tutorial that introduces analysis, well-tempered metadynamics, and multiple-restraints umbrella sampling.
http://www.youtube.com/watch?v=iDvZmbWE5ps A short video introduction to the use of multicolvars in PLUMED 2
http://www.youtube.com/watch?v=PxJP16qNCYs A short video introduction to the syntax of the PLUMED 2 input file
http://en.wikipedia.org/wiki/Metadynamics A wikipedia article on metadynamics

Belfast tutorial: Analyzing CVs

1
2
3
4
5
6
#my first plumed input:

DISTANCE ATOMS=2,253 LABEL=e2edist
PRINT ARG=e2edist STRIDE=1 FILE=COLVAR 
ENDPLUMED
here I can write what I want it won't be read.
1
plumed driver --plumed plumed.dat --ixyz trajectory-short.xyz --length-units 0.1
1
2
3
4
5
6
7
8
9
10
first: CENTER ATOMS=1,2,3,4,5,6
last: CENTER ATOMS=251-256

e2edist: DISTANCE ATOMS=2,253
comdist: DISTANCE ATOMS=first,last
cvtor: TORSION ATOMS=first,102,138,last

PRINT ARG=e2edist,comdist,cvtor STRIDE=1 FILE=COLVAR 

ENDPLUMED
1
2
3
4
5
6
7
8
# 分析整条链上所有CA碳之间的距离

ca: GROUP ATOMS=9,16,31,55,69,90,102,114,124,138,160,174,194,208,224,238
dd: DISTANCES GROUP=ca MEAN MIN={BETA=50} MAX={BETA=0.02} MOMENTS=2

PRINT ARG=dd.mean,dd.min,dd.max,dd.moment-2 STRIDE=1 FILE=COLVAR 

ENDPLUMED
1
2
3
4
5
6
MOLINFO STRUCTURE=template.pdb
abeta: ANTIBETARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} STRANDS_CUTOFF=1

PRINT ARG=abeta.lessthan STRIDE=1 FILE=COLVAR 

ENDPLUMED

Analysis of Collective Variables

计算自由能
CV通常被用来描述系统的自由能，给定一个固定原子数目的系统在固定温度、固定体积下，它将以如下概率采样到某构象(q为微观坐标)：

通过cv(也就是s(q))进行分析:

得到自由能:

也就是通过cv的概率分布即可得到系统在cv这个维度的自由能。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
MOLINFO STRUCTURE=template.pdb
abeta: ANTIBETARMSD RESIDUES=all TYPE=DRMSD LESS_THAN={RATIONAL R_0=0.08 NN=8 MM=12} STRANDS_CUTOFF=1
ca: GROUP ATOMS=9,16,31,55,69,90,102,114,124,138,160,174,194,208,224,238
DISTANCES ...
GROUP=ca MEAN MIN={BETA=50} MAX={BETA=0.02} MOMENTS=2 LABEL=dd
... DISTANCES 

PRINT ARG=abeta.lessthan,dd.mean,dd.min,dd.max,dd.moment-2 STRIDE=1 FILE=COLVAR 

HISTOGRAM ...
ARG=abeta.lessthan,dd.mean
LABEL=hh
KERNEL=DISCRETE
GRID_MIN=0,0.8
GRID_MAX=4,1.2
GRID_BIN=40,40
... HISTOGRAM

DUMPGRID GRID=hh FILE=histo

ENDPLUMED
1
2
3
4
5
6
7
8
9
10
11
12
HISTOGRAM ...
LABEL=hh
ARG=abeta.lessthan,dd.mean
GRID_MIN=0,0.8
GRID_MAX=4,1.2
GRID_SPACING=0.04,0.004
BANDWIDTH=0.08,0.008
... HISTOGRAM

DUMPGRID GRID=hh FILE=histo

ENDPLUMED