LAMMPS分析固液体系中溶液离子MSD起始时刻MSD不为0的问题

Mingxin

体系是上层为溶液，下层基底为固体，利用LAMMPS想要计算溶液中离子自运行100wstep后的MSD，因此lammps中的输入文件为运行100w步后的restart.文件，LAMMPS manual关于msd命令说是“The displacement of an atom is from its reference position. This is normally the original position at the time the compute command was issued, unless the average keyword is set to yes. ”将命令发出时刻的原子位置作为起始参考位置，所以MSD应该起始时刻应该是0，但我计算后的结果起始时刻的MSD不为零，并且此时的MSD数值对应的是把最初0step时的位置作为起始参考位置所得到的值，而不是把我输入的100wstep文件作为起始位置，想问下大家是什么原因呢?
以下是我的lammps code，我是将溶液中的离子按照距离表面的距离分到了不同的group里分别计算MSD，最后再附上其中一个group的MSD计算结果图:

1. #Basic setup
   
2. units                real
   
3. boundary             p p p
   
4. atom_style           full
   
5. bond_style           harmonic
   
6. angle_style          harmonic
   
7. pair_style           lj/cut/coul/long 12.5
   
8. timestep             1   
   
9. neighbor             2 bin
   
10. neigh_modify         every 2 delay 0 check yes one 3000
    
11. 
    
12. #Structure
    
13. read_restart         restart.1000000
    
14. pair_modify          mix arithmetic
    
15. kspace_style         ewald  1.0e-4   
    
16. 
    
17. 
    
18. #Solution atom group
    
19. group                Na  type 5
    
20. group                O*  type 7
    
21. group                h*  type 8
    
22. group                Cu  type 10
    
23. group                Cl  type 9
    
24. 
    
25. #N-A-S-H group
    
26. group                Al  type 2
    
27. group                Si  type 1
    
28. group                O   type 3 4
    
29. group                H   type 6                 
    
30. group                Geopolymer union Al Si O H Na
    
31. group                Solution subtract all Geopolymer
    
32. 
    
33. #divide the system into surface and solution part
    
34. region              surface_part block INF INF INF INF 22 30 units box
    
35. group               Surface_atom region surface_part
    
36. 
    
37. region              solution_part_1 block INF INF INF INF 30 35 units box
    
38. group               Solution_atom_1 region solution_part_1
    
39. 
    
40. region              solution_part_2 block INF INF INF INF 35 40 units box
    
41. group               Solution_atom_2 region solution_part_2
    
42. 
    
43. region              solution_part_3 block INF INF INF INF 40 45 units box
    
44. group               Solution_atom_3 region solution_part_3
    
45. 
    
46. region              solution_part_4 block INF INF INF INF 45 50 units box
    
47. group               Solution_atom_4 region solution_part_4
    
48. 
    
49. #puts the Cu atoms into different part
    
50. group               Cu_surface intersect Surface_atom Cu
    
51. group               Cu_solution_1 intersect Solution_atom_1 Cu
    
52. group               Cu_solution_2 intersect Solution_atom_2 Cu
    
53. group               Cu_solution_3 intersect Solution_atom_3 Cu
    
54. group               Cu_solution_4 intersect Solution_atom_4 Cu
    
55. 
    
56. #initial temperature
    
57. velocity             all create 300 39849 mom yes rot yes dist gaussian
    
58. 
    
59. #Thermo outputs
    
60. thermo               100
    
61. thermo_style         custom step temp pe ke
    
62. 
    
63. #rising temperature
    
64. fix                  1 all nvt temp 300 300 1
    
65. 
    
66. #计算Cu离子的surface part MSD
    
67. #modified code
    
68. compute              MSD1 Cu_surface msd com yes
    
69. fix                  2 Cu_surface ave/time 100 5 500 c_MSD1[1] c_MSD1[2] c_MSD1[3] c_MSD1[4] file Cu_msd_surface.data mode scalar
    
70. 
    
71. #计算Cu离子的solution part MSD
    
72. compute              MSD2 Cu_solution_1 msd com yes
    
73. fix                  3 Cu_solution_1 ave/time 100 5 500 c_MSD2[1] c_MSD2[2] c_MSD2[3] c_MSD2[4] file Cu_msd_solution_1.data mode scalar
    
74. 
    
75. #计算Cu离子的solution part MSD
    
76. compute              MSD3 Cu_solution_2 msd com yes
    
77. fix                  4 Cu_solution_2 ave/time 100 5 500 c_MSD3[1] c_MSD3[2] c_MSD3[3] c_MSD3[4] file Cu_msd_solution_2.data mode scalar
    
78. 
    
79. #计算Cu离子的solution part MSD
    
80. compute              MSD4 Cu_solution_3 msd com yes
    
81. fix                  5 Cu_solution_3 ave/time 100 5 500 c_MSD4[1] c_MSD4[2] c_MSD4[3] c_MSD4[4] file Cu_msd_solution_3.data mode scalar
    
82. 
    
83. #计算Cu离子的solution part MSD
    
84. compute              MSD5 Cu_solution_4 msd com yes
    
85. fix                  6 Cu_solution_4 ave/time 100 5 500 c_MSD5[1] c_MSD5[2] c_MSD5[3] c_MSD5[4] file Cu_msd_solution_4.data mode scalar
    
86. 
    
87. #dump
    
88. dump                 1 all custom 100000 A1.sputter id type x y z
    
89. run                  1000000
    
90. write_restart        restart.2000000
    

复制代码

Graphite

restart文件只是为了中断续算的，里面的step定义是比较模糊的，特别是一次模拟中进行了多个过程时，建议用data文件。
最靠谱的算各种轨迹分析/能量分析...的方法是：
1、正常做一次模拟，读入data文件而不是restart文件
2、跑的时候不要compute MSD或者compute RDF等等分析，而是较高频次地输出轨迹
3、跑完了把轨迹rerun一遍，这不会动力学走步，而是对每步重算能量和结构、以及各种分析

假如500步输出一次轨迹，那么事后rerun的基础时间成本仅仅是1/500。加了各种分析之后，也只是把该耗的那些时间补上。
这样做的好处是简单化，避免“路径依赖”的歧义，也避免因为各种分析中途出错而导致耗时很长的主采样挂掉。

Mingxin

Graphite 发表于 2023-9-13 15:23
restart文件只是为了中断续算的，里面的step定义是比较模糊的，特别是一次模拟中进行了多个过程时，建议用d ...

谢谢您的解答，看来之前是对restart文件的应用理解有偏差，更感谢您建议的分析流程！我之前的分析过程确实不规范，学习到了！

zhaoxiaollong

Graphite 发表于 2023-9-13 15:23
restart文件只是为了中断续算的，里面的step定义是比较模糊的，特别是一次模拟中进行了多个过程时，建议用d ...

您好，请问一下假如第一次输出轨迹的时候没有计算MSD、RDF等一些参数，那么在第二次在不输出轨迹文件以及不进行rerun的情况下专门去计算这些参数，这样的结果正确吗？

Graphite

zhaoxiaollong 发表于 2023-10-19 22:48
您好，请问一下假如第一次输出轨迹的时候没有计算MSD、RDF等一些参数，那么在第二次在不输出轨迹文件以及 ...

可以，和保存轨迹再算是一样的。

wilesyu

大佬，您的data文件是咋建立的？