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

fhh2626

海棠依旧cdut 发表于 2025-1-26 11:24
付博士您好，我想问一下您，用colvars进行伞形采样时，谐波力的计算方法是forceconstant/width^2吗？也就是 ...

是的
用distanceXY作为CV施加约束。但是修改axis选择YZ平面

海棠依旧cdut

付博士您好，不好意思再次打扰您。我在用wham计算伞形采样PMF时候出现了free=inf的问题。以下图是我施加K=100kcal后拉伸的离子的概率密度，范围是45.2-48.0埃（膜处于48埃，我模拟的是离子从体溶液中45埃拉伸到膜孔的PMF），间距为0.2埃，一共15个窗口，一个窗口运行2ns。从图中看到，各处的概率密度都有较好的重叠，这代表采样已经充分。然后我将lammps输出的文件放入wham分析，分析结果出现inf，据我所知，出现inf的原因是因为该处没有采样导致lg0=inf，可是从概率密度图中可以看到实际上每个位置都是采样了的，请问是什么原因造成的inf呢？

https://lammpstutorials.github.io/sphinx/build/html/tutorials/level3/free-energy-calculation.html。这个是我的代码参考网站

wham命令：./wham 45.2 48 15 1e-10 300 0 metadata.dat PMF.dat

fhh2626

海棠依旧cdut 发表于 2025-2-3 22:40
付博士您好，不好意思再次打扰您。我在用wham计算伞形采样PMF时候出现了free=inf的问题。以下图是我施加K=1 ...

这是wham代码的问题吧，我不清楚

qmlearner

海棠依旧cdut 发表于 2025-1-26 11:24
付博士您好，我想问一下您，用colvars进行伞形采样时，谐波力的计算方法是forceconstant/width^2吗？也就是 ...

您好，我们也在算这种体系的PMF，发现即使离子能过孔，但得到的轨迹中没有离子处于孔中的构型，所以不确定这是否意味着采样不充分，我们用的是ABF方法，您用伞状采样方法您的体系有没有这种情况呢？期待交流

xukw

po390 发表于 2024-9-27 18:08
付老师，您好
这个meta-eABF直接应用于NAMD的QM/MM中是否可行，例如我有涉及化学反应的多个原子直接的五个 ...

我用2.14版本的NAMD调用ORCA，然后ORCA调用xtb，可以min/heat，但是加入colvars约束之后，无法继续进行？把QMforce关闭之后，又可以正常计算。请问这个是版本之间不兼容吗？
colvars: Updating NAMD interface:
colvars: updating atomic data (0 atoms).
colvars: updating group data (6 scalable groups, 6 atoms in total).
colvars: updating grid object data (0 grid objects in total).
colvars: Re-initialized atom group for variable "rc1":0/0. 1 atoms: total mass = 1.008, total charge = 0.
colvars: Re-initialized atom group for variable "rc1":0/1. 1 atoms: total mass = 16, total charge = 0.
colvars: Re-initialized atom group for variable "rc1":1/0. 1 atoms: total mass = 1.008, total charge = 0.
colvars: Re-initialized atom group for variable "rc1":1/1. 1 atoms: total mass = 12.01, total charge = 0.
colvars: Re-initialized atom group for variable "rc1":2/0. 1 atoms: total mass = 12.01, total charge = 0.
colvars: Re-initialized atom group for variable "rc1":2/1. 1 atoms: total mass = 12.01, total charge = 0.
colvars: Re-initialized atom group for variable "rc2":0/0. 1 atoms: total mass = 12.01, total charge = 0.
colvars: Re-initialized atom group for variable "rc2":0/1. 1 atoms: total mass = 12.01, total charge = 0.
colvars: The restart output state file will be "meta.restart.colvars.state".
colvars: The final output state file will be "meta.colvars.state".
colvars:   Prepared sample and gradient buffers at step 0.
colvars: Synchronizing (emptying the buffer of) trajectory file "meta.colvars.traj".
Info: List of ranks running QM simulations: 0.
QMENERGY:       0         1.0000    -66787.7034    -65698.6764

PRESSURE: 0 -4721.45 -32.0357 -88.0976 -41.1786 -4885.46 -312.823 -89.4787 -305.058 -4273.71
GPRESSURE: 0 -4676.27 -16.1009 -97.309 -49.8086 -4839.46 -320.3 -103.88 -305.36 -4192.85
ETITLE:      TS           BOND          ANGLE          DIHED          IMPRP               ELECT            VDW       BOUNDARY           MISC        KINETIC               TOTAL           TEMP      POTENTIAL         TOTAL3        TEMPAVG            PRESSURE      GPRESSURE         VOLUME       PRESSAVG      GPRESSAVG

QMETITLE:      TS           QMID         ENERGY  PMECORRENERGY

ENERGY:       0       341.5874      1274.5170      5463.4951         0.0000        -232807.6619     25449.7334         0.0000    -65698.6764      2294.3551        -263682.6502        20.0582   -265977.0053   -263682.1328        20.0582          -4626.8727     -4569.5245    604504.0000     -4626.8727     -4569.5245

OPENING EXTENDED SYSTEM TRAJECTORY FILE
QMENERGY:       1         1.0000    -66787.3567

Segmentation fault (core dumped)

xukw

xukw 发表于 2025-2-17 22:51
我用2.14版本的NAMD调用ORCA，然后ORCA调用xtb，可以min/heat，但是加入colvars约束之后，无法继续进行？ ...

换成了AMD的CPU，然后报这个错OwnerBox::close() - close called, but closeCount 0 openCount -1

海棠依旧cdut

qmlearner 发表于 2025-2-17 11:37
您好，我们也在算这种体系的PMF，发现即使离子能过孔，但得到的轨迹中没有离子处于孔中的构型，所以不确 ...

好像就是因为ABF没有办法获得某个不利位置的PMF，我才选用伞型采样的(小白理解，如果有更好的解释请大佬指点)。实际上伞型采样在膜孔处的采样结果也是不太理想的(从输出的直方图来看)，参考这篇文章《Molecular Dynamics Study of Mg2+/Li+ Separation via Biomimetic Graphene-based Nanopores: The Role of Dehydration in Second Shell 》，但是学术上好像普遍认为histogram有交叉就是采样均匀，。所以对于我这种二维膜体系，我选择的方法是只测膜一侧的PMF，也就是-5到0这么一个区间。用x=0处的PMF来代表离子处于膜孔的能量

JCenter

付老师，我使用metaDynamics做聚合物中酯键的碱性氢氧根水解的自由能计算(CP2K中FPMD增强采样过程),想请教您2个关于化学反应中集合变量设置的问题：
1.集合变量只考虑断键和成键的原子对间的距离是不是不足以描述这个化学过程，是不是还要加入氢氧根离子周围水分子的配位数、成键的角度这两个呢，我能想到的就这些了，聚合物的构型这个变量这个我想给约束掉。
2.关于这个化学反应的过程，是不是要对聚合物进行位置约束，才能消除聚合物的运动这个变量对自由能的影响呢。希望付老师百忙之中有时间给学生解疑答惑一番，范感觉化学反应的集合变量考虑还蛮复杂的，不知道自己的CVs是否合理

fhh2626

JCenter 发表于 2025-5-6 00:18
付老师，我使用metaDynamics做聚合物中酯键的碱性氢氧根水解的自由能计算(CP2K中FPMD增强采样过程),想请教 ...

选反应坐标需要和体系有很大的关系，我也说不好

如果你只是距离、角度这样的变量的话，平动是没有影响的

JCenter

fhh2626 发表于 2025-5-6 09:42
选反应坐标需要和体系有很大的关系，我也说不好

如果你只是距离、角度这样的变量的话，平动是没有影响 ...

感谢付老师

JCenter

fhh2626 发表于 2025-5-6 09:42
选反应坐标需要和体系有很大的关系，我也说不好

如果你只是距离、角度这样的变量的话，平动是没有影响 ...

付老师，我这几天在细品您回复的内容，您应该是指考虑了成键、断键、两个原子间的角度等CVs后，不需要考虑约束聚合物的运动了吧，但是我一直没想明白，聚合物的运动也会造成自由能的变化吧，这样不约束聚合物运动，自由能景观图中岂不是增加了聚合物运动这一部分的影响呢。

fhh2626

JCenter 发表于 2025-5-9 20:17
付老师，我这几天在细品您回复的内容，您应该是指考虑了成键、断键、两个原子间的角度等CVs后，不需要考 ...

如果是整体的平动，是不会造成自由能的变化的

JCenter

fhh2626 发表于 2025-5-10 12:54
如果是整体的平动，是不会造成自由能的变化的

我明白了，将两物种的角度作为cv,算是将转动考虑进去了。但是成键的角度应该不等于两者质心间的相对角度，因此，我要么加入质心间的角度作为一个新的cv，要么约束两者间的转动。

kaiyuan

付老师你好，我在尝试用gromacs 2025.1+colvars计算两个疏水分子（石墨烯纳米带加侧链）的二聚化自由能，CV选择的是分子中石墨烯纳米带（重原子+H）之间的质心距离，以及两个分子之间的二面角，用的是WTM-eABF方法，尝试了不同的参数，但是在运行几到几十ns后总是报错：terminate called after throwing an instance of 'gmx::InternalError'
  what():  Freeing of the device buffer failed. CUDA error #700 (cudaErrorIllegalAddress): an illegal memory access was encountered.
唯一一次能够运行时间长一点的是将模拟步长由2fs改成1fs，请问是因为我的CV设置不合理吗还是我设置WTM-eABF的参数有问题
下面是我的colvars设置文件：
colvarsTrajFrequency         5000
colvarsRestartFrequency      25000000
colvar {
    name                     dcom
    width                    0.05
    lowerboundary            0.15
    upperboundary            2.8
    extendedlagrangian       on
    extendedFluctuation      0.05
    extendedTimeConstant     200
    subtractAppliedForce     on
    expandboundaries         on
    distance {
        group1 {
            name plane1
            atomNumbersRange 1-80  #C atoms
            atomNumbersRange 147-170 #H atoms
        }
        group2 {
            name plane2
            atomNumbersRange 271-350
            atomNumbersRange 417-440
        }
    }   
}

colvar {
    name                     tor
    width                    5.0
    lowerboundary            -180
    upperboundary            180
    extendedlagrangian       on
    extendedFluctuation      5.0
    extendedTimeConstant     200
    subtractAppliedForce     on
    #expandboundaries         on
    dihedral {
        group1 { atomNumbers 1 }
        group2 { atomsOfGroup plane1 }
        group3 { atomsOfGroup plane2 }
        group4 { atomNumbers 271 }
    }
}

harmonicWalls {
    name walls
    colvars dcom
    upperWalls 2.5
    forceConstant 2.0
}

abf {
    colvars                 dcom tor
    fullSamples            1000
    historyfreq            25000000
    writeCZARwindowFile
}

metadynamics {
    colvars                 dcom tor
    hillwidth              5
    hillweight             0.5
    welltempered           on
    biastemperature        2000
}

fhh2626

kaiyuan 发表于 2025-5-15 13:34
付老师你好，我在尝试用gromacs 2025.1+colvars计算两个疏水分子（石墨烯纳米带加侧链）的二聚化自由能，CV ...

看看你的轨迹，CVs是否能描述你预期的运动