rdf, msd, group/group, com 计算脚本

看到有大神发关于这几个命令的计算代码，自己收藏！

1 rdf

rdf又称g( r )，rdf的计算与输出通过以下命令计算输出：

 

　　1 compute ID group-ID rdf Nbin itype1 jtype1 itype2 jtype2 ...　　

　　2 fix ID group-ID ave/time Nevery Nrepeat Nfreq value1 value2 ... keyword args ... 

 

　　其中，Nbin类似于切片数，数也多，可能曲线会更连续，光滑，但过多会影响计算效率。一般200-500足够。itype1是径向分布函数的中心原子， jtype1是被统计的原子。fix ave/time命令能够对任何全局量进行时间平均，

由用户指定一个或多个需要进行时间平均的物理量。

　　Nevery，Nrepeat和Nfreq参数指定在什么时间步上使用输入值以贡献平均值。 最终的平均每Nfreq生成产生一次。 由相隔Nevery个时间不到的Nrepeat个值得到、因此，Nfreq必须是Nevery的整数倍，Nevery必须是非零

（即使Nrepeat为1，Nevery也必须为非零）。且用于计算的相邻平均值的时间步长不能重叠，即Nrepeat *Nevery每个值都不能超过Nfreq。例如，如果Nevery = 2，Nrepeat = 6和Nfreq = 100，则时间步90,92,94,96,98,100的

值将为用于计算时间步长100的最终平均值。类似地，时间步长190,192,194, 196,198,200的值将用于计算时间步长200的值。如果Nrepeat = 1且Nfreq = 100，则不进行时间平均，命令产生的值是时间步100，200的值。

value=c_ID, c_ID[N], f_ID,f_ID[N], v_name。其中，如果是通过compute命令计算的结果，则用“c_ID”或者添加后缀输出。如果是通过fix命令计算的结果，则通过"f_ID"或者添加后缀输出。这里rdf是通过compute计算得到，

后面输出则为c_ID。mode可选scalar或者vector。这里rdf输出为全局矢量，故选mode vector。输出结果按列依次为序号，c_ID[1]，c_ID[2]，c_ID[3]。c_ID[1]为group到原点的距离，即x轴，c_ID[2]为g( r )，c_ID[3]为coor( r )，

即 g( r )的积分。例如以下命令可用于输出每500000步的rdf值，在截断半径范围内，距离被分成100段：

 

　　1 compute rd1 hbond1 rdf 100 2 6　　#rdf of OH-Ow

　　2 fix hb1 hbond1 ave/time 500000 1 500000 c_rd1[1] c_rd1[2] c_rd1[3] file rdf.hb1 mode vector

 

2 msd

　　msd计算材料的均方位移，msd计算与输出通过以下命令计算输出：

1 compute ID group-ID msd keyword values ... 
2 fix ID group-ID ave/time Nevery Nrepeat Nfreq value1 value2 ... keyword args ...  

　　其中keyword = com或者average，值可以为yes或者no。

　　输出结果按列依次为c_ID[1]，c_ID[2]，c_ID[3]，c_ID[4]。前3行一次为xyz方向上的msd值，第4行为总的msd值。例如以下命令可用于输出每500000步的msd值：

　　1 compute md1 w_03A msd
　　2 fix msd01 w_03A ave/time 500000 1 500000 c_md1[*] file w_03A.msd mode vector

3 group/group

　　group/group计算材料的相互作用能，由以下命令计算输出：

　　1 compute ID group1-ID group/group group2-ID keyword value ... 
　　2 thermo_style style args 
　　3 thermo Nsteps

　　其中，关键词 pair value = yes (or no)；kspace value = no(or yes)；boundary value = yes(or no)，通常按默认的来即可。thermo_style按照自己需要的格式输出，具体可以见手册。

　　相互作用能为全局标量，可通过thermo_style输出，输出结果用c_ID。例如以下命令可用于输出每5000步的msd值：

　　1 compute tot sio2_l group/group sio2_r #interaction
　　2 thermo_style custom step temp etotal c_tot
　　3 thermo 5000 

4 com

　　com用来计算材料的质心，由一下命令输出：

　　1 compute ID group-ID com
　　2 thermo_style style args 
　　3 thermo Nsteps

　　相互作用能为全局矢量，可通过thermo_style输出，输出结果用c_ID[1]，c_ID[2]，c_ID[3]，分别代表质心的三个坐标。默认输出未折叠的坐标，如果粒子出了盒子，按照实际坐标输出，而非镜像粒子的坐标。

　　例如以下命令可用于输出每5000步的msd值：

　　1 compute centl sio2_l com
　　2 thermo_style custom c_centl[1] c_centl[2] c_centl[3] 
　　3 thermo 5000