|
2018-10-07 22:23:47
生成分子初始构型的时候, 经常需要对分子进行随机旋转. 这并不是一个困难的工作, 因为早就有成熟的方法了. 在这里我展示几种不同的实现方法, 供参考.
下面是四种不同方法分别得到的球面上的随机点
Fig.1
Fig.2
Fig.3
Fig.4
前面三种是通常所谓的随机点, 最后的这种点可用于一些特殊的情况. 这种点涉及到的知识很多, 这里就不多说了. 如果需要请参考后面的资料.
bash | |
---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 | echo "1molH 0 0 5" | awk 'BEGIN{ Pi=4*atan2(1,1); TwoPi=2*Pi method=3; Nmax=1000}{ Natm=$1; getline for(i=1; i<=Natm; i++) { getline; Satm[i]=$1; Xatm[i]=$2; Yatm[i]=$3; Zatm[i]=$4 } print Nmax*Natm"\nRot. Method: "method for(n=1; n<=Nmax; n++) { isang=getAxi(method, axi, n, Nmax) if(isang) { tht=Pi x=.5*axi[1] y=.5*axi[2] z=.5*(axi[3]+1) } else { tht=axi[0] x=axi[1] y=axi[2] z=axi[3] } AxiRotMat(isang, tht, x, y, z) for(i=1; i<=Natm; i++) { V[1]=Xatm[i]; V[2]=Yatm[i]; V[3]=Zatm[i] RotVect(RotMat, V) printf "%s %9.6f %9.6f %9.6f\n", Satm[i], V[1], V[2], V[3] } } exit}function acos(x) { return atan2(sqrt(1-x*x),x) }function getAxi(method, axi, Ipnt, Npnt, i,r1,r2,tht1,tht2) { isang=1; ;axi[0]=0 if(method==1) { # Marsaglia Method while(1) { r1=2*rand()-1 r2=2*rand()-1 if(r1*r1+r2*r2<1) break } tht1=r1*r1+r2*r2 axi[1]=2*r1*sqrt(1-tht1) axi[2]=2*r2*sqrt(1-tht1) axi[3]=1-2*tht1 } else if(method==2) { # theta distributuion Method tht1=acos(2*rand()-1) # tht1=rand()*Pi # can NOT use tht2=rand()*TwoPi axi[1]=sin(tht1)*cos(tht2) axi[2]=sin(tht1)*sin(tht2) axi[3]=cos(tht1) } else if(method==3) { isang=0 # 四元数方法 tht1=rand()*TwoPi tht2=rand()*TwoPi r1=sqrt(1-rand()) r2=sqrt(1-r1*r1) axi[0]=sin(tht1)*r1 axi[1]=cos(tht1)*r1 axi[2]=sin(tht2)*r2 axi[3]=cos(tht2)*r2 } else if(method==4) { # Sprial Method, for regular equidistributed Ipnt-- r1 = -1+2*(Ipnt+0.5)/Npnt # +0.5 表示在厚度的中点处取点 tht2 = Pi*(sqrt(5)-1)*Ipnt r2 = sqrt(1-r1*r1) axi[1] = r2*cos(tht2) axi[2] = r2*sin(tht2) axi[3] = r1 } return isang}function RotVect(RotMat, V) { # 按照矩阵RotMat旋转矢量V Vx = V[1]; Vy = V[2]; Vz = V[3] V[1]= RotMat[1,1]*Vx + RotMat[1,2]*Vy + RotMat[1,3]*Vz V[2]= RotMat[2,1]*Vx + RotMat[2,2]*Vy + RotMat[2,3]*Vz V[3]= RotMat[3,1]*Vx + RotMat[3,2]*Vy + RotMat[3,3]*Vz}function AxiRotMat(isang, tht, x, y, z, w,rsq) { # 绕矢量(x,y,z)旋转tht弧度,四元数方法 if(isang) { # w 为角度, 使用 轴-角 计算旋转矩阵 w=cos(tht/2) rsq=sqrt(x*x+y*y+z*z); if(rsq>0) rsq = sin(tht/2)/rsq x *= rsq; y *= rsq; z *= rsq } else { # 否则直接使用传入的四元数构建 w=tht rsq=sqrt(x*x+y*y+z*z+w*w); if(rsq>0) rsq = 1/rsq x *= rsq; y *= rsq; z *= rsq; w *= rsq } RotMat[1,1]= 1-2*(y*y+z*z); RotMat[1,2]= 2*(x*y-w*z); RotMat[1,3]= 2*(x*z+w*y); RotMat[2,1]= 2*(x*y+w*z); RotMat[2,2]= 1-2*(x*x+z*z); RotMat[2,3]= 2*(y*z-w*x); RotMat[3,1]= 2*(x*z-w*y); RotMat[3,2]= 2*(y*z+w*x); RotMat[3,3]= 1-2*(x*x+y*y); if(x*x+y*y+z*z+w*w==0) { RotMat[1,1]=-1; RotMat[2,2]=-1; RotMat[3,3]=-1 }}' |
uniform random rotations, Ken Shoemake 1993; Graphics Gems 3
algorithms - uniform distribution of points on the surface of a sphere
When Random Numbers Are Too Random: Low Discrepancy Sequences