可计算The residence time correlation function的程序或脚本

Lacrimosa

请教一下各位老师们，有什么程序/轨迹分析插件/脚本可以统计MD轨迹中某个组分的The residence time correlation function呢

Lacrimosa

在此提供一个计算tip5p水分子转动自相关函数的计算代码，通过修改可以实现对residence time correlation function的计算，个人编程水平低下，欢迎指教。将以下内容保存为racf.f90, linux端执行gfortran racf.f90 -o racf.x即可生成可执行文件racf.x. 用法：./racf.x tip5p.gro >racf.dat
以下为程序代码:

1. program main
   
2. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
   
3. character fname*216, dum*216
   
4. real(16) :: mu(2000,10000,3), O(3), H1(3), H2(3) !vector mu(nframe,nmolecule,xyz)
   
5. real(16) :: dohh, ddot
   
6. integer :: fr
   
7. integer :: dt
   
8. integer :: atom
   
9. 
   
10. 
    
11. call getarg(1,fname)
    
12. open(10,file=fname,status='OLD',form='FORMATTED',action='READ')
    
13. 
    
14. 
    
15. dohh=0.058588
    
16. nl=0
    
17. do
    
18.    read(10,*,end=900)
    
19.    nl=nl+1
    
20.    read(10,*) num
    
21.    num = num/5
    
22.    do i=1,num
    
23.      read(10,"(A20,3F8.3)") dum, O(1), O(2), O(3)
    
24.      read(10,"(A20,3F8.3)") dum, H1(1),H1(2),H1(3)
    
25.      read(10,"(A20,3F8.3)") dum, H2(1),H2(2),H2(3)
    
26.      read(10,*) dum
    
27.      read(10,*) dum
    
28.      mu(nl,i,:) = (O(:)-(H1(:)+H2(:))/2)/dohh
    
29.    enddo
    
30.    read(10,*) bx,by,bz
    
31. enddo
    
32. 900  continue
    
33. 
    
34. 
    
35. n_sample=INT(nl/2-1)
    
36. do dt=1,n_sample
    
37.    d_sum=0
    
38.    do fr=1,nl-dt
    
39.      do atom=1, num
    
40.        ddot=mu(fr,atom,1)*mu(fr+dt,atom,1)+mu(fr,atom,2)*mu(fr+dt,atom,2)+mu(fr,atom,3)*mu(fr+dt,atom,3)
    
41. !   write(6,*) ddot
    
42.        d_sum=d_sum+ddot
    
43.      enddo
    
44.    enddo
    
45.    d_sum=d_sum/(nl-dt)/num
    
46.    write(6,"(I8,F16.3)") dt, d_sum
    
47. enddo
    
48. return
    
49. end
    

复制代码

southeastBig

楼主解决了吗，同想求离子周围水分子的residence time

Lacrimosa

比如要计算离子周围水的residence time correlation function ,在vmd中 先在0ps选取离子周围所有的水分子, 记录下resid和水分子数量N, 这个时候R(0) = 1
然后过1ps再选取所有离子周围的水分子,对比resid, 如果resid与初始的相同就取1 ,不同就取0 这样把所有的水分子考察一遍之后求和除以N 结果即为R(1)
以此类推

我按照以上思路写完脚本后遇到两个问题:
1. 每次选取的离子周围的水分子数量可能略有差别, 那么N的值怎么选取, 是按照最初的还是当前的呢? 无论按照哪个,当水分子数不同的时候相差的几个都会算做是0, 这样还准确么?
2.用一段2 ns的轨迹运行之后发现离子周围的水分子基本都是最初的, 没什么变化, 所以求出来的R(t)就在一个极其接近1的区间波动

在此请各位老师指点一下

Lacrimosa

一条君 发表于 2020-12-17 01:41
请问楼主解决了吗，同样想得到A-B的residence time

没有解决呢, 我尝试过根据公式自己写, 但是一直没成功. 应该是我对公式理解有问题

一条君

请问楼主解决了吗，同样想得到A-B的residence time