msd计算方法
compute ID group-ID msd keyword values ...
计算一个四元的向量，是global vector,前三项分别是
$dx^2$… 最后一项是加和ID 可以是数字也可以是字符串
group-ID 组的名字，需要提前定义一个组
msd compute的style
keyword =com or average
keyword = com values = yes or no
默认值为com =no average=no
keword = average values = yes or no
com yes 表示 计算当前选定的组的原子的质心，并在组内去除质心漂移的效果
avergae yes 表示 计算msd的参考点设置为这个原子的平均位置 只适合没有原子扩散的时候，比如热传导
example 五洲 /public/home/liz/workspace/4_error_bar_for_na/1/800/test/test2/msd_test
lammps计算时默认是固定所有原子的整体的质心的模式来做的。
以Na5YSi4O12为例，分别计算all 和Na分组的msd ,com 分别设定为yes和no。预期结果是all 的yes和no 相同，因为模拟是固定了质心的，所以all 没有质心漂移。而na分组的yes和no的结果不同 ，因为na的质心没有固定.做这个测试是为了说明在分组计算msd的时候，不能用com=yes，因为会导致不正确，只去除了扩散原子的质心，但是扩散原子质心的变化并不是真正的质心漂移
正确的应该是以不动的原子为基准，确立一个质心，及其后续模拟中的质心漂移，用扩散的原子的msd减去质心的漂移得到真正的扩散原子的msd，如果在模拟晶胞时，原子扩散的取向是各向同性的，那Na分组的质心应该也不会有太大变化。
以下是结果，和预测相同，all没有变化，但是Na有，不能使用Na com=yes, 但是na com=no又没有去除质心的漂移
质心漂移、关于fix recenter
https://tinevez.github.io/msdanalyzer/tutorial/MSDTuto_drift.html
由于在进行分子动力学模拟时，系统的动量为零，因此质心的速度为零（质心速度公式与总动量相关），当一个离子向前扩散时，此时离子有向前的速度，那么系统整体的除了这个原子的总动量质心有向后的速度。这里有两种观看视角，一种是在质心的位置，观察到的是一个离子向前运动，其他的离子向后运动，一种是在除了这个扩散离子的其他的离子整体的位置看，离子有一个位移。计算相对质心位置的msd，那就是lammps现在做的事情，但是msd应该计算相对除了这个扩散离子的其他的整体的位移，这需要固定住除了这个扩散离子后的部分
这是模拟的限制而不是现实
见ppt L-2-2024 4_17.ppts
动量守恒与质心不变是互易的，因此，既可以按移动原子后移动质心回到原点的图像，也可以按在原子有动量时其他原子有相反的动量的图像来考虑物理的过程。这两种思考的统一结论是，msd的原子应该是相对于原来位置的移动距离而不是相对质心的移动距离。从模拟出的图像上看，是没有问题的，但是轨迹线是有问题的
如果能够保证扩散在各个方向上是各向同性的，那轨迹线是没问题的，但是如果是各向异性的，那就不合理。各向同性要求晶胞足够大。
一个问题是为什么会发生原子的‘聚集’扩散现象
如何计算正确的Na的msd ，把不动的那些原子的质心固定住，用他们作为一个基准，这样，计算出来的钠的msd就是正确的
fix recenter是与fix nve/nvt/ 这些命令同时进行的，应该放在这些会改变原子信息的命令之后，fix recenter的执行顺序与它的位置无关，即便在之前有计算质心的命令，输出的质心还是改变之后的。
fix recenter 中shift 的使用
不加shift就只移动recenter的原子群
加了shift all 才移动所有的原子
不加shift all 是不对的，会破坏原子间相互位置关系
fix 1 ysio recenter INIT INIT INIT shift all
https://matsci.org/t/fix-recenter-question/15131/6
在计算中发现的不动原子的msd的增加是正常的，是上述质心不动的原因，发现的不动离子的移动也是可以理解的了。
对于计算msd，问题在于以什么为基点，是计算离子相对质心的运动还是相对初始位置的运动（初始位置由与其他原子的位置关系确定）。在其他原子组成的质心没有发生太大变化时，这两者是一致的
现在，对于msd的理解，直接使用相对质心的位移，并且不能计算所有的原子的总均方位移。使用相对真实位置的位移是真实的物理环境，因为在kysio中，使用相对质心的位移，那么从结果看，ysio也在运动，这是不对的，但是在模拟时无法做到计算相对真实物理世界的位置的。
考虑液体时，某一类离子向不同方向的扩散是各向同性的，那么这类离子不会带来其他离子的位置的改变，进而这类离子的扩散是真实的。
非常重要
一旦在模拟中出现某个方向上各向异性的扩散或者移动，那么此时的模拟图像是准确的，但是msd计算就是不准确的了，这也是分子动力学模拟大胞的必要性，可以消除这种各向异性
在某些特定时间段内，离子向一个方向的聚集扩散没有理由认为是错误的（极简时可以考虑一个体系中只有一个足够重的间隙离子在运动）
如何判断模拟的尺寸是否足够大？看计算的不同种类的离子的质心是否保持稳定。在计算分子动力学过程中，理论上要求所计算的不同种类的原子的质心都保持稳定
小胞需要模拟足够长的时间才能消除这种各向异性，在某些特定时间段内，离子向一个方向的聚集扩散没有理由认为是错误的（极简时可以考虑一个体系中只有一个足够重的间隙离子在运动），如果为了计算msd，这种局部的各向异性需要计算足够长的时间来消除。在足够长时间的模拟消除各向异性后，计算msd就可以使用相对质心的计算来代表了，各向异性带来波动
如何判断模拟的时间是否足够长，如何判断模拟的尺寸是否足够大？看计算的不同种类的离子的质心是否保持稳定。在计算分子动力学过程中，理论上要求所计算的不同种类的原子的质心都保持稳定，模拟大的晶胞是更划算的
以质心为标定点，对于一个离子扩散出晶胞，在模拟中这个离子的msd的倾向是小于真实的值，但是对于一个离子向晶胞扩散，在模拟中这个离子的msd的倾向是大于向内扩散的距离
离子向外扩散时，真实是B，按质心计算是A，但是离子向内扩散时，质心计算是A，实际是B。不论向内还是向外，各向异性的扩散在使用质心作为标定点时，计算的msd总是偏小的，这种误差也许是可以被忽略的，因为以质心看，除了扩散离子的晶体的其他部分的位移（两部分，一部分扩散原子，一部分非扩散原子，其中非扩散原子多且重于扩散原子的情形）应该是比较小的，因为小部分扩散原子带不动大部分原子。
模拟和真实是完全不同的两种情况
当需要处理两帧之间的nvt系综确定的量时，这种问题就会出现
lammps的NPT模拟输出的unwraped轨迹不能用于计算扩散系数！！，因为随着距离中心原子的距离增加，其轨迹是失真的（见Unwrapping NPT Simulations to Calculate Diffusion Coefficients）
应该尽量使用NVT而不是NPT
计算时固定质心和折叠轨迹
固定质心后输出折叠后的轨迹（x y z)，计算质心会发现质心变化，但这种变化是折叠轨迹导致的，而不是质心真的发生了变化，计算仍然是固定了质心，需要反折叠(xu yu zu)才能看到真正的质心
xu yu zu和x y z的区别
xu yu zu是未折叠的轨迹，
group stay type 2 3 4 fix 7 stay recenter INIT INIT INIT
能够固定不动原子的质心，把跑出去的轨迹拉到使得不动原子保持质心不动的位置，优先级高于默认的使得整个体系质心保持不动，变相地做到了只考虑运动原子的扩散
4. 计算MSD
msd图 ，msd斜率收敛图 温度压力能量图三者不可缺
1.对比lammps自带的msd，自己脚本的，和vmd的
所用机器：华为/home/jildxwlxyljlstdui/cccs-share02/lijx/liz/2-Na-X-Si-O/2_NaXSiO_big_model/2_Sb_test/4_test_formsd
自带
脚本
VMD
vmd和自己脚本算的基本一致，lammps增加了不动的部分原子，最后计算出来也和自己脚本相差无几 （VMD算出来的额数字需要先平方才能比较）
YSIO不动，导致lammps算出的总msd有一个上移的步骤，不论如何，不应该考虑这部分向上平移。最开始算不出直接的平移是因为单种原子本身的均方根位移也为0，还没有到平衡位置，间隔一段时间后才能显示出平移的效应
在去掉平移的一部分后，两个获得的斜率相差很小，因此，可以用lammps自带的msd，但是注意这时候的N为所有原子总数包括Na和非钠的所有元素
关于TMSD和MSD和D
电导率最终看的是TMSD，只需要保证TMSD对时间的斜率最终相同就不会错
需要注意的是，扩散系数的值是和包括的不动的原子的数目相关的，包括的越多，扩散系数越小，因此，这里的扩散系数并不是真正的扩散系数，而是修改版的（因为除以的原子数中有一部分没有动）。真正的扩散系数应该只看运动的那些
对任何一个系统，不管算哪些原子，扩散系数肯定是不变的，变化的只是原子数，但是只有这些原子贡献了斜率，所以只需要求出所有的TMSD，求出斜率后换算到对应的原子数就得到了准确的扩散系数
严格来说，均方根位移应该是对所有移动的原子做平均均方，所以我现在做的不算是均方根位移，因为包括了那些不移动的原子
2.python脚本
三处需要修改，1处可能需要改
import numpy as np
import time
def timer(func):
def wrapper(*args, **kwargs):
start_time = time.time()
result = func(*args, **kwargs)
end_time = time.time()
print(f"Time elapsed: {end_time - start_time:.4f} seconds")
return result
return wrapper
@timer
def read_lmp_trj(file_name: str):
f = open(file_name, 'r')
"""data segment"""
n_atoms_list = []
lat = []
pos = []
num1 = 0
#跳过文件的行数，需要修改
num_lines_to_skip = (132000+9)*100
current_line_num = 0
"""temp"""
# n_atoms = None
while True:
line = f.readline()
if not line:
break
current_line_num += 1
if current_line_num <= num_lines_to_skip:
continue
if "NUMBER OF ATOMS" in line:
line = f.readline()
n_atoms = int(line)
n_atoms_list.append(n_atoms)
f.readline()
"""lat"""
for i in range(3):
lines = f.readline()
lat.append(lines)
f.readline()
"""pos"""
for i in range(n_atoms):
pos.append(f.readline())
#每隔1000个时间间隔读取一次 需要改
hhh = (n_atoms+9)*999
for m in range(hhh):
lines2 = f.readline()
num1 +=1
print(num1)
"""post treatment"""
#固定一个三维数组的两个维度，第三个维度也固定了，numpy的逻辑应该是3-2-1，第三个维度先排布，依次向前，第一个是深度
lat = np.array([i.split() for i in lat]).astype(float).reshape(-1, 3, 3)
#需要看这里的8需不需要改
pos = np.array([i.split() for i in pos]).astype(float).reshape(len(lat), -1, 8)
# idx = pos[:, :, 0].argsort(1)
#lexsort的顺序是从最后一个传入的参数开始，依次到第一个传入的参数
idx =np.lexsort((pos[:,:,0],pos[:,:,1]))
sort_pos = pos[np.arange(len(pos))[:, None], idx, :]
start3 = sort_pos[:,:,0:2].reshape(-1,2)
n_atoms_list = np.array(n_atoms_list)
return n_atoms_list, lat, sort_pos[:, :, [0, 2, 3, 4]],start3
if __name__ == '__main__':
file_name = "A.lammpstrj"
print("Loading lmp trajectory")
n_atoms, _, x,stc = read_lmp_trj(file_name=file_name)
x = x[:, :, 1:]
#np.savetxt('first_element.txt', x[0])
#np.savetxt('end_element.txt',x[-1])
print("X shape : ", x.shape)
print("N frames: ", n_atoms.shape)
#30这里需要修改，是移动的原子个数
res = np.square(x - x[0])[1:,:30,:].sum((1, 2))
res = res / 30
# print(res)
np.save("./res2", res)
arr = np.load('res2.npy')
#print(arr)
rows = np.arange(arr.shape[0])
rows += 1
#print(rows)
stac = np.stack([rows, arr], axis=1)
#print(stac)
np.savetxt('newdata.txt',stac,fmt='%.5f',delimiter=',')
np.savetxt('stc.txt',stc,fmt='%d',delimiter=',')
转载请注明来源 有问题可通过github提交issue