三维空间随机旋转

类别:    标签: 数理   阅读次数:   版权: (CC) BY-NC-SA

生成分子初始构型的时候, 经常需要对分子进行随机旋转. 这并不是一个困难的工作, 因为早就有成熟的方法了. 在这里我展示几种不同的实现方法, 供参考.

下面是四种不同方法分别得到的球面上的随机点

Marsaglia方法


视图: 投影 正交    速度:
模型: 球棍 范德华球 棍状 线框 线型   名称
左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动

Fig.1

球坐标分布方法: 生成正确分布的球坐标


视图: 投影 正交    速度:
模型: 球棍 范德华球 棍状 线框 线型   名称
左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动

Fig.2

四元数方法: 生成随机的单位四元数


视图: 投影 正交    速度:
模型: 球棍 范德华球 棍状 线框 线型   名称
左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动

Fig.3

螺线方法: 近似均匀分布的球面点


视图: 投影 正交    速度:
模型: 球棍 范德华球 棍状 线框 线型   名称
左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动

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 "1
mol
H 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 }
}
'

参考资料

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: gnuplot显示优化过程
后一篇: AMBER 2015伦敦研讨会 - 组合聚类分析示例

访问人次(2015年7月 9日起): | 最后更新: 2018-10-14 22:15:41 UTC | 版权所有 © 2008 - 2018 Jerkwin