gromacs在有机分子-磷脂DPPC的体系做伞形采样过程出现异常

chengyuseu

各位老师好，我参考之前的教程，想自己做一个有机小分子在穿模过程的一个伞形采样，正常运行到了构型生成这一步（对应于教程中的第五步），但是当我检查summary_distances.dat中帧和对应有机分子与DPPC之间的距离时发现，会出现多个平台。我想问一下各位老师，这样的平台的出现是不是合理的？距离出现这样大的突越是否合理？还有就是我给的md_pull.mdp是不是合理的呢？谢谢各位老师。



md_pull.mdp如下：
title                = Umbrella pulling simulation
define               = -DPOSRES  ; 对DPPC限制，不对有机分子限制

; Run control
integrator             = md     
nsteps                     = 5000000  
dt                     = 0.002   
tinit                = 0        
nstcomm              = 10     

; Output control
nstxout              = 5000   
nstvout              = 5000     
nstfout              = 5000     
nstenergy            = 5000   
nstxout-compressed   = 5000   
nstlog                     = 5000  

; Bonds
constraint_algorithm = lincs   
constraints          = all-bonds
continuation         = yes      

; Neighbor searching
cutoff-scheme             = Verlet
nstlist              = 10   
ns_type              = grid
rlist                = 1.4
pbc                  = xyz

; Electrostatics
coulombtype          = PME
rcoulomb             = 1.4

; Van der Waals
rvdw                 = 1.4
DispCorr             = EnerPres

; Ewald
fourierspacing       = 0.12
fourier_nx           = 0
fourier_ny           = 0
fourier_nz           = 0
pme_order            = 4
ewald_rtol           = 1e-5

; Temperature coupling
tcoupl               = Nose-Hoover
tc_grps              = DPP_KDZ Water_and_ions
tau_t                = 0.5       0.5
ref_t                = 323       323

; Pressure coupling
pcoupl               = Parrinello-Rahman
pcoupltype           = semiisotropic
tau_p                = 5.0      
ref_p                     = 1.0        1.0                  
compressibility      = 4.5e-5        4.5e-5        
refcoord_scaling     = com

; Velocity generation
gen_vel              = no

; COM pulling
pull                 = yes
pull_ngroups         = 2
pull_ncoords         = 1
pull_group1_name     = KDZ
pull_group2_name     = DPP
pull_coord1_type     = umbrella      ; harmonic biasing force
pull_coord1_geometry = distance      ; simple distance increase
pull_coord1_groups   = 1 2
pull_coord1_dim      = N N Y
pull_coord1_rate     = -0.001          ; 0.001 nm per ps = 1 nm per ns
pull_coord1_k        = 1000          ; kJ mol^-1 nm^-2
pull_coord1_start    = yes           ; define initial COM distance > 0

wyb

查看一下轨迹文件，是否跨过了周期性边界盒子

chengyuseu

麻烦各位老师看下，谢谢。

chengyuseu

wyb 发表于 2020-7-17 08:43
查看一下轨迹文件，是否跨过了周期性边界盒子

谢谢老师的回复，我检查了第500帧和第1000帧的gro文件，体系中的有机小分子都在体系的中间。感觉500到1000帧之间，有机小分子没有被牵引，没有移动。老师，是不是我设置的参数有问题。

wyb

但你的图片是显示400到500之间出现了问题，你再看一下，不过按你的描述，500-1000小分子没动那肯定是你没有设置好体系

chengyuseu

wyb 发表于 2020-7-20 09:08
但你的图片是显示400到500之间出现了问题，你再看一下，不过按你的描述，500-1000小分子没动那肯定是你没有 ...

好的，谢谢老师，我再检查下。

byqyb

chengyuseu 发表于 2020-7-22 14:36
好的，谢谢老师，我再检查下。

您好，想咨询一下您伞取样的问题，我也是做有机分子穿过磷脂双分子层。我是在umbrella-sampling过程中提示When the maximum distance from a pull group reference atom to other atoms in the group is larger than 0.5 times half the box size a centrally placed atom should be chosen as pbcatom. Pull group 1 is larger than that and does not have a specific atom selected as reference atom.但实际上我的box长度一半远大于间隔距离。
以下是我的pull code：
; Pull code
pull                    = yes
pull_ncoords            = 1         ; only one reaction coordinate
pull_ngroups            = 2         ; two groups defining one reaction coordinate
pull_group1_name        = POPC
pull_group2_name        = DMP
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
想问一下您这个该怎么解决呢？或者我能参考一下您的pull code吗？

Kamala

这教程  楼主可以分享以下吗

黄舒伟

byqyb 发表于 2021-5-11 16:29
您好，想咨询一下您伞取样的问题，我也是做有机分子穿过磷脂双分子层。我是在umbrella-sampling过程中提 ...

可以参考：https://www.researchgate.net/pos ... ion_in_GROMACS_2020

byqyb

黄舒伟 发表于 2022-6-2 16:27
可以参考：https://www.researchgate.net/pos ... ion_in_GROMACS_2020

好的，谢谢你的答复。这个问题我之前已经解决了。下面是我修改过的pull code
; Pull code
pull                    = yes
pull_ncoords            = 1         ; only one reaction coordinate
pull_ngroups            = 2         ; two groups defining one reaction coordinate
pull_group1_name        = dehp_benzan
pull_group2_name        = POPC
pull_group2_pbcatom     = 8576
pull-pbc-ref-prev-step-com = yes
pull_coord1_k           = 1000      ; kJ mol^-1 nm^-2
pull-coord1-vec      = 0.0 0.0 -1.0
pull_coord1_type     = umbrella
pull_coord1_geometry = direction
pull_coord1_groups   = 1 2
pull_coord1_dim      = N N Y
pull_coord1_rate     = -0.0001      ; 0.0001 nm per ps = 0.1 nm per ns
pull_coord1_start    = yes

dimoo

请问楼主为什么要添加如图所示的pull参数

tangzi_via

同问楼主，这个8576是你第二个组的原子数吗？

NIiko

dimoo 发表于 2022-10-5 16:34
请问楼主为什么要添加如图所示的pull参数

当不指定此选项时，gromacs 只会从列表中搜索群中的中心原子，即，如果您有一组 10 个原子，则 gromacs 使用第五个原子，这不一定与您的群的几何中心相对应。这通常不是问题，但是当组大于盒子的一半时，使用 PBC 计算距离的选择可能并不理想，因此您必须指定一个更好的原子供 gromacs 使用。