• 2018-04-25 16:24:25

The following is my notes about the GROMACS Soft-Core Potential for Free Energy Calculation

Theory

The soft-core potential between two states A and B is defined as

Vsc(r)=(1−λ)VA(rA)+λVB(rB)=(1−λ)VA(r)+λVB(r)if ασ=0rA=(ασ6Aλp+r6)1/6rB=(ασ6B(1−λ)p+r6)1/6Vsc(r)=(1−λ)VA(rA)+λVB(rB)=(1−λ)VA(r)+λVB(r)if ασ=0rA=(ασA6λp+r6)1/6rB=(ασB6(1−λ)p+r6)1/6

When α=0α=0 or σ=0σ=0, the soft-core potential degenerates to the original potential. αα is usually about 0.5 and p=1p=1 or p=2p=2. The recommended combination is α=0.5,p=1α=0.5,p=1.

For LJ potential,

V=C12r12−C6r6Vsc=C12[1−λ(ασ6λp+r6)2+λ(ασ6(1−λ)6+r6)2]−C6[1−λασ6λp+r6+λασ6(1−λ)6+r6]V=C12r12−C6r6Vsc=C12[1−λ(ασ6λp+r6)2+λ(ασ6(1−λ)6+r6)2]−C6[1−λασ6λp+r6+λασ6(1−λ)6+r6]

Some special cases

Vsc=0⇒V(rA)=0V(rB)=0rA=rB⇒λ=0.5Vsc=0⇒V(rA)=0V(rB)=0rA=rB⇒λ=0.5

If A state is totally decoupled, which means VA=0VA=0, in this case

Vsc(r)=λVB(rB)=λ{C12(ασ6(1−λ)p+r6)2−C6ασ6(1−λ)p+r6}Vsc(r)=λVB(rB)=λ{C12(ασ6(1−λ)p+r6)2−C6ασ6(1−λ)p+r6}

Define

s6=ασ6rsc=(r6+(1−λ)ps6)1/6minrsc=rwhen λ=1maxrsc=(r6+s6)1/6when λ=0s6=ασ6rsc=(r6+(1−λ)ps6)1/6minrsc=rwhen λ=1maxrsc=(r6+s6)1/6when λ=0

Plot of rscrsc with λλ

rsc/r=(1+(1−λ)p(s/r)6)1/6rsc/r=(1+(1−λ)p(s/r)6)1/6

Plot of maxrscmaxrsc with λλ

maxrsc/r=(1+(s/r)6)1/6maxrsc/r=(1+(s/r)6)1/6

Plot of Vsc(r)Vsc(r)

% <![CDATA[
\alg
V_{sc}(r, s, \l) &= \lambda \left\\{ {C_{12} \over \left(r^6+(1-\lambda)^p s^6\right)^2} -{C_6 \over r^6+(1-\lambda)^p s^6} \right\\} \\
&=4\e \lambda \left\\{ {\s^{12} \over \left(r^6+(1-\lambda)^p s^6\right)^2} -{\s^6 \over r^6+(1-\lambda)^p s^6} \right\\}
\ealg %]]>% <![CDATA[\algV_{sc}(r, s, \l) &= \lambda \left\\{ {C_{12} \over \left(r^6+(1-\lambda)^p s^6\right)^2} -{C_6 \over r^6+(1-\lambda)^p s^6} \right\\} \\ &=4\e \lambda \left\\{ {\s^{12} \over \left(r^6+(1-\lambda)^p s^6\right)^2} -{\s^6 \over r^6+(1-\lambda)^p s^6} \right\\}\ealg %]]>

GMX Options

In GMX, α,σ,pα,σ,p is corresponding to sc_alpha, sc_sigma, sc_power, respectively. The manual said

• sc-alpha: (0) the soft-core alpha parameter, a value of 0 results in linear interpolation of the LJ and Coulomb interactions

• sc-r-power: (6) the power of the radial term in the soft-core equation. Possible values are 6 and 48. 6 is more standard, and is the default. When 48 is used, then sc-alpha should generally be much lower (between 0.001 and 0.003).

• sc-power: (0) the power for lambda in the soft-core function, only the values 1 and 2 are supported

• sc-sigma: (0.3) [nm] the soft-core sigma for particles which have a C6 or C12 parameter equal to zero or a sigma smaller than sc-sigma

The description of sc-sigma is not correct, actually

σ={σ∗=(C12/C6)1/6sc_sigma if C6×C12=0σ={σ∗=(C12/C6)1/6sc_sigma if C6×C12=0

See the discussion

Testing

Script

Figures

Conclusion

• GROMACS DOES support table potential for free energy calculations

• GROMACS calculates the interaction only If r_sc less than rvdw

• The r_sc for any distance should be in the range of the table files, otherwise GROMACS will use a unkown value

• If DispCorr is used, be sure C6 is the correct value

• If sc-alpha is not zero, be sure the σσ GROMACS used is the correct value

• If DispCorr not used, put total interaction in dispersion or repulsion column will force GROMACS to use sc-sigma

• If DispCorr is used and sc-alpha is not zero, be sure C6 and C12 are both correct

◆本文地址: https://jerkwin.github.io/2018/04/25/GROMACS自由能计算中软核势的使用/, 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆