NAMD二聚体的单链返回镜像，周期性边界条件设置求助

难破船

用namd跑一个二聚体蛋白（~300个氨基*2），水盒子大小为蛋白向外延申10a，进行了5000步能量最小化后进行模拟。设置了 periodic boundary conditions（pbc），模拟进行到一定时间，当其中一条单链达到pbc界限时，就会返回镜像。理想结果是整个二聚体达到pbc界限才返回镜像，而非其中的单链。图1-2大概表述了该问题。

尝试：
1）扩大水盒子大小，增大到二聚体的两倍，约向外延申了25A；
2）怀疑是构建的结构有问题，于是用实际测得的晶体结构（模板蛋白，相似度~80%），进行了相同的模拟；
均失败，总是出现单链返回镜像。
3）现在的笨办法是找出这些错误的帧，直接把单链移动回来。

哪位能给我些建议，如何设置pbc，才能使二聚体作为一个整体返回镜像吗？万分感谢！


###############################################
## self-defined variable
set struct         dimer_wbtwice_n
set last_run       dimer_wbtwice_n_min
set this_run       dimer_wbtwice_n_md_57
set fts            0 ;# firsttimestep
set nmd            10000000 ;# 20ns
set tst            2.0 ;# 2fs/timestep
set tep            330.15 ;# temperature
set freq           1000 ;# 1000 steps = every 2ps


## input and output filename prefix
structure           ../common/$struct.psf
coordinates         $last_run.coor
extendedsystem      $last_run.xsc
outputname          $this_run
binaryoutput        no


## force field
paratypecharmm      on
parameters          ../common/par_all36_prot.prm
parameters          ../common/par_all36_carb.prm
parameters          ../common/par_all36_cgenff.prm
parameters          ../common/par_all36_lipid.prm
parameters          ../common/par_all36_na.prm
parameters          ../common/toppar_water_ions_namd.str
exclude             scaled1-4
1-4scaling          1.0
cutoff              12.0
switching           on
switchdist          10.0 ;# cutoff - 2.0
pairlistdist        14.0 ;# cutoff + 2.0
margin              2.0

## pbc amd pme

wrapall               on
cellbasisvector1      114.56        0.0        0.0
cellbasisvector2      0.0        158.22        0.0
cellbasisvector3      0.0        0.0        152.76
cellorigin            21.45        35.03        14.39

pme                   yes
pmegridspacing        1.0

## ensemble: npt

temperature          $tep
commotion            no
timestep             $tst
rigidbonds           all
firsttimestep        $fts

fullelectfrequency   1
nonbondedfreq        1
stepspercycle        20

## constant temperature control
tcouple               on
tcoupletemp           $tep

## constant pressure control
langevinpiston        on
langevinpistontarget  1.01325 ;# in bar -> 1 atm
langevinpistonperiod  100.0
langevinpistondecay   50.0
langevinpistontemp    $tep

usegrouppressure      yes ;# needed for rigidbonds
useflexiblecell       no
useconstantarea       no

## output frequency
restartfreq           $freq
dcdfreq               $freq
xstfreq               $freq
outputenergies        $freq
outputpressure        $freq

## execution
reinitvels            $tep
run                   $nmd
#####################################

FIG. 1 frame0 为初始结构；frame337 时，单链A（黄色）达到了 PBC 界限；frame338 时，链A 返回了镜像。这导致了 FIG2 中的 RMSD 从 2A 突然激增。
理想结果：frame338 链A 不返回镜像，而是二聚体作为整体达到 PBC 界限才返回镜像。


FIG.1


FIG.2

霜晨月

这个可能只是显示的问题，对实际计算过程和计算结果没有影响。要先修正PBC再计算RMSD，才有意义。

fhh2626

你把轨迹读进去，pbc wrap一下就可以了
https://www.ks.uiuc.edu/Research/vmd/plugins/pbctools/

（不要问我参数怎么设，wrap是门玄学，自己多试试，不过我还没见过wrap不成功的）

难破船

霜晨月 发表于 2018-12-25 12:20
这个可能只是显示的问题，对实际计算过程和计算结果没有影响。要先修正PBC再计算RMSD，才有意义。

Hi，感谢回复。

修正PBC可以用命令，或者写个小脚本实现。
我疑惑的是进行的是多聚体的模拟，如果单链返回镜像，就像拆成了多条链的同时模拟，这在计算链与链之间的一些参数时，不会影响到嘛。
所以想着有没有可能，在模拟过程中就把多聚体作为一个整体，避免出现单链返回镜像。

yl1088

rmsd是align之后才计算的，不是直接就计算

难破船

yl1088 发表于 2018-12-25 13:38
rmsd是align之后才计算的，不是直接就计算

是的，要先align，我知道的。
就是想知道有没有方法可以在模拟时，就直接把多聚体作为一个整体，而不是单链返回镜像。

难破船

fhh2626 发表于 2018-12-25 12:47
你把轨迹读进去，pbc wrap一下就可以了
https://www.ks.uiuc.edu/Research/vmd/plugins/pbctools/

pbc wrap 在弄啦，参数设置还不太会，摸索ing。
谢谢你~

霜晨月

这种情况，只是看上去两条链分开了，其实计算时还是按照两条链在一起来计算的，不用担心。你算一下frame338两条链间的作用能，与前后帧对比一下就知道了。

难破船

霜晨月 发表于 2018-12-26 00:35
这种情况，只是看上去两条链分开了，其实计算时还是按照两条链在一起来计算的，不用担心。你算一下frame338 ...

谢谢你😀