NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算

fhh2626

秋心 发表于 2023-5-9 11:11
付老师，您好，我最近在学习lammps+colvars使用US计算一个有机分子从气相到液相的自由能变化，我已经实现了 ...

harmonic {  
  colvars        dist
  outputEnergy   yes
  forceConstant  50.0
  centers        1.95915304014070e+01
}

秋心

fhh2626 发表于 2023-5-9 12:57
harmonic {  
  colvars        dist
  outputEnergy   yes

收到，非常感谢老师您的解答

YZH_jesse

付老师您好：


我在尝试LAMMPS+colvars计算钾离子（绿色小球）通过石墨烯孔过程中的迁移势垒（当前体系已达到平衡态）：我选定了一个位于孔中心正上方的钾离子作为main group，选定了孔周边一圈原子作为ref group，以下是我使用的colvars脚本

1. # collective variable for PMF
   
2. indexFile Graphene_pore.ndx
   
3. 
   
4. colvar {
   
5.   name dist
   
6.   lowerBoundary   79.144   #钾离子关于石墨烯膜对称的目标位置z轴坐标
   
7.   upperBoundary  135.144  #钾离子z轴坐标
   
8.   extendedlagrangian     on
   
9.   extendedFluctuation    1.0
   
10.   extendedTimeConstant   200
    
11.   extendedTemp           300.0
    
12.   distanceZ {
    
13.     main {
    
14.       atomNumbers 56528
    
15.     }
    
16.     ref {
    
17.       indexGroup G_pore
    
18.     }
    
19.         axis (0.0,0.0,1.0)
    
20.   }
    
21. }
    
22. 
    
23. harmonicWalls {
    
24.   name walls
    
25.   colvars dist
    
26.   lowerWalls 33.4124 #下石墨烯挡板坐标
    
27.   upperWalls 181.359 #上石墨烯挡板坐标
    
28.   forceConstant 1.0
    
29. }
    
30. 
    
31. abf {
    
32. colvars dist
    
33. fullSamples 500
    
34. historyfreq  1000000
    
35. writeCZARwindowFile #此处参考了付老师在此贴首页的教程
    
36. }

复制代码

目前遇到的问题是：在fix colvars计算过程中钾离子会卡在中间石墨烯膜的非孔位置，无法迁移到另一侧的溶液。计算时间尺度为4ns，时间步长2fs。
想请教一下付老师：
1. colvars脚本的编写是否有问题
2. 如何使钾离子在fix colvars计算中迁移到溶液另一侧

非常感谢老师！

秋心

付老师，您好，lammps+colvars使用US方法计算分子从气相到液相的自由能变化时，在每个窗口的采样模拟中按照下面的方法进行设置（以其中一个窗口为例），得到的us.colvars.traj中距离是负值，很明显不在我设置的墙里面，对此感到疑惑。希望得到老师您的指导。非常感谢。
colvarsTrajFrequency        100
colvarsRestartFrequency     100

colvar {     #设置研究的cv
  name   dist
  width  0.0005
  lowerBoundary 1.80871854278497e+01
  upperBoundary 2.25915304014070e+01
  lowerWallConstant 30.0
  upperWallConstant 30.0
  distanceZ {
    main { atomNumbersRange 2401-2410 }
    ref { dummyAtom (11.5, 11.5, 79.5) }
  }
}

harmonic {  
  colvars        dist
  outputEnergy   yes
  forceConstant  1.0
  centers        2.0e+01
}

histogram {  #用来记录设置cv 的贡献，用于后面的分析
  colvars  dist
}

fhh2626

YZH_jesse 发表于 2023-5-9 22:46
付老师您好：

添加一个distanceXY的约束（用harmonicWall），将离子在xy平面的位置约束在一个圆柱体内

fhh2626

秋心 发表于 2023-5-10 09:41
付老师，您好，lammps+colvars使用US方法计算分子从气相到液相的自由能变化时，在每个窗口的采样模拟中按照 ...

这个看不出来，你自己在检查一下吧

YZH_jesse

fhh2626 发表于 2023-5-10 11:19
添加一个distanceXY的约束（用harmonicWall），将离子在xy平面的位置约束在一个圆柱体内

收到！非常感谢老师解答！

秋心

fhh2626 发表于 2023-5-10 11:21
这个看不出来，你自己在检查一下吧

好的，谢谢付老师您的帮助。另外，还有一个问题我想问一下，我看到您之前的一个帖子http://bbs.keinsci.com/forum.php ... 12240&fromuid=30893里面提到了限制水板质心，是如何实现限制水板质心的？我翻看了Colvars手册并没有找到相应的命令，于是我在in文件中使用了fix 2  spce spring/self 100.0 z命令，对水板中的所有水分子在Z方向上添加了大小为100 kcal/mol的弹簧力，限制住水板在z方向上的移动，再去做伞形采样模拟，这样做可以吗？我这样做得到图一的结果，感觉不太正常，但是又不知道哪里出错了。感谢老师您的指导！
colvarsTrajFrequency        100
colvarsRestartFrequency     100

colvar {     #设置研究的cv
  name   dist
  width  0.0005
  lowerBoundary 1.80871854278497e+01
  upperBoundary 2.25915304014070e+01
  lowerWallConstant 30.0
  upperWallConstant 30.0
  distanceZ {
    main { atomNumbersRange 2401-2410 }
    ref { atomNumbersRange 1-2400 }
    forceNoPBC yes
  }
}

harmonic {  
  colvars        dist
  outputEnergy   yes
  forceConstant  2
  centers        1.95e+01
}

histogram {  #用来记录设置cv 的贡献，用于后面的分析
  colvars  dist
}

YZH_jesse

fhh2626 发表于 2023-5-10 11:19
添加一个distanceXY的约束（用harmonicWall），将离子在xy平面的位置约束在一个圆柱体内

付老师您好：

再次感谢您的回答！

以下是我提交计算后3ns内的结果：

初始状态（已平衡）钾离子在孔的z轴上方，在很短时间内钾离子就跑到了孔的下方，在Z轴位置就几乎保持不变了；从采样上来看，孔上方空间的采样几乎为零。从PMF的数据来看Z大于70的自由能随时间几乎没有变动。

付老师，我想问一下：
1. 当前现状会随计算的进行有所改善吗（我在继续跑后面的计算中）？
2. 为什么初始钾离子的坐标变化会变化这么快（我的width选取为1，是否偏大）？

PS：以下是我计算过程中使用的colvars脚本：

1. # collective variable for PMF
   
2. colvarsTrajFrequency    25000
   
3. colvarsRestartFrequency 1000000
   
4. 
   
5. colvar {
   
6.   name dist
   
7.   lowerBoundary   45.744
   
8.   upperBoundary  101.744
   
9.   extendedlagrangian     on
   
10.   extendedFluctuation    1.0
    
11.   extendedTimeConstant   200
    
12.   extendedTemp           300.0
    
13.   distanceZ {
    
14.     main {
    
15.       atomNumbers 56528
    
16.     }
    
17.     ref {
    
18.       atomNumbers {35 36 39 40 44 46}
    
19.     }
    
20.         axis (0.0,0.0,1.0)
    
21.         forceNoPBC yes
    
22.   }
    
23. }
    
24. 
    
25. colvar {
    
26.   name distxy
    
27.   upperBoundary  4.5
    
28.   distanceXY {
    
29.     main {
    
30.       atomNumbers 56528
    
31.     }
    
32.     ref {
    
33.       atomNumbers {35 36 39 40 44 46}
    
34.     }
    
35.         axis (0.0,0.0,1.0)
    
36.   }
    
37. }
    
38. 
    
39. harmonicWalls {
    
40.   name walls
    
41.   colvars dist
    
42.   upperWalls 101.744
    
43.   lowerWalls  45.744
    
44.   forceConstant 1.0
    
45. }
    
46. 
    
47. harmonicWalls {
    
48.   name wallsxy
    
49.   colvars distxy
    
50.   upperWalls 4.5
    
51.   forceConstant 1.0
    
52. }
    
53. 
    
54. abf {
    
55. colvars dist
    
56. fullSamples 500
    
57. historyfreq  1000000
    
58. writeCZARwindowFile
    
59. }

复制代码

fhh2626

YZH_jesse 发表于 2023-5-11 12:08
付老师您好：

再次感谢您的回答！

你dist的范围是正数，不就是代表只在一侧采样吗？

fhh2626

秋心 发表于 2023-5-10 18:48
好的，谢谢付老师您的帮助。另外，还有一个问题我想问一下，我看到您之前的一个帖子http://bbs.keinsci.c ...

在Colvars限制一个原子集合的质心也是用harmonic约束distance，其实US本质就是在限制质心的情况下进行模拟

秋心

fhh2626 发表于 2023-5-11 13:44
在Colvars限制一个原子集合的质心也是用harmonic约束distance，其实US本质就是在限制质心的情况下进行模 ...

收到，感谢老师您的再次解答！

YZH_jesse

fhh2626 发表于 2023-5-11 13:42
你dist的范围是正数，不就是代表只在一侧采样吗？

在现在的计算里面我是将最下层的石墨烯挡板作为了ref group，这个dist范围对应于石墨烯孔上方和下方各28A的范围（总长56A），但是钾离子只倾向于在孔下方28A范围内运动，孔上方空间的几乎没有采样。

感谢老师指正！

fhh2626

YZH_jesse 发表于 2023-5-11 19:39
在现在的计算里面我是将最下层的石墨烯挡板作为了ref group，这个dist范围对应于石墨烯孔上方和下方各28A ...

有可能是时间太短，再等等

另外一个可能是你这个定义方式（距离超过盒子长度的一半）会导致Colvars内部自动进行了wrap，试试ForceNoPBC

秋心

付老师，再次打扰您一下，按照文献中的描述（文献中使用的程序也是lammps+colvars），使用1kcal/mol 的弹簧力对一个有机分子（我研究的也是有机分子）在液相中进行伞形采样，每个窗口的采样区间应该是在1.5 埃左右，而我按照下面的设置使用lammps+colvars去跑伞形采样的时候每个窗口的采样区间只能达到0.5 埃左右（图一和图二是其中一个窗口的采样结果），并且发现300 ps和1 ns的采样区间区别很小(图一和图二)，为了解决这个问题我尝试过修改弹簧力大小和增加模拟时长至5 ns，弹簧力大小为0.5 kcal/mol 和2 kcal/mol，模拟时长均为5 ns时，会使采样区间变大，但是采样区间最多能达到0.9 埃左右。我觉着问题大概率不是弹簧力大小和模拟时长的问题，应该是我COLVARS中设置或者lammps中设置存在的一些问题，我也不太确定对不对，麻烦老师能否帮我检查一下COLVARS的设置和lammps的设置（附件中的in.complex）是否有问题。谢谢老师！
colvarsTrajFrequency        100
colvarsRestartFrequency     100

colvar {     #设置研究的cv
  name   dist
  width  0.1
  lowerBoundary 18.5
  upperBoundary 21.5
  lowerWallConstant 30.0
  upperWallConstant 30.0
  distanceZ {
    main { atomNumbersRange 2401-2410 }
    ref { atomNumbersRange 1-2400 }
    forceNoPBC yes
  }
}

harmonic {  
  colvars        dist
  outputEnergy   yes
  forceConstant  1
  centers        19.5
}

histogram {  #用来记录设置cv 的贡献，用于后面的分析
  colvars  dist
}