CP2K中metadynamics输入文件编辑(DISTANCE)

Keith Lv2

随着最近alphafold-latest和RFAA等模型的出现,给人一种just wait的感觉,所以日后也会尝试去更新一些AI指导的分子生成文章,由于最近在学习CP2K,同时又觉得使用CP2K去进行metadynamics的过程或者教程不是特别详细或者是比较乱?因此本守旧派想要在这里对metadynamics在CP2K中的一些options作一下简要解释并同时给出一些较为常见的集合变量定义的模板,目前暂时先介绍一下关于DISTANCE的一些集合变量在CP2K输入文件中的定义:

1.Introduction

​ 关于metadynamics元动力学的理论基础与基本概念,在知乎上面以及有很多人介绍过了,在这里我也简单的介绍一下,首先Mtd本身是属于分子动力学中增强采样的其中一种,也是一种偏置势方法,也就是说当体系自由能处于局部最小值的状态的时候,会导致在你期望的反应坐标空间的采样不充分的问题,比如在蛋白质构象问题的探究上,当陷入自由能的局部极小值中会导致蛋白质的构象采样不充分无法探索到你自己感兴趣的构象状态,所以在metadynamics中可以通过每隔一段时间添加偏置势的方法从而促进体系更早的脱离自由能局部极小值点,这里这个操作在CP2K的**&METADYN模块中被定义为NT_HILLS**,使用者可以在后面添加自定义的整数integer为N,即每隔N个step,便对体系施加一个偏置势,进而更加充分的探索化学空间,从而也就对应上了之前所提到的增强采样的目的,之后在根据计算得到的数据并根据自定义的反应坐标绘制自由能景观图,可以较为直观的观察到整个过程的自由能变化。

​ 个人觉得MtD方法在CP2K中需要自定义的参数是非常多的,也就是说其计算精度在很大层面上会较为依赖于偏置势参数,在CP2K中即为NT_HILLS(定义每隔多少步给体系添加偏置势)HILL_TAIL_CUTOFF(设置此值可以使得高斯尾更快的到达阻尼0值),**&COLVAR(集合变量)模块中的SCALE(可以认为是WIDTH参数项),LAMBDA(指定拓展拉格朗日方案中集体变量与系统坐标耦合的参数),MASS(质量参数)等参数,在这一些参数的定义与选择上针对于不同的体系,往往不同,也就是缺乏普适性,可以被认为是metadynamics这一类方法的不足之一。额外在简单介绍一下,Mtd中,定义的反应坐标的所对应的自由能可以通过如下公式进行计算:

在此公式中
A(z)为自由能项,z即为反应坐标,则是MtD方法中偏置势的表示,有人证明过当模拟时间t趋向于无穷大时,偏置势并不能收敛到一个确定值,并且由于MtD在计算的过程中还需要累积之前计算过的高斯峰,所以说随着高斯峰的数量的增长,模拟每一步的时间都会变慢,从而在Mtd的计算上所需要的资源可能就会略显昂贵和耗时的显著增加,对于CP2K来说更是如此(个人感觉需要50000**步向上才能得到一个较好的结果,如果体系的大小在300-400左右计算的速度就更慢了,如果课题组非常有钱的话也可以当我没说哈哈哈哈)。

​ 目前为止,在CP2K中的**&FREE_ENERGY**模块中可以涉及到的有关自由能计算的软件自带的子模块section主要有下面四个:

  1. ALCHEMICAL_CHANGE
  2. FREE_ENERGY_INFO
  3. METADYN
  4. UMBRELLA_INTEGRATION

目前为止,CP2K在自带的模块中基本只集成了metadynamics这一种较为完善的方法,如果需要使用别的增强采样的算法例如ABF,US等可以将CP2K与PLUMED进行联用,在CP2K的**&FREE_ENERGY**模块中可以设置导入plumed.dat文件,前提是对于plumed的各项设置都较为熟悉。

2.CP2K中的主体部分关键词注解

首先关于集合变量的概念我在我之前的一篇文章中已经提到过了,所以我在此处就暂时不做详细的解释了,在CP2K中&COLVAR模块的定义只能存在于&SUBSYS中,而后续的&FREE_ENERGY模块及其关键词则应该存在于&MOTION之中,**&SUBSYS**模块本身用于定义系统中的几何,电子结构或者原子间的相互作用等信息。在CP2K的输入文件的一般开头我们都会遇到下面几个关键词:

  1. &CELL:用于设置晶胞的一些性质,比如盒子的大小以及体系是否为周期性等等。
  2. ABC:在**&CELL**下会有ABC三个字母单独成行的信息,一般为你的体系的长宽高的设置,由于很多体系需要使用VMD在Z轴或者XY轴方向添加真空层从而避免在模拟过程中镜像的影响,所以为了确保自己添加的真空层是否正确可以在这三行中对应的数值进行检查和判断体系盒子构建的正确性。
  3. PERIODIC:定义体系的周期性,在什么边界条件上产生周期性,可供选择的选项有XYZ,XY,X,Y,Z等。

由于本身CP2K的关键词非常之多,所以后续有时间之后会在上述段落慢慢更新一下几个关键选项的具体含义。

3.CP2K中有关DISTANCE的集合变量的注解

  • &DISTANCE:&DISTANCE应该是CP2K中涉及到DISTANCE距离的集合变量中最容易进行定义的一个,只需要完成定义两个原子的编号以及沿着哪一个方向进行距离计算和考察即可(X | XY | XYZ etc.),关于&DISTANCE的大致输入模板如下所示:
1
2
3
4
5
6
7
&COLVAR
&DISTANCE
AXIS XYZ
ATOMS 1 2
SIGN T
&END DISTANCE
&END COLVAR

关于上述的模板文件中需要对以下的一些变量作以下基本的解释:

  1. AXIS:主要是考察沿着那一个方向进行距离上的计算与考察,在这里如果不是研究特殊的类似于石墨烯的体系仅考虑单一平面上的距离变化的话,建议将AXIS设置为XYZ较为合适。
  2. ATOMS:在后面加上的往往是原子的索引编号,即{integer} {integer}的格式,在这里唯一需要注意的点就是不需要在两个整数之间添加”;” ,”,”等分隔符号。
  3. SIGN:此选项个人感觉可加可不加,如果设置AXIS为单个X,Y,Z的方向则此选项设置了也没有什么用,且SIGN默认为T。
  • &DISTANCE_FUNCTION:&DISTANCE_FUNCTION这个集合变量主要是用于整合两个distance,将其合并为一个distance类型的;集合变量,个人感觉如果在体系中你希望考察多个距离的变化过程,可以采用此CV,从而减少了由于通过第一种方法创建了过多的与距离有关的CV从而导致的模拟问题。其模板文件大致下所示:
1
2
3
4
5
6
7
&COLVAR
&DISTANCE_FUNCTION
ATOMS 1 2 3 4
COEFICIENT -1.0 / 1.0 (another)
PBC T
&END DISTANCE_FUNCTION
&END COLVAR

关于上述的模板文件中我简要的作出如下的注解:

其实这个FUNCTION在公式化之后大致是这样的:

  1. COEFICIENT:所以这里唯一可能需要解释的就是COEFICIENT这一个选项了,但是我个人觉得不需要过多解释,当=1的时候,即考察你所定义的两个距离之和,当=-1的时候则是考察你所定义的两个距离之差。
  2. ATOMS:这里的ATOMS和上面&DISTANCE中的ATOMS存在一定的区别,由于&DISTANCE_FUNCTION主要围绕两个距离建立。可以将此项普适的表示为:ATOMS {integer1} {integer2} {integer3} {integer4},其中则是原子index为 {integer1} 和 {integer2}之间的距离,则可以按照上述逻辑以此类推。
  • &DISTANCE_POINT_PLANE:&DISTANCE_POINT_PLANE用于定义一个点到平面的距离作为一个CV集合变量,概念较为简单,下面则是定义此模块的一个基本示例:
1
2
3
4
5
6
&COLVAR
&DISTANCE_POINT_PLANE
ATOM_PLANE {INTEGER1} {INTEGER2} {INTEGER3}
ATOM_POINT {INTEGER4}
&END DISTANCE_POINT_PLANE
&END COLVAR

上述这样子的模板文件针对与&POINT即&ATOM_POINT较为明确的情况下,如果想要定义一堆原子,将其的质心作为**&POINT**,则可以像下面这样进行定义:

1
2
3
4
5
&POINT
ATOMS {INTEGER1} {INTEGER2} {INTEGER3} {INTEGER4}
TYPE GEO_CENTER
WEIGHTS {REAL1} {REAL2} {REAL3} {REAL4}
&END POINT

基于上述的模板文件,我作出如下注释:

  1. &POINT:这个其实主要用于定义和点相关的性质比如下面会提到的几何中心和权重分配等等
  2. &ATOM_PLANE:用于确定平面PLANE,此处需要定义三个原子编号index,三点确立一个平面的具体位置。
  3. &ATOM_POINT:通过原子编号定义了点是哪一个,以及点的位置在那里。
  4. WEIGHTS:四个原子所占的权重分配,默认的integer数值为1,可以进行自我调节,范围在 ,[0,1] 即可。
  5. TYPE:这里其实主要有两个参数选项,一个是GEO_CENTER,通过所确立四个原子来确定你想要考察的几何中心,当然这里TYPE也可以被改为FIX_POINT,但是我个人觉得没什么必要用这个选项,当然各位如果存在异议的话也欢迎在评论区提问。

最后再摆上一个我最近跑METAD的一个示例文件吧,包括了**&COLVAR&FREE_ENERGY**这两个部分:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
&COLVAR
&DISTANCE
AXIS XYZ
ATOMS 382 381
&END DISTANCE
&END COLVAR

&COLVAR
&DISTANCE
AXIS XYZ
ATOMS 382 380
&END DISTANCE
&END COLVAR

&COLVAR
&DISTANCE
AXIS XYZ
ATOMS 382 377
&END DISTANCE
&END COLVAR

下面的是有关于**&FREE_ENERGY**部分:

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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
&FREE_ENERGY
&METADYN
DO_HILLS T
LANGEVIN
NT_HILLS 40
SLOW_GROWTH
WW 3.0e-3
TEMPERATURE 298.15
TEMP_TOL 10
HILL_TAIL_CUTOFF 2
P_EXPONENT 8
Q_EXPONENT 20
WELL_TEMPERED T
WTGAMMA 4000

&METAVAR
SCALE 0.1
COLVAR 1
&WALL
POSITION 1.0
TYPE QUARTIC
&QUARTIC
DIRECTION WALL_MINUS
K 5.0
&END QUARTIC
&END WALL
&END METAVAR

&METAVAR
SCALE 0.1
COLVAR 2
&WALL
POSITION 1.0
TYPE QUARTIC
&QUARTIC
DIRECTION WALL_MINUS
K 5.0
&END QUARTIC
&END WALL
&END METAVAR

&METAVAR
SCALE 0.1
COLVAR 3
&WALL
POSITION 1.0
TYPE QUARTIC
&QUARTIC
DIRECTION WALL_MINUS
K 5.0
&END QUARTIC
&END WALL
&END METAVAR

&PRINT
&COLVAR
COMMON_ITERATION_LEVELS 3
&EACH
MD 1
&END
&END
&HILLS
COMMON_ITERATION_LEVELS 3
&EACH
MD 1
&END
&END
&END
&END METADYN
&END FREE_ENERGY

如有相关建议和文章中错误的地方:欢迎在评论区评论或者联系我的微信:condainfoenvs :)

  • Title: CP2K中metadynamics输入文件编辑(DISTANCE)
  • Author: Keith
  • Created at : 2024-03-20 20:46:47
  • Updated at : 2024-03-20 20:59:51
  • Link: https://redefine.ohevan.com/2024/03/20/CP2K中metadynamics输入文件编辑-DISTANCE/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments