調節平面分子間距的方法
調節平面分子間距的方法
文/Sobereva@北京科音 2013-Feb-3
計算化學研究中有時需要調節兩個平面分子間距離，比如研究不同間距下的相互作用能。如果兩個平面分子的坐標正好對著，或者分子平面平行于某個笛卡爾平面，那么改變間距很容易。前者就用gview的設定鍵長大小的工具就能調，后者就自行對笛卡爾坐標統一加減一定的數值就可以。比較麻煩的是分子既不平行于笛卡爾平面，兩個分子也不是正對著，如下圖這種情況
這里提供兩個辦法設定這兩個分子平面的（垂直）間距。
方法1：
使用protoplane程序（/usr/uploads/file/20150610/20150610055901_52653.rar），得到一個分子上某原子在另一個分子平面上的（垂直）投影點坐標，并在此處設一個虛原子，然后用gview的鍵長設定工具調節這兩個原子的間距（設定窗口中Atom1和Atom2都選為Translate group模式），兩個分子就會相應地整體移動，從而調整了分子平面的間距。
protoplane程序啟動后先依次輸入三個原子的坐標來定義平面方程。然后再輸入要投影的原子的坐標，就會輸出投影后的坐標，以及此原子到平面的垂直距離。
輸入的時候為了省事，可以直接把原子笛卡爾坐標從Gaussian輸入文件里復制到DOS窗口內（先在文本編輯器里復制坐標，然后在DOS窗口標題欄上點右鍵，選編輯-粘貼）。之后，可以將生成的坐標直接拷貝到Gaussian輸入文件里（DOS窗口標題欄上點右鍵，選編輯-復制，然后圈中要復制的區域，點回車，然后粘貼到文本編輯器里），并定義為X原子。
如下所示，將9號原子投影到27、35、41號原子定義的平面上就是53X，之后在gview里調節9C和53X的距離即可改變分子平面的間距了
protoplane程序源代碼如下：
program prjtoplane
implicit real*8 (a-h,o-z)
write(*,*) "prjtoplane: Get the projected position of a point on a specific plane"
write(*,*) "Programmed by Sobereva (sobereva@sina.com), 2013-Feb-2"
write(*,*)
write(*,*) "You will input three point to define the plane"
write(*,*) "Input x,y,z coordinate of point 1 e.g. 2.3,1.0,-4.2 or 2.3 1.0 -4.2"
read(*,*) x1,y1,z1
write(*,*) "Input x,y,z coordinate of point 2"
read(*,*) x2,y2,z2
write(*,*) "Input x,y,z coordinate of point 3"
read(*,*) x3,y3,z3
do while(.true.)
write(*,*) "Input x,y,z coordinate the point you want to project"
read(*,*) x0,y0,z0
call pointprjple(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,prjx,prjy,prjz)
write(*,*)
write(*,*) "The x,y,z coordinate of the projected position is"
write(*,"(3f15.8)") prjx,prjy,prjz
write(*,"(a,f15.8)") " The vertical distance from the point to the plane is",dsqrt((prjx-x0)**2+(prjy-y0)**2+(prjz-z0)**2)
write(*,*)
write(*,*) "Now you can input another point you want to project, or press Ctrl+C to exit"
end do
end program
subroutine pointprjple(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,prjx,prjy,prjz)
real*8 x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,prjx,prjy,prjz,A,B,C,D,t
call pointABCD(x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D)
t=(D+A*x0+B*y0+C*z0)/(A**2+B**2+C**2)
prjx=x0-t*A
prjy=y0-t*B
prjz=z0-t*C
end subroutine
subroutine pointABCD(x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D)
real*8 v1x,v1y,v1z,v2x,v2y,v2z,x1,y1,z1,x2,y2,z2,x3,y3,z3,A,B,C,D
v1x=x2-x1
v1y=y2-y1
v1z=z2-z1
v2x=x3-x1
v2y=y3-y1
v2z=z3-z1
A=v1y*v2z-v1z*v2y
B=-(v1x*v2z-v1z*v2x)
C=v1x*v2y-v1y*v2x
D=A*(-x1)+B*(-y1)+C*(-z1)
end subroutine
方法2：
這個方法是先讓指定的分子平面平行于XY平面，然后自行對分子的Z坐標統一進行適當加減，就改變了分子平面間距。
2021-Aug-16注：下文用VMD的做法已經不推薦使用了！因為使用Multiwfn來實現方便得多，見《Multiwfn中非常實用的幾何操作和坐標變換功能介紹》（http://www.shanxitv.org/610）。讓體系的一個平面平行于某笛卡爾平面的做法在此文2.4節專門有示例。在此文介紹的Multiwfn的子功能300的子功能7里用選項1還可以方便地讓體系根據特定矢量進行平移。
這里提供一個VMD程序的Tcl腳本用于實現此目的。先將以下內容復制到VMD的命令行窗口執行
proc alignplane {ind1 ind2 ind3} {
set atm1 [atomselect top "serial $ind1"]
set atm2 [atomselect top "serial $ind2"]
set atm3 [atomselect top "serial $ind3"]
set vec1x [expr [$atm2 get x] - [$atm1 get x]]
set vec1y [expr [$atm2 get y] - [$atm1 get y]]
set vec1z [expr [$atm2 get z] - [$atm1 get z]]
set vec2x [expr [$atm3 get x] - [$atm1 get x]]
set vec2y [expr [$atm3 get y] - [$atm1 get y]]
set vec2z [expr [$atm3 get z] - [$atm1 get z]]
set sel [atomselect top all]
$sel move [transvecinv [veccross "$vec1x $vec1y $vec1z" "$vec2x $vec2y $vec2z"]]
$sel move [transaxis y 90]
}
然后在命令行窗口運行諸如alignplane 27 35 41就可以讓27、35、41號原子定義的平面平行于XY平面，如下所示
現在可以利用file-save coordinates將當前結構保存出來，然后用諸如Excel/Origin之類或自寫程序來平移分子。但是在VMD里也可以很方便地整體平移分子。
VMD會將自動將每個分子以fragment x描述。諸如前例，就有fragment 0和fragment 1代表這兩個分子。假設要平移的是fragment 1，就先執行
set sel [atomselect top "fragment 1"]
然后比如向Z軸正方向平移0.5埃，即平移向量是0 0 0.5，就運行$sel move [transoffset {0 0 0.5}]