金属-水溶液体系，使用isotropic或semiisotropic控压，要么有真空，要么溶液就飞了

牧生

一、建立了α-Fe晶胞

二、扩展盒子为长方体，将α-Fe置于盒子中间（铁参数来自于10.1038/s41524-020-00478-1，文献中特别说明，使用α-Fe进行模拟时，atomic positions should be fixed）

三、使用insert命令，向空白部分加入50个表面活性剂分子（使用ATB服务器得到的itp）

四、使用solvate命令加满水，加入50个CL离子

五、能量最小化
define = -DFLEXIBLE
integrator = steep
nsteps = 10000
emtol  = 100.0
emstep = 0.01
;
nstxout   = 100
nstlog    = 50
nstenergy = 50
;
pbc = xyz
cutoff-scheme            = Verlet
coulombtype              = PME
rcoulomb                 = 1.0
vdwtype                  = Cut-off
rvdw                        = 1
DispCorr                 = EnerPres
;
constraints              = none
freezegrps              = FE ；文献中提到的atomic positions should be fixed，我就在这里冻结了所有的铁原子。。在itp里面的参数，都是大写的FE。 如果这一步不冻结铁原子，那么em以后，铁块就有点变形。
freezedim               = Y Y Y

em可以正常进行，em结束后打开得到的gro文件看起来也正常，vmd没有提示不合理的接触，溶液区没有乱跑，铁块完整不变形。


六、进行NPT模拟，但是此时就有问题了。
  使用isotropic的时候，铁块不变形，说明冻结铁原子是合理的，但就会出现真空区，如图右上角和右下角


使用semiisotropic的时候，右边一半的溶液就跑到了很远很远的地方。



NPT内容如下：

define =
integrator = md
dt         = 0.002  ; ps
nsteps     = 50000000
comm-grps  = system
energygrps =
;
nstxout = 0
nstvout = 0
nstfout = 0
nstlog  =
nstenergy = 50000
nstxout-compressed = 50000
compressed-x-grps  = system
;
annealing = single single
annealing_npoints = 2 2
annealing_time = 0 100 0 100
annealing_temp = 0 298.15 0 298.15
;
pbc = xyz
cutoff-scheme = Verlet
coulombtype   = PME
rcoulomb      = 1
vdwtype       = cut-off
rvdw          = 1
DispCorr      = EnerPres
;
Tcoupl  = V-rescale
tau_t   = 0.2 0.2
tc_grps = SOL non-water
ref_t   = 298.15 298.15
;
Pcoupl     = Berendsen
pcoupltype = isotropic     ; 用 isotropic或者semiisotropic，结果不一样，我可以理解。但是isotropic进行npt模拟，出现真空区，我就不理解了。
tau_p = 0.5
ref_p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5
;
gen_vel  = no
gen_temp = 298.15
gen_seed = -1
;
freezegrps  = FE            
freezedim   = Y Y Y
constraints = hbonds

chuxuedexiaobai

牧生 发表于 2024-4-23 15:59
目前情况下，用gmx跑的话，除非你把铁原子都冻结，才能跑出一个符合预期的合理结果。这样做的话，潜在风 ...

感觉要换课题了，确实我试了好多次，铁面就是会变形，哎，感谢您的解答

喵星大佬

chuxuedexiaobai 发表于 2024-4-23 15:43
老师，如果要用铁的话，跑的话要用什么跑才行呢

gmx就没法做，换力场换软件

牧生

chuxuedexiaobai 发表于 2024-4-23 15:46
如何确定自己的金属参数是否适合用来模拟呀

目前情况下，用gmx跑的话，除非你把铁原子都冻结，才能跑出一个符合预期的合理结果。这样做的话，潜在风险是被质疑冻结的合理性。

至于确定参数的合理性，跑一下MD，铁胞变不变形等，自然就看出来了。

chuxuedexiaobai

牧生 发表于 2024-4-22 21:42
最终实际没解决。。因为文献中的α铁参数不完美，需要用限制才行，但限制后仍然有轻微变形，冻结的话倒是 ...

如何确定自己的金属参数是否适合用来模拟呀

chuxuedexiaobai

牧生 发表于 2024-4-22 21:42
最终实际没解决。。因为文献中的α铁参数不完美，需要用限制才行，但限制后仍然有轻微变形，冻结的话倒是 ...

我现在就是这个铁也是一直有问题导致我课题干不下去，但是我看那个论文里面他就是用的铁但是他是110切面，是不是说切面就可以用呢

chuxuedexiaobai

喵星大佬 发表于 2024-4-22 21:50
α铁是BCC的晶胞，用LJ势绝对不可能做到稳定晶胞是BCC的，会自发变成FCC的，所以就不该用LJ势跑这种东西

老师，如果要用铁的话，跑的话要用什么跑才行呢

喵星大佬

牧生 发表于 2024-4-22 21:42
最终实际没解决。。因为文献中的α铁参数不完美，需要用限制才行，但限制后仍然有轻微变形，冻结的话倒是 ...

α铁是BCC的晶胞，用LJ势绝对不可能做到稳定晶胞是BCC的，会自发变成FCC的，所以就不该用LJ势跑这种东西

牧生

chuxuedexiaobai 发表于 2024-4-22 20:34
想请教一下楼主，您铁面固定问题解决了吗，我用MS直接构建了铁的超晶胞结构然后用位置限制跑完模拟铁就散 ...

最终实际没解决。。因为文献中的α铁参数不完美，需要用限制才行，但限制后仍然有轻微变形，冻结的话倒是纹丝不动。我觉得用了限制或冻结的话，那么参数的意义就不是很大了，所以我放弃α铁的模拟了。。

但如果做别的有完美参数的金属，比如金，铝板，完全不用任何限制或者冻结，模拟结果就很符合预期。

chuxuedexiaobai

牧生 发表于 2021-11-11 09:07
已经解决。

想请教一下楼主，您铁面固定问题解决了吗，我用MS直接构建了铁的超晶胞结构然后用位置限制跑完模拟铁就散乱了，我以为是因为我只写了一个铁原子限制有问题然后位置限制文件将所以有的铁原子都写进去但也出现了你这同样的问题，想问最后这个铁板固定您是如何解决的呢

cheviax

xiaoruoqiu 发表于 2022-8-17 17:13
有个笨办法供参考，搭建体系时候先分三个小的体系（即水、铁、水）做npt，xy向压缩度设为0。npt平衡之后 ...

水在铁表面会富集，分开NPT然后拼起来的话溶液密度整体还是会偏小

xiaoruoqiu

牧生 发表于 2021-11-11 09:07
已经解决。

有个笨办法供参考，搭建体系时候先分三个小的体系（即水、铁、水）做npt，xy向压缩度设为0。npt平衡之后再拼接起来，搭建好最终体系。

牧生

sobereva 发表于 2021-10-9 17:49
用semiisotropic，Fe用位置限制固定而不要用freeze，否则结合控压的时候容易出问题。平行于界面方向的可压 ...

已经解决。



使用gmx genrestr -f FE.pdb -o posre_FE.itp得到了铁块限制的itp文件，FE.pdb原子编号为1-8000，gro里面，FE的编号也是从1-8000，得到posre_FE.itp编号也是1-8000

然后在top里面使用了该限制itp，内容为：
#include "forcefield.itp"
#include "FE.itp"
; Include Position restraint file
#ifdef POSRES         //问题出在这里，改成#ifdef POSRES_FE就不报错了   ，但知其然，不知其所以然，仍需要指导
#include "posre_FE.itp"
#endif

在npt.mdp里面，第一行写上
define = -DPOSRES

结果报错
ERROR 1 [file posre_FE.itp, line 6]:
  Atom index (2) in position_restraints out of bounds (1-1).
  This probably means that you have inserted topology section
  "position_restraints"
  in a part belonging to a different molecule than you intended to.
  In that case move the "position_restraints" section to the right molecule.

npt以后，铁块还是有很轻微的变形

牧生

nianbin 发表于 2021-11-7 18:21
你Fe的晶格参数以及力场参数贴上来

铁参数来自于10.1038/s41524-020-00478-1，文献中特别说明，使用α-Fe进行模拟时，atomic positions should be fixed

针对gromos力场
FE      26    0            0      A   0.01514608     2.2873928e-06


针对opls力场和amber力场
FE              FE           26     55.845        0.0000    A      0.230742768      25.104

nianbin

你Fe的晶格参数以及力场参数贴上来