计算化学公社

标题: 金属-水溶液体系，使用isotropic或semiisotropic控压，要么有真空，要么溶液就飞了 [打印本页]

作者
Author: 牧生 时间: 2021-10-8 16:15
标题: 金属-水溶液体系，使用isotropic或semiisotropic控压，要么有真空，要么溶液就飞了
本帖最后由牧生于 2021-10-8 16:36 编辑

一、建立了α-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模拟，但是此时就有问题了。
(, 下载次数 Times of downloads: 36)   使用isotropic的时候，铁块不变形，说明冻结铁原子是合理的，但就会出现真空区，如图右上角和右下角

(, 下载次数 Times of downloads: 35) 使用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

作者
Author: sobereva 时间: 2021-10-9 17:49
用semiisotropic，Fe用位置限制固定而不要用freeze，否则结合控压的时候容易出问题。平行于界面方向的可压缩系数设为0

作者
Author: supersix 时间: 2021-10-10 16:36
我是这么控压的，金属原子也是freeze，好像没遇到这些问题
pcoupl = Berendsen
pcoupltype = semiisotropic
tau_p = 5.0
ref_p = 1.0 1.0
compressibility = 0 4.5e-5

作者
Author: 牧生 时间: 2021-10-10 17:03
本帖最后由牧生于 2021-10-10 17:40 编辑

supersix 发表于 2021-10-10 16:36
我是这么控压的，金属原子也是freeze，好像没遇到这些问题
pcoupl = Berendsen
pcoupltype = semiiso ...

因为铁的参数比较特殊，目前文献中尚无α铁晶胞的参数。文献中只给出了γ-铁的参数，如果直接使用γ-铁的晶胞进行模拟，那么，既不需要冻结，也不需要限制，金属部分也丝毫不变，一切一切都很简单，操作都是常规操作，结果也都非常符合预期。

但是如果使用α-铁晶胞，文献中说明了，γ-铁的参数也是可用的，但需要将铁原子fix。我试过sob说的限制，用了gmx genrestr -f FE.gro -o posre_FE.itp，将铁限制了，且在top中引用了该限制文件itp，但是，发现NPT以后，铁胞还是有少许变形。如果使用冻结，那么，α-铁倒是一点也不会变形，但是水溶液区域就有真空了。感觉限制不如冻结那么严格，或者我的操作有误，目前正在考虑。
所以，我现在就在考虑，文献中那个fix，到底应该是冻结，还是应该限制？？我自己觉得应该是冻结比较合适。好像我距离成功很近了，但就是还差一点。我明天去服务器上试试你这个参数。

作者
Author: 牧生 时间: 2021-11-4 10:56
本帖最后由牧生于 2021-11-11 08:43 编辑

2021.11.4补充：
即使使用solvate -maxsol 命令，虽然看起来加满水，但npt以后，自发调整了，还是会有真空。
现在的解决办法：
在npt以后有真空区的盒子里面，再用solvate命令加水，再em，再npt。反复操作几次，这样盒子里面基本就正常了，也没有明显真空了。似乎总有哪里不对。
在我反复做了好多次的结果上，我始终认为冻结比限制好。。冻结可以保持铁原子纹丝不动，而限制，铁块始终会轻微变形。

2021.11.8补充，使用冻结命令，虽然铁不变形，但是长时间MD，盒子会变化，溶液区的大小也变化，从而还是出现真空。

作者
Author: nianbin 时间: 2021-11-7 18:21
你Fe的晶格参数以及力场参数贴上来

作者
Author: 牧生 时间: 2021-11-7 19:20

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

作者
Author: 牧生 时间: 2021-11-11 09:07
本帖最后由牧生于 2021-11-11 16:52 编辑

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以后，铁块还是有很轻微的变形

作者
Author: xiaoruoqiu 时间: 2022-8-17 17:13

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

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

作者
Author: cheviax 时间: 2023-7-1 17:51

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

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

作者
Author: chuxuedexiaobai 时间: 2024-4-22 20:34

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

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

作者
Author: 牧生 时间: 2024-4-22 21:42

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

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

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

作者
Author: 喵星大佬 时间: 2024-4-22 21:50

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

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

作者
Author: chuxuedexiaobai 时间: 2024-4-23 15:43

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

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

作者
Author: chuxuedexiaobai 时间: 2024-4-23 15:45

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

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

作者
Author: chuxuedexiaobai 时间: 2024-4-23 15:46

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

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

作者
Author: 牧生 时间: 2024-4-23 15:59
本帖最后由牧生于 2024-4-23 16:13 编辑

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

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

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

作者
Author: 喵星大佬 时间: 2024-4-23 17:54

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

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

作者
Author: chuxuedexiaobai 时间: 2024-4-24 15:11

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

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

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