计算化学公社

标题: 可计算The residence time correlation function的程序或脚本 [打印本页]
作者
Author: Lacrimosa    时间: 2020-10-15 15:49
标题: 可计算The residence time correlation function的程序或脚本
请教一下各位老师们，有什么程序/轨迹分析插件/脚本可以统计MD轨迹中某个组分的The residence time correlation function呢
作者
Author: 一条君    时间: 2020-12-17 01:41
请问楼主解决了吗，同样想得到A-B的residence time
作者
Author: Lacrimosa    时间: 2021-3-14 23:20
一条君 发表于 2020-12-17 01:41
请问楼主解决了吗，同样想得到A-B的residence time

没有解决呢, 我尝试过根据公式自己写, 但是一直没成功. 应该是我对公式理解有问题
作者
Author: Lacrimosa    时间: 2021-3-14 23:38


比如要计算离子周围水的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的区间波动

在此请各位老师指点一下
作者
Author: southeastBig    时间: 2021-9-3 17:11
楼主解决了吗，同想求离子周围水分子的residence time
作者
Author: Lacrimosa    时间: 2022-12-10 15:21
在此提供一个计算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


  11. call getarg(1,fname)
  12. open(10,file=fname,status='OLD',form='FORMATTED',action='READ')


  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


  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
复制代码




欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3