通过colvars结合lammps以eABF方法计算平行板石墨烯之间PMF的若干问题

dumsky2019928

大家好，我想请教一下关于计算PMF的若干问题。
我的目的是计算浸泡在lj流体（即非键力仅存在范德华相互作用，没有库仑力）中的平行石墨烯之间的平均力势，体系如图1所示。其中石墨烯的xy方向通过fix spring/self施加了限制势，以维持质心的连线与z轴平行。反应路径定义为石墨烯质心之间的距离，数值为1埃到31埃。我参照帖子http://bbs.keinsci.com/thread-10805-1-1.html设置了以eABF方法计算上述反应路径下PMF的colvars脚本，对计算结果以及一些细节存在疑惑：
1.在colvars的输出文件中，有out.pmf与out.czar.pmf，其中哪个才是eABF计算得到的PMF？
2.out.pmf与out.czar.pmf的结果如图2和图3，不论是out.pmf还是out.czar.pmf，在曲线的最左侧，即石墨烯质心之间的距离很小时，均存在一个看似不合理“平台”，这里ΔG的数值本应很大，应该如何解决这一错误呢？
3.PMF曲线是否应该减去限制势对石墨烯做的功？
colvars脚本和lammps输入文件也附于图片之后，期待dalao的回答。

dumsky2019928

colvars的输入脚本如下：



colvarsTrajFrequency    1000
colvarsRestartFrequency 1000
units real
indexFile gra.ndx
colvar {
name r
width 0.1
lowerboundary 1
upperboundary 31
extendedlagrangian   on
extendedFluctuation  0.1
extendedTimeConstant   200
distance {
   group1 {
      indexGroup graphene1
    }
   group2 {
      indexGroup graphene2
    }
   
  }
}
harmonicWalls {
name mywalls
colvars r
lowerWalls 1 1
upperWalls 31 31
forceConstant 20
}

abf {
colvars r
fullSamples 1000
historyfreq  1000000
writeCZARwindowFile
}

dumsky2019928

lammps的in文件如下：



# ----------------- Init Section -----------------

dimension    3
units real
boundary p p p
neighbor 2 multi
neigh_modify every 1 delay 0 check yes# page 50000 one 5000
atom_style   full
pair_style lj/charmm/coul/charmm 9.0 10.0
pair_modify mix arithmetic
bond_style harmonic
angle_style harmonic
dihedral_style fourier
improper_style cvff
special_bonds amber
#kspace_style pppm 1.0e-4

# ----------------- Atom Definition Section -----------------

#read_data "system.data"
read_restart "xiafeng.restart"

# ----------------- Settings Section -----------------
pair_coeff 1 1 0.03377 3.91
bond_coeff 1 461.1 1.3984
angle_coeff 1 66.6 120.0
dihedral_coeff 1 1 3.625 2 180.0
improper_coeff 1 1.1 -1 2
pair_coeff 2 2 0.03377 3.91
bond_coeff 2 461.1 1.3984
angle_coeff 2 66.6 120.0
dihedral_coeff 2 1 3.625 2 180.0
improper_coeff 2 1.1 -1 2
pair_coeff 3 3 0.08882803 3.5

group liquid type 3
group graphene1 type 1
group graphene2 type 2
group graphene type <> 1 2
#fix freeze graphene setforce 0.0 0.0 0.0

#region lwallsy block EDGE EDGE -11.88 11.88 -22.5 -2.5 open 1 open 2 open 5 open 6 units box
#region lwallsz block EDGE EDGE -11.88 11.88 -22.5 -2.5 open 1 open 2 open 3 open 4 open 6 units box
#
#region rwallsy block EDGE EDGE -11.88 11.88 2.5 22.5 open 1 open 2 open 5 open 6 units box
#region rwallsz block EDGE EDGE -11.88 11.88 2.5 22.5 open 1 open 2 open 3 open 4 open 5 units box
#
#fix lywall liquid wall/region lwallsy lj93 0.03377 3.91 10
#fix lzwall liquid wall/region lwallsz lj93 0.03377 3.91 10
#
#fix rywall liquid wall/region rwallsy lj93 0.03377 3.91 10
#fix rzwall liquid wall/region rwallsz lj93 0.03377 3.91 10                                                     

#fix zrest graphene spring/self 10000.0
fix zrest1 graphene1 spring/self 10.0 xy
fix zrest2 graphene2 spring/self 10.0 xy
fix_modify zrest1 energy yes
fix_modify zrest2 energy yes
thermo 1000
reset_timestep 0
timestep 1
fix int all nve
fix lt liquid temp/csvr 298 298     100 1919
fix gt1 graphene1 temp/csvr 298 298 100 810
fix gt2 graphene2 temp/csvr 298 298 100 114              
dump eqtraj all custom 5000 eq.dump id type x y z
#fix pmf1 graphene1 ave/time 1 100000 100000 f_zrest1 file pmf1.txt mode scalar
#fix pmf2 graphene2 ave/time 1 100000 100000 f_zrest2 file pmf2.txt mode scalar
fix eabf all colvars colvars.in seed 114514 tstat gt1
#fix pmf all ave/time 1 50000 50000 f_zrest file pmf.txt mode scalar
compute g1temp graphene1 temp
compute g2temp graphene2 temp
compute ltemp liquid temp
#min_style quickmin
#minimize  1.0e-4 1.0e-6 10000000 1000000000
thermo_style custom time c_ltemp c_g1temp c_g2temp f_zrest1 f_zrest2 press
restart 100 xiafeng.restart xiafeng.restart
run 50000000
#group2ndx gra.ndx graphene1 graphene2
write_restart xiafeng.restart
write_data xiafeng.data

fhh2626

1、看out.czar.pmf
2、那个区域是没采到样，可以直接截掉（当然也可以耐心等那个区域采到样，如果你想的话）
3、如果是在XY方向上的限制的话，应该对你的反应坐标方向不做功

dumsky2019928

fhh2626 发表于 2023-7-27 21:32
1、看out.czar.pmf
2、那个区域是没采到样，可以直接截掉（当然也可以耐心等那个区域采到样，如果你想的话 ...

谢谢大佬！也就是说，如果跑的时间再长一些，比如跑100ns，这样或许就能使平台的地方充分采样吧？
还有就是反应路径应该选为距离在z方向的投影才更合理

fhh2626

dumsky2019928 发表于 2023-7-29 11:05
谢谢大佬！也就是说，如果跑的时间再长一些，比如跑100ns，这样或许就能使平台的地方充分采样吧？
还有 ...

这个不一定，如果范围设置得特别不合理的话可能跑多久都采不到样。如果合理的话可能过一段时间就会有采样

如果完全限制（constrain）XY方向的运动的话，选distance和distanceZ结果应该是一样的

1227962972

dumsky2019928 发表于 2023-7-27 20:06
colvars的输入脚本如下：

请问有没有一份完整的lammps.in文件，colvars文件和indexFile gra.ndx文件样本，我想拿来学习一下怎么写。