在此提供一个计算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
    

复制代码

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

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

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

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