xyz2DRC (xyz 座標からDRC計算の入力データの作るプログラム) の改良
回転と推進を加えた方向を作る方法
・重心をもとめる (xx,yy,zz)
2原子間の距離をもとめ、その2原子の重さの比より 2点の重心をもとめ、繰り返す。
ソース
xx=x(nn1)
yy=y(nn1)
zz=z(nn1)
aaa=aa(nn1)
c
do i =nn1,nn2-1
bbb = aaa +aa(i+1)
xx= aaa/bbb * xx + aa(i+1)/bbb * x(i+1)
yy= aaa/bbb * yy + aa(i+1)/bbb * y(i+1)
zz= aaa/bbb * zz + aa(i+1)/bbb * z(i+1)
aaa= bbb
end do
・回転方向を決める。
適当な値を自分で入力( ax,ay,az )
・回転方向の速度を決める。(bx,by,bz)
ソース
do i = nn1,nn2
t=((x(i)-xx)*ax+(y(i)-yy)*ay+(z(i)-zz)*az)/(ax**2+ay**2+az**2)
tx(i) = ax * t + xx ← 原子の座標から回転軸への垂線と
ty(i) = ay * t + yy ← 回転軸の交点の座標
tz(i) = az * t + zz ←
r(i) = sqrt((tx(i)-x(i))**2+(ty(i)-y(i))**2+(tz(i)-z(i))**2)
if(r(i).gt.0.05) then ↑ 回転軸との距離
bx(i) = ay * z(i) - az * y(i)
by(i) = az * x(i) - ax * z(i) ← 回転軸ベクトルと原子座標との外積
bz(i) = ax * y(i) - ay * x(i) ( 回転方向のベクトルが求まる )
else
bx(i) = 0.0d0 半径の小さい物は省略
by(i) = 0.0d0
bz(i) = 0.0d0
end if
bx(i) = bx(i) * r(i) 半径をかけて、回転の大きさの比を決める。
by(i) = by(i) * r(i) ( のちに正規化し、最大値を1にする )
bz(i) = bz(i) * r(i)
end do
・正規化
半径が一番大きい原子の方向ベクトルの大きさを 1にする。
ソース
rmax = r(nn1)
nn3 = nn1
do i = nn1+1,nn2
if (r(i).gt.rmax) then
rmax = r(i) 半径の最大値を求める
nn3 = i その時の "i"
end if
end do
c
ttt=sqrt(bx(nn3)**2+by(nn3)**2+bz(nn3)**2)
do i=nn1,nn2
bx(i) = bx(i) / ttt
by(i) = by(i) / ttt
bz(i) = bz(i) / ttt
end do
c
・動かさない原子の振動
乱数を与え、大きさを"0.01" にする
do i = 1,nn1-1
bx(i) = rand()
by(i) = rand()
bz(i) = rand()
ttt=sqrt(bx(i)**2+by(i)**2+bz(i)**2)
ttt = ttt / 0.01d0
bx(i) = bx(i) / ttt
by(i) = by(i) / ttt
bz(i) = bz(i) / ttt
end do
先週は上のプログラム等を作成した。
次の課題
入力するデータは、"KINETIC" の値でなく原子の温度を入力し、
その場合の原子の速度を計算し出力させる。
再記入(〜10/15) tex ファイルにまとめ中
mosic なら表示されるのですが.....