计算化学公社

标题: Gromacs如何实现SMD？ [打印本页]

作者
Author: FrancisCho 时间: 2024-6-13 19:21
标题: Gromacs如何实现SMD？
本帖最后由 FrancisCho 于 2024-6-13 19:21 编辑

GMX的mdp指令中如何实现对一个原子/组施加一个力，使其穿过碳纳米管？
我在复刻完GMX官方的Umbrella Sampling教程（http://www.mdtutorials.com/gmx/umbrella/05_pull.html）后想将其迁移至自己研究的体系中。
针对教程中给的相关文件，我对其进行修改，在md_pull.mdp中使用如下指令：
; Pull code
pull                   = yes
pull_ncoords          = 1       ; only one reaction coordinate
pull_ngroups          = 2       ; two groups defining one reaction coordinate
pull_group1_name       = NA
pull_group2_name       = CNT
pull_coord1_type       = umbrella  ; harmonic potential
pull_coord1_geometry = distance  ; simple distance increase
pull_coord1_dim       = N N Y
pull_coord1_groups    = 1 2
pull_coord1_start    = yes    ; define initial COM distance > 0
pull_coord1_rate       = 0.01    ; 0.01 nm per ps = 10 nm per ns
pull_coord1_k          = 20    ; kJ mol^-1 nm^-2

进行500ps模拟后使用VMD观看轨迹，发现钠原子并不往碳纳米管方向移动，而是沿着反方向移动。
使用 pull_coord1_type  = constant-force 后，钠原子只会被拉动至碳纳米管的质心处，而不会穿过整个碳纳米管。
在lammps中，使用fix smd cation smd cvel 5 0.001 tether NULL NULL -30 0.0 就可以让钠离子穿过碳纳米管，我想请教一下：
1、在GMX中如何实现拉动（最好是umbrella）钠离子穿过碳纳米管？
2、GMX中的pull指令总是会因为拉动距离超过盒子对应方向距离的0.49倍而报错，这一点除使用gmx editconf -f CNT10x10_3.gro -o CNT_box.gro扩大盒子尺寸外有解决方法吗？

另外，在使用pull_coord1_type  = constant-force后，钠离子被拉动至碳纳米管质心附近，使用gmx trjconv -s pull.tpr -f pull.xtc -o conf.gro -sep从轨迹提取帧后使用：title    = Umbrella pulling simulation
define    = -DPOSRES_B
; Run parameters
integrator  = md
dt       = 0.002
tinit    = 0
nsteps    = 50000    ; 100 ps
nstcomm    = 10
; Output parameters
nstvout          = 5000  ; every 10 ps
nstfout          = 5000
nstxout-compressed  = 5000
nstenergy          = 5000
; Bond parameters
constraint_algorithm = lincs
constraints          = all-bonds
continuation          = no
; Single-range cutoff scheme
cutoff-scheme = Verlet
nstlist       = 20
ns_type       = grid
rlist          = 1.4
rcoulomb       = 1.4
rvdw          = 1.4
; PME electrostatics parameters
coulombtype    = PME
fourierspacing  = 0.12
fourier_nx    = 0
fourier_ny    = 0
fourier_nz    = 0
pme_order    = 4
ewald_rtol    = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups
Tcoupl    = Berendsen
tc_grps    = Protein Non-Protein
tau_t    = 0.5    0.5
ref_t    = 310    310
; Pressure coupling is on
;Pcoupl       = Berendsen
;pcoupltype    = isotropic
;tau_p          = 1.0
;compressibility = 4.5e-5
;ref_p          = 1.0
;refcoord_scaling = com
freezegrps = CNT
freezedim = N N Y
; Generate velocities is on
gen_vel    = yes
gen_temp = 310
; Periodic boundary conditions are on in all directions
pbc    = xyz
; Long-range dispersion correction
DispCorr = EnerPres
; Pull code
pull                   = yes
pull_ncoords          = 1       ; only one reaction coordinate
pull_ngroups          = 2       ; two groups defining one reaction coordinate
pull_group1_name       = NA
pull_group2_name       = CNT
pull_coord1_type       = umbrella  ; harmonic potential
pull_coord1_geometry = distance  ; simple distance increase
pull_coord1_dim       = N N Y
pull_coord1_groups    = 1 2
pull_coord1_start    = yes    ; define initial COM distance > 0
pull_coord1_rate       = 0.0    ; restrain in place
pull_coord1_k          = 1000    ; kJ mol^-1 nm^-2

进行NPT，然后使用：
title    = Umbrella pulling simulation
; Run parameters
integrator  = md
dt       = 0.002
tinit    = 0
nsteps    = 5000000 ; 10 ns
nstcomm    = 10
; Output parameters
nstxout-compressed  = 5000    ; every 10 ps
nstenergy          = 5000
; Bond parameters
constraint_algorithm = lincs
constraints          = all-bonds
continuation          = yes
; Single-range cutoff scheme
cutoff-scheme = Verlet
nstlist       = 20
ns_type       = grid
rlist          = 1.4
rcoulomb       = 1.4
rvdw          = 1.4
; PME electrostatics parameters
coulombtype    = PME
fourierspacing  = 0.12
fourier_nx    = 0
fourier_ny    = 0
fourier_nz    = 0
pme_order    = 4
ewald_rtol    = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups
Tcoupl    = Nose-Hoover
tc_grps    = Protein
tau_t    = 1.0
ref_t    = 310
; Pressure coupling is on
;Pcoupl       = Parrinello-Rahman
;pcoupltype    = isotropic
;tau_p          = 1.0
;compressibility = 4.5e-5
;ref_p          = 1.0
;refcoord_scaling = com
freezegrps = CNT
freezedim = N N Y
; Generate velocities is off
gen_vel    = no
; Periodic boundary conditions are on in all directions
pbc    = xyz
; Long-range dispersion correction
DispCorr = EnerPres

; Pull code
pull                   = yes
pull_ncoords          = 1       ; only one reaction coordinate
pull_ngroups          = 2       ; two groups defining one reaction coordinate
pull_group1_name       = NA
pull_group2_name       = CNT
pull_coord1_type       = umbrella  ; harmonic potential
pull_coord1_geometry = distance  ; simple distance increase
pull_coord1_dim       = N N Y
pull_coord1_groups    = 1 2
pull_coord1_start    = yes    ; define initial COM distance > 0
pull_coord1_rate       = 0.0    ; restrain in place
pull_coord1_k          = 1000    ; kJ mol^-1 nm^-2

使用gmx wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kCal进行伞式采样，发现PMF曲线与预期不符
查看histo.xvg曲线发现并没有交联。
请问该如何解决？

作者
Author: lucky1999 时间: 2024-6-14 08:38
若需要通过gmx实现SMD，个人建议利用plumed与gmx联动完成SMD，这样的例子可以在plumed的官网上可以找到(https://www.plumed.org/doc-v2.8/user-doc/html/belfast-5.html)。gmx对粒子进行牵引有某些限制(例如，牵引距离需要小于盒子牵引方向所在轴长度的二分之一)，个人感觉不是很好用，而若使用plumed则无这样的限制。

作者
Author: Dempey 时间: 2024-6-14 08:59
1. GROMACS的pull我一般设置是

pull-coord1-type=umbrella
pull-coord1-geometry=direction
pull-coord1-vec= 0 0 1

复制代码

geometry用direction可以指定拉动方向的矢量，很好用
2. 要想拉动距离超过盒子最小长的一半，geometry设置为direction-periodic，但是不能控压
3. 拉动产生多个采样窗口时可以调小拉动速度(pull-coord1-rate)，增大拉力(pull-coord1-k)，这样产生的构型比较均匀
4. 你的采样很不充分，可以减小采样窗口的距离或者减小采样时的拉力(pull-coord1-k)，前者让histo.xvg中的曲线更密，后者让曲线更胖
用GROMACS自己做伞形采样要不断尝试，优化出一个比较好的参数，同时也很消耗资源，也建议你使用PLUMED或COLVARS做自由能计算，很省事

作者
Author: FrancisCho 时间: 2024-6-14 09:17

lucky1999 发表于 2024-6-14 08:38
若需要通过gmx实现SMD，个人建议利用plumed与gmx联动完成SMD，这样的例子可以在plumed的官网上可以找到(htt ...

非常感谢，我之前在PLUMED上看了不少有关gmx结合PLUMED算PMF的教程，但是GMX安装PLUMED需要在GMX安装之前先编译PLUMED，后来使用plumed driver去读取gmx轨迹文件，但是运行-ixtc报错。我现在在尝试用LAMMPS+Colvars计算PMF，但是我发现在运行Colvars前使用
fix 1 pore setforce 0.0 0.0 0.0
来固定碳纳米管，但是在Colvars运行过程中整个碳纳米管变形，并且移动，想向您请教一下

作者
Author: lucky1999 时间: 2024-6-14 09:37

FrancisCho 发表于 2024-6-14 09:17
非常感谢，我之前在PLUMED上看了不少有关gmx结合PLUMED算PMF的教程，但是GMX安装PLUMED需要在GMX安装之前 ...

对于前一个问题我可以试着跟你解答一下，对于后一个问题，由于我较少使用LAMMPS和Colvars，所以无法给你解答。这里应该用不到plumed中的driver -ixtc命令吧。drive -ixtc命令是你已经通过gmx得到了xtc文件，而后续使用plumed对得到的xtc文件进行相关处理。显然，你这里需要在gmx模拟的过程中就调用plumed进行SMD。对于此，你需要写一个plumed的输入文件(文件里需要写啥可以看我上面给你发的那个网站)，然后在使用gmx mdrun命令时采用-plumed参数指定plumed输入文件。这样，你就可以在gmx模拟的过程中同时调用plumed进行SMD了。

作者
Author: FrancisCho 时间: 2024-6-14 09:46

lucky1999 发表于 2024-6-14 09:37
对于前一个问题我可以试着跟你解答一下，对于后一个问题，由于我较少使用LAMMPS和Colvars，所以无法给你 ...

我想请教一下如果我先安装Gromacs，后安装plumed，Gromacs是不是无法被后安装的plumed patched？

作者
Author: lucky1999 时间: 2024-6-14 09:55

FrancisCho 发表于 2024-6-14 09:46
我想请教一下如果我先安装Gromacs，后安装plumed，Gromacs是不是无法被后安装的plumed patched？

得使用打了plumed补丁的gmx，这样才能在gmx模拟的过程中调用plumed。否则，plumed只能用于gmx模拟完成后的后处理操作。

作者
Author: Daniel_Arndt 时间: 2024-6-14 10:56

FrancisCho 发表于 2024-6-14 09:46
我想请教一下如果我先安装Gromacs，后安装plumed，Gromacs是不是无法被后安装的plumed patched？

plumed的Documentation中会提及适用的gromacs，如plumed 2.9适用的gromacs可以看 https://www.plumed.org/doc-v2.9/user-doc/html/index.html#codes 。正常的流程是先编译plumed，再对gromacs代码进行patch，最后编译patch过的gromacs，可以看 https://www.plumed.org/doc-v2.9/user-doc/html/_installation.html 。

作者
Author: FrancisCho 时间: 2024-6-14 14:45

Dempey 发表于 2024-6-14 08:59
1. GROMACS的pull我一般设置是

geometry用direction可以指定拉动方向的矢量，很好用

您好，我能请教一下有使用Colvars或PLUMED结合GMX或LAMMPS计算PMF的教程实例吗？最好是有给出相应代码的，谢谢。

作者
Author: Dempey 时间: 2024-6-14 15:33

FrancisCho 发表于 2024-6-14 14:45
您好，我能请教一下有使用Colvars或PLUMED结合GMX或LAMMPS计算PMF的教程实例吗？最好是有给出相应代码的 ...

PLUMED官网就有：https://www.plumed.org/doc-v2.8/user-doc/html/_m_e_t_a_d.html
COLVARS公社就有：http://bbs.keinsci.com/thread-10805-1-1.html

作者
Author: QYQ 时间: 2024-10-8 20:06
兄弟，你的问题和我的有点相似，请问你有解决方案了吗。我的是pull_coord1_type = constant-force可以把小分子拉出来，而且很快，但改成Umbrella不行了或者是往蛋白的另一端靠，其它参数都没变

作者
Author: QYQ 时间: 2024-10-8 20:10

Dempey 发表于 2024-6-14 08:59
1. GROMACS的pull我一般设置是

geometry用direction可以指定拉动方向的矢量，很好用

兄弟，我的设置一开始是：
pull-coord1-type=constant-force
pull-coord1-geometry=direction
pull-coord1-vec= 0 0 1
pull_coord1_dim = N N Y
后面把力的类型改了：
pull-coord1-type=umbrella
pull-coord1-geometry=direction
pull-coord1-vec= 0 0 1
pull_coord1_dim = N N Y
但为啥小分子移动的方向变了，前面设置的是可以解离出去的，换伞形后配体往蛋白中心走了
请问兄台你遇到过这种情况吗

欢迎光临计算化学公社 (http://bbs.keinsci.com/)