任意の温度を入力し、その温度でナノチューブの構造最適化を実行するプログラム。
このプログラムでは、結果が出力されるファイルはあらかじめ決まっており(result.xyz)、 アニメーション表示に対応したデータが出力される。
c
c TersoffpotentialforCarbonsystem
c (1998/Jun/8) By R.Matsuo.
c (2000/Nov/8) restoredbyT.Nakamura
c 入力ファイル.xyz 出力ファイルresult.xyz
c このプログラムは別枠のパラメーターファイルを必要とする。
c
c 温度を任意のステップごとに入力し、計算させるプログラム。
PROGRAMMD
INCLUDE'PARAMETER'
C 原子座標のFILE名変数
CHARACTER*50FILENAME
CHARACTER*50FILENAME2
C 原子数の保持変数
INTEGERN
C loop用変数
INTEGERI,J,K
C 温度の変数(k)
REAL*8H,G
REAL*8TEMP1
REAL*8TEMP2
c 繰り返し変数(steps)
INTEGERA,B,M
C ALPHA、VV2の宣言vv2(2)1:C602:nanotube
REAL*8ALPHA1,ALPHA2
DIMENSIONVV2(2)
REAL*8VV2
REAL*8SQRT
C 座標用配列
REAL*8 ZAHYO(MAX_ATOM3)
CHARACTER*8 AN(MAX_ATOM)
C 近接リスト配列
REAL*8LIST(0:MAX_LIST_N2,MAX_ATOM)
C 第二2近接リスト配列
REAL*8LIST2(0:MAX_LIST2_N,MAX_ATOM)
C 近接リスト配列
REAL*8LIST_VDW(0:MAX_LIST_VN2,MAX_ATOM)
C 近接データ配列IDXでリストから参照される
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
REAL*8 DIFF1(MAX_ATOM3)
REAL*8 VELO(MAX_ATOM3)
REAL*8 DD(MAX_ATOM3)
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
C TERSOFF関数の型
REAL*8TERSOFF
C 系のエネルギー
REAL*8ENERGY,ENERGY_VDW
INTEGERIDX
c itest=1debug
itest=0
C FILE名を取得する
CALLGETFILENAME(FILENAME)
open(STATUS='OLD',unit=20,file='tempmd.log')
write(20,25)FILENAME
25 format(a50)
C XYZ座標を読み込む
CALLREAD_ZAHYO(FILENAME,N,I,ZAHYO,AN)
write(*,*)'I=',I
DO29J=1,MAX_ATOM3
VELO(J)=0.0D0
29 CONTINUE
IDX=0
CALLMAKE_LIST(N,ZAHYO,LIST,KYORI,CUTOFF,IDX)
CALLMAKE_LIST2(N,LIST,LIST2)
CALLMAKE_LIST_VDW(N,ZAHYO,LIST_VDW,KYORI,CUTOFF,IDX,LIST,
#LIST2,AN)
ENERGY=TERSOFF(N,ZAHYO,LIST,KYORI,CUTOFF)
ENERGY_VDW=VDW(N,ZAHYO,LIST_VDW,KYORI,CUTOFF)
ENERGY=ENERGY+ENERGY_VDW
WRITE(*,*)'何回繰り返す?A='
READ(*,*)A
WRITE(*,*)'TEMPARATUREofC_60(Real)?H(K)='
READ(*,*)H
WRITE(*,*)'TEMPARATUREofNANO(Real)?G(K)='
READ(*,*)G
WRITE(*,*)'N*3',N*3
WRITE(*,*)'C_60 幾つ入ってるB=?'
READ(*,*)B
WRITE(*,*)'何回ごとに温度入力M=?'
READ(*,*)M
K=0
DO41I=1,A
write(*,*)'I=',I
C 近接リストを生成
IDX=0
CALLMAKE_LIST(N,ZAHYO,LIST,KYORI,CUTOFF,IDX)
C 第二近接リストを生成
CALLMAKE_LIST2(N,LIST,LIST2)
c WRITE(*,*)IDX
c WRITE(*,*)'make_list_vdw'
CALLMAKE_LIST_VDW(N,ZAHYO,LIST_VDW,KYORI,CUTOFF,IDX,LIST,
#LIST2,AN)
c WRITE(*,*)IDX
c WRITE(*,*)'make_newkyori'
CALLMAKE_NEWKYORI(N,ZAHYO,LIST,KYORI,CUTOFF)
c WRITE(*,*)'make_newkyoriv'
CALLMAKE_NEWKYORIV(N,ZAHYO,LIST_VDW,KYORI,CUTOFF)
c WRITE(*,*)'calc_diff'
CALLCALC_DIFF1(N,ZAHYO,LIST,KYORI,CUTOFF,DIFF1)
c WRITE(*,*)'calc_diff_vdw'
CALLCALC_DIFF1_VDW(N,ZAHYO,LIST_VDW,KYORI,CUTOFF,DIFF1)
c 共役勾配法
c DO38J=1,MAX_ATOM3
c DIFF2(0,J)=0
c DD(J)=0.0D0
c38 CONTINUE
c CALLCALC_DIFF2(N,ZAHYO,LIST,LIST2,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
c CALLCALC_DIFF2_VDW
C # (N,ZAHYO,LIST_VDW,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
c CALLCONJGRAD(N,ZAHYO,DIFF1,DIFF2,DIFF2_DATA)
c
c 分子動力学法
DO39J=1,N*3
VELO(J)=VELO(J)+DIFF1(J)*CONV_VELO
ZAHYO(J)=ZAHYO(J)+VELO(J)
39 CONTINUE
ccc
c 温度入力
ccc
c
c c_60とnanotubeを別けてそれぞれに温度を振り分ける。
c 読み込む ファイル の上から60番目までがc_60
c JJ=1C60,JJ=2NANOTUBE
c c_60
c
c ここでM回ごとにこのループに入るように条件IF文を書く。
c IF文とMOD構文を用いる。
c
IF(MOD(I,M).EQ.0) THEN
JJ=1
NNC60=60*B
NNC603=NNC60*3
VV2(JJ)=0.0
DOJ=1,NNC603
VV2(JJ)=VV2(JJ)+VELO(J)**2
ENDDO
WRITE(*,*)'VV2OFC60=',VV2(JJ)
C
ALPHA1=SQRT(2.08D-7*REAL(NNC60)*H/VV2(JJ))
C ---計算から求めたalphaをここで入力する。
---C ここではC_60用の入力画面。
C
DOJ=1,NNC603
VELO(J)=ALPHA1*VELO(J)
ENDDO
C
C次はnanotube
C
JJ=2
VV2(JJ)=0.0
DOJ=NNC603+1,N*3
VV2(JJ)=VV2(JJ)+VELO(J)**2
ENDDO
WRITE(*,*)'VV2ofnanotube=',VV2(JJ)
ALPHA2=SQRT(2.08D-7*REAL(N-NNC60)*G/VV2(JJ))
C
C ---nanotube用にここでalphaを入力する。
---C
DOJ=NNC603+1,N*3
VELO(J)=ALPHA2*VELO(J)
ENDDO
c
cALPHA=SQRT(3*N*Kb*T/m*VV2)
c =SQRT(3*N*1.38D-23*T/1.99D-26*VV2)
c =SQRT(2.08D3*N*T/VV2)
c m/s とA/fsを比較して、係数を10D-10倍する。
c
cALPHAの画面表示
WRITE(*,*)'alpha1,2',alpha1,alpha2
TEMP1=0.48D7*VV2(1)*ALPHA1/REAL(NNC60)
TEMP2=0.48D7*VV2(2)*ALPHA2/REAL(N-NNC60)
WRITE(*,*)'TEMP1,2=',TEMP1,TEMP2
ELSE
ENDIF
cファイルへの出力(.log)
OPEN(STATUS='UNKNOWN',UNIT=30,FILE='tempC_60.log')
write(30,*)TEMP1
OPEN(STATUS='UNKNOWN',UNIT=31,FILE='temp_NANO.log')
write(31,*)TEMP2
OPEN(STATUS='UNKNOWN',UNIT=32,FILE='alphaC_60.log')
write(32,*)alpha1
OPEN(STATUS='UNKNOWN',UNIT=33,FILE='alpha_NANO.log')
write(33,*)alpha2
OPEN(STATUS='UNKNOWN',UNIT=34,FILE='VV2C_60.log')
write(34,*)VV2(1)
OPEN(STATUS='UNKNOWN',UNIT=35,FILE='VV2_NANO.log')
write(35,*)VV2(2)
ENERGY=TERSOFF(N,ZAHYO,LIST,KYORI,CUTOFF)
ENERGY_VDW=VDW(N,ZAHYO,LIST_VDW,KYORI,CUTOFF)
ENERGY=ENERGY+ ENERGY_VDW
WRITE(20,40)I,ENERGY
40 format(i6,f20.7)
CALLPRINT_ZAHYO(N,J,ZAHYO,I,ENERGY,DIFF1,AN)
41 CONTINUE
C温度読み込みの確認画面表示
WRITE(*,*)'TEMP=',T
close(20)
STOP
END
c
cThisprogramisformdtorunonDecFortran
c
CC
C 座標データのFILE名を取得する サブルーチン
CC
101 SUBROUTINEGETFILENAME(FILENAME)
CHARACTER*50FILENAME
INTEGERI
I=IARGC()
IF(I.NE.1)then
WRITE(*,*)'Usage:opthoge.xyz'
GOTO110
ENDIF
CALLGETARG(1,FILENAME)
RETURN
110 STOP
END
121 SUBROUTINEGETARGU(LINE,I,ARGU)
CHARACTER*80LINE
CHARACTER*80DUM
CHARACTER*80ARGU
INTEGERJ
INTEGERK
INTEGERlenline
DUM=LINE
do128K=1,I
lenline=len(DUM)
DO122J=1,lenline
IF(DUM(J:J).EQ.'')THEN
GOTO122
ENDIF
IF(DUM(J:J).EQ.char(9) )THEN
GOTO122
ENDIF
IF(DUM(J:J).EQ.char(0) )THEN
GOTO122
ENDIF
GOTO123
122 CONTINUE
123 DUM=DUM(J:lenline)
lenline=len(DUM)
DO124J=2,lenline
IF(DUM(J:J).EQ.'')THEN
GOTO125
ENDIF
IF(DUM(J:J).EQ.char(9) )THEN
GOTO125
ENDIF
IF(DUM(J:J).EQ.char(0) )THEN
GOTO125
ENDIF
124 CONTINUE
125 ARGU=DUM(1:J-1)
DUM=DUM(J:lenline)
128 continue
RETURN
c129 STOP
END
130 REAL*8FUNCTIONATOR(STR)
CHARACTER*80STR
INTEGERI
INTEGERJ
INTEGERS
INTEGERD
INTEGERDCU
INTEGER*4R
INTEGER*4E
I=1
S=1
R=0
EC=0
E=0
IF(STR(1:1).EQ.'+')THEN
S=1
I=I+1
ENDIF
IF(STR(1:1).EQ.'-')THEN
S=-1
I=I+1
ENDIF
DO140J=I,80
IF(STR(J:J).EQ.'.')THEN
EC=1
GOTO140
ENDIF
D=ICHAR(STR(J:J))-ICHAR('0')
IF((D.GT.9).OR.(D.LT.0))THEN
GOTO145
ENDIF
R=R*10+D
E=E+EC
IF(R.GT.9999999)THEN
GOTO145
140 CONTINUE
145 ATOR=1.0D0 *S*R/(10.0D0**E)
RETURN
c149 STOP
END
CC
C 座標データ取り込み サブルーチン
CC
201 SUBROUTINEREAD_ZAHYO(FILENAME,N,I,ZAHYO,AN)
INCLUDE'PARAMETER'
CHARACTER*50FILENAME
C 原子数を保持する変数
INTEGERN
INTEGERN3
C FILEIOチェック用変数
INTEGERIOCHECK
C loopvariable
INTEGERI
C 炭素原子の座標用配列
REAL*8ZAHYO(MAX_ATOM3)
C 元素名用変数
CHARACTER*8 AN(MAX_ATOM)
CHARACTER*80LINE,ARGU
C FILEOPEN
OPEN(60,FILE=FILENAME,STATUS='OLD',IOSTAT=IOCHECK)
C FILEOPENのエラーチェック
IF(IOCHECK)THEN
WRITE(*,*)'Fileopenerror.'
GOTO299
ENDIF
C 原子数の制限チェック
READ(60,*)N
IF(MAX_ATOM.LT.N)THEN
WRITE(*,*)"Toomanyatoms."
GOTO299
ENDIF
N3=N*3
READ(60,209)LINE
208 format(f14.8)
209 format(a80)
C 座標読み込み 単位はÅ
DO210I=1,N3,3
READ(60,209)LINE
CALLGETARGU(LINE,1,ARGU)
AN(I/3+1)=ARGU
CALLGETARGU(LINE,2,ARGU)
ZAHYO(I)=ATOR(ARGU)
CALLGETARGU(LINE,3,ARGU)
ZAHYO(I+1)=ATOR(ARGU)
CALLGETARGU(LINE,4,ARGU)
ZAHYO(I+2)=ATOR(ARGU)
C 読み込んだ座標の画面表示（確認）
itest=0
if(itest.eq.1)then
WRITE(*,*)'X,Y,Z=',ZAHYO(I),ZAHYO(I+1),ZAHYO(I+3)
endif
C READ(60,*)AN(I/3+1),ZAHYO(I),ZAHYO(I+1),ZAHYO(I+2)
210 CONTINUE
CLOSE(60,STATUS='KEEP')
RETURN
299 CLOSE(60,STATUS='KEEP')
STOP
END
cc
C 結果出力サブルーチン1
cc
INCLUDE'PARAMETER'
INTEGERN,J,N3
REAL*8TEMP
REAL*8ZAHYO(MAX_ATOM3)
CHARACTER*8 AN(MAX_ATOM)
REAL*8DIFF1(MAX_ATOM3)
REAL*8ENERGY
INTEGERI
c結果のファイルへの出力(result.xyz)
cresult.xyzは一回ずつの原子の座標をプロットする。アニメーションで確認できる。
OPEN(1,STATUS='OLD',FILE='result.xyz')
WRITE(1,*)N
C WRITE(1,*)'TIME=',T*TIME_DIV*1.0D15,'fs ENERGY=',ENERGY
c WRITE(1,*)'fs ENERGY=',ENERGY,'TEMP(K)=',TEMP
WRITE(1,304)'fs ENERGY(eV)=',ENERGY
N3=N*3
DO311I=1,N3,3
WRITE(1,305)AN(I/3+1),ZAHYO(I),ZAHYO(I+1),ZAHYO(I+2)
# ,DIFF1(I)*10.0,DIFF1(I+1)*10.0,DIFF1(I+2)*10.0
304 FORMAT(a20,f15.2)
305 FORMAT(a5,f15.7,f15.7,f15.7,f15.7,f15.7,f15.7)
311 CONTINUE
C
C結果出力サブルーチン2
C
c
c302 SUBROUTINEPRINT_ZAHYO2(N,ZAHYO,T,ENERGY,DIFF1,AN)
c INCLUDE'PARAMETER'
c INTEGERN,N3
c INTEGERT
c REAL*8ZAHYO(MAX_ATOM3)
c CHARACTER*8AN(MAX_ATOM)
c REAL*8DIFF1(MAX_ATOM3)
c REAL*8ENERGY
c INTEGERI
c
c結果のファイルへの出力(result1.xyz)
c OPEN(STATUS='OLD',UNIT=N,FILE='result1.xyz')
c WRITE(N,*)N
C WRITE(N,*)'TIME=',T*TIME_DIV*1.0D15,'fs ENERGY=',ENERGY
c WRITE(N,*)'fs ENERGY=',ENERGY,'TEMP(K)=',T
c
C
N3=N*3
DO312I=1,N3,3
WRITE(N,306)AN(I/3+1),ZAHYO(I),ZAHYO(I+1),ZAHYO(I+2)
# ,DIFF1(I)*10.0,DIFF1(I+1)*10.0,DIFF1(I+2)*10.0
306 FORMAT(a5,f12.7,f12.7,f12.7,f12.7,f12.7,f12.7)
312 CONTINUE
398 RETURN
c399 STOP
END
CC
C 最近接原子のLISTを作成するサブルーチン
CC
C
C LIST(0,N):原子番号Nの近接原子の数
C LIST(I,N):原子番号NのI番目の近接原子インデックス
C LIST(I+MAX_LIST_N,N):原子番号Nの I番目の近接原子番号
401 SUBROUTINEMAKE_LIST(N,ZAHYO,LIST,KYORI,CUTOFF,IDX)
INCLUDE'PARAMETER'
INTEGERN,N3
REAL*8 ZAHYO(MAX_ATOM3)
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
INTEGERIDX
REAL*8 MAX_X,MAX_Y,MAX_Z
REAL*8 MIN_X,MIN_Y,MIN_Z
INTEGERI,K,J
INTEGERE,F,G
INTEGERL,M
INTEGERMESH_X,MESH_Y,MESH_Z
INTEGERTX,TY,TZ
INTEGERTI,TJ,TI3,TJ3
INTEGERMESH(0:MAX_MESH_ATOM,
# 0:MAX_MESH_X+1,0:MAX_MESH_Y+1,0:MAX_MESH_Z+1)
REAL*8TMP_R,TMP_CUTOFF
REAL*8SIN
N3=N*3
C リスト配列の初期化
DO405I=1,MAX_ATOM
LIST(0,I)=0
405 CONTINUE
C MESH配列の初期化
DO406I=0,MAX_MESH_X+1
DO406J=0,MAX_MESH_Y+1
DO406K=0,MAX_MESH_Z+1
MESH(0,I,J,K)=0
406 CONTINUE
C 系の座標の最大値-最小値を検索する
MAX_X=ZAHYO(1)
MAX_Y=ZAHYO(2)
MAX_Z=ZAHYO(3)
MIN_X=ZAHYO(1)
MIN_Y=ZAHYO(2)
MIN_Z=ZAHYO(3)
DO410I=4,N3,3
IF(MAX_X.LT.ZAHYO(I))THEN
MAX_X=ZAHYO(I)
ENDIF
IF(MAX_Y.LT.ZAHYO(I+1))THEN
MAX_Y=ZAHYO(I+1)
ENDIF
IF(MAX_Z.LT.ZAHYO(I+2))THEN
MAX_Z=ZAHYO(I+2)
ENDIF
IF(MIN_X.GT.ZAHYO(I))THEN
MIN_X=ZAHYO(I)
ENDIF
IF(MIN_Y.GT.ZAHYO(I+1))THEN
MIN_Y=ZAHYO(I+1)
ENDIF
IF(MIN_Z.GT.ZAHYO(I+2))THEN
MIN_Z=ZAHYO(I+2)
ENDIF
410 CONTINUE
c 最大、最小の画面表示（確認）
WRITE(*,*)'MAX_',MAX_X,MAX_Y,MAX_Z
WRITE(*,*)'MIN_',MIN_X,MIN_Y,MIN_Z
WRITE(*,*)'D_.',REAL(MESH_D)
C 一度、大まかに原子をグルーピングする。
C 最大値最小値から、メッシュの個数を決める
MESH_X=1+((MAX_X-MIN_X)/REAL(MESH_D))
MESH_Y=1+((MAX_Y-MIN_Y)/REAL(MESH_D))
MESH_Z=1+((MAX_Z-MIN_Z)/REAL(MESH_D))
C MESH用の配列 に収まらない場合STOP
IF(MESH_X.LE.MAX_MESH_X)THEN
GOTO415
ENDIF
IF(MESH_Y.LE.MAX_MESH_Y)THEN
GOTO415
ENDIF
IF(MESH_Z.LE.MAX_MESH_Z)THEN
GOTO415
ENDIF
WRITE(*,*)'MESHisoverflow:MAX_MESH.', MESH_X,MESH_Y,MESH_Z
WRITE(*,*)'MESHisoverflow:MAX_.',MAX_X,MAX_Y,MAX_Z
WRITE(*,*)'MESHisoverflow:MIN_.',MIN_X,MIN_Y,MIN_Z
STOP
C 全ての原子をMESHに グルーピングする。
415 DO420I=1,N
I3=I*3-2
TX=1+((ZAHYO(I3)-MIN_X)/MESH_D)
TY=1+((ZAHYO(I3+1)-MIN_Y)/MESH_D)
C MESH配列があるれる場合STOP
IF(MESH(0,TX,TY,TZ).EQ.MAX_MESH_ATOM)THEN
WRITE(*,*)'MESH_ATOMisoverflow:MAX_MESH_ATOM.'
STOP
ENDIF
MESH(0,TX,TY,TZ)= MESH(0,TX,TY,TZ)+1
MESH(MESH(0,TX,TY,TZ),TX,TY,TZ)=I
420 CONTINUE
C MESHのブロックの 近接内(27ブロック)
C で 距離を計算し 近接LISTを作成する
C ついでに カットオフ関数の値も計算する
DO443I=1,MESH_X
DO442J=1,MESH_Y
DO441K=1,MESH_Z
IF(MESH(0,I,J,K).EQ.0)THEN
GOTO441
ENDIF
DO433 E=I-1 ,I+1
DO432F=J-1, J+1
DO431G=K-1 ,K+1
IF(MESH(0,E,F,G).EQ.0)THEN
GOTO431
ENDIF
DO425 L=1,MESH(0,I,J,K)
DO424M=1,MESH(0,E,F,G)
TI=MESH(L,I,J,K)
TJ=MESH(M,E,F,G)
IF(TI.LE.TJ)THEN
GOTO424
ENDIF
TI3=TI*3-2
TJ3=TJ*3-2
C
C 原子間距離の計算
C
TMP_R=
# ((ZAHYO(TJ3 )-ZAHYO(TI3 ))**2+
# (ZAHYO(TJ3+1)-ZAHYO(TI3+1))**2+
# (ZAHYO(TJ3+2)-ZAHYO(TI3+2))**2 )**0.5d0
C
C 近接でない場合場合 処理を飛ばす
C
IF(TMP_R.GT.(PRM_CR+PRM_CD) )THEN
GOTO424
ENDIF
C LISTの配列があふれる時STOP
IF((LIST(0,TI).EQ.MAX_LIST_N).OR.
# (LIST(0,TJ).EQ.MAX_LIST_N))THEN
WRITE(*,*)'LIST_Nisoverflow:MAX_LIST_N.'
STOP
ENDIF
IDX=IDX+1
KYORI(IDX)=TMP_R
LIST(0,TI)=LIST(0,TI)+1
LIST(LIST(0,TI),TI )=IDX
LIST(LIST(0,TI)+MAX_LIST_N,TI)=TJ
LIST(0,TJ)=LIST(0,TJ)+1
LIST(LIST(0,TJ),TJ )=IDX
LIST(LIST(0,TJ)+MAX_LIST_N,TJ)=TI
C
C カットオフ関数を計算
C
TMP_CUTOFF=1.0D0
IF(TMP_R.GE.(PRM_CR-PRM_CD))THEN
TMP_CUTOFF=0.5D0*
# (1.0D0-SIN(PI*(TMP_R-PRM_CR)/(2.0D0*PRM_CD)))
ENDIF
IF(TMP_R.GE.(PRM_CR+PRM_CD))THEN
TMP_CUTOFF=0.0D0
ENDIF
424 CONTINUE
425 CONTINUE
431 CONTINUE
432 CONTINUE
433 CONTINUE
441 CONTINUE
442 CONTINUE
443 CONTINUE
C DO452I=1,N
C write(*,*)"--",I
C DO451J=1,LIST(0,I)
C WRITE(*,*)LIST(J+MAX_LIST_N,I)
C451 CONTINUE
C452 CONTINUE
469 RETURN
c499 STOP
END
C
C i原子に関するTersoffポテンシャル
C
501 REAL*8FUNCTIONTERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
INCLUDE'PARAMETER'
INTEGERI
REAL*8 ZAHYO(MAX_ATOM3)
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
INTEGERJ,K,L,M
INTEGERI3,J3,K3
INTEGERIDXJ,IDXK
REAL*8G,GCOS
REAL*8ZETA
REAL*8I_X,I_Y,I_Z
REAL*8J_X,J_Y,J_Z
REAL*8K_X,K_Y,K_Z
REAL*8Bij
REAL*8EXP
REAL*8TMP
I3=I*3-2
I_X=ZAHYO(I3)
I_Y=ZAHYO(I3+1)
I_Z=ZAHYO(I3+2)
TMP=0.0D0
DO520L=1,LIST(0,I)
IDXJ=LIST(L,I)
J=LIST(L+MAX_LIST_N,I)
J3=J*3-2
J_X=ZAHYO(J3)
J_Y=ZAHYO(J3+1)
J_Z=ZAHYO(J3+2)
ZETA=0.0D0
DO510M=1,LIST(0,I)
IF(L.EQ.M)THEN
GOTO510
ENDIF
IDXK=LIST(M,I)
K=LIST(M+MAX_LIST_N,I)
K3=K*3-2
K_X=ZAHYO(K3)
K_Y=ZAHYO(K3+1)
K_Z=ZAHYO(K3+2)
GCOS=
# (
# (J_X-I_X)*(K_X-I_X)+
# (J_Y-I_Y)*(K_Y-I_Y)+
# (J_Z-I_Z)*(K_Z-I_Z)
# )/KYORI(IDXJ)/KYORI(IDXK)
IF(KYORI(IDXK).EQ.0.0D0)THEN
WRITE(*,*)KYORI(IDXK),I,J,K
G=(PRM_C*PRM_C)/(PRM_D*PRM_D)
-# (PRM_C*PRM_C)/
# (PRM_D*PRM_D+(PRM_H-GCOS)*(PRM_H-GCOS))
G=PRM_A*(G+1.0D0)
ZETA=ZETA+CUTOFF(IDXK)*G
510 CONTINUE
Bij=(1.0+(ZETA)**(PRM_ETA))**(-PRM_DELTA)
TMP=TMP+
# CUTOFF(IDXJ)*
# (
# PRM_EA*EXP(-PRM_LUMBDA1*KYORI(IDXJ))+
# (PRM_EB)*Bij*EXP(-PRM_LUMBDA2*KYORI(IDXJ))
# )
520 CONTINUE
TERSOFF_I=TMP
RETURN
c599 STOP
END
601 REAL*8FUNCTIONTERSOFF(N,ZAHYO,LIST,KYORI,CUTOFF)
INCLUDE'PARAMETER'
INTEGERI
INTEGERN
REAL*8 ZAHYO(MAX_ATOM3)
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
REAL*8 TERSOFF_I
TERSOFF=0.0D0
DO610I=1,N
TERSOFF=TERSOFF+
# TERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
610 CONTINUE
TERSOFF=TERSOFF/2.0D0
RETURN
c699 STOP
END
701 SUBROUTINECALC_DIFF1(N,ZAHYO,LIST,KYORI,CUTOFF,DIFF1)
INCLUDE'PARAMETER'
INTEGERN,N3
REAL*8 ZAHYO(MAX_ATOM3)
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
REAL*8 DIFF1(MAX_ATOM3)
INTEGERI,ID3
REAL*8 ZAHYO_ORG
REAL*8TERSOFF_LI
REAL*8TMP
N3=N*3
DO710I=1,N3
ID3=(I+2)/3
ZAHYO_ORG=ZAHYO(I)
ZAHYO(I)=ZAHYO_ORG+EPS1
CALLRECALC(ID3,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TERSOFF_LI(ID3,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I)=ZAHYO_ORG-EPS1
CALLRECALC(ID3,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
-# TERSOFF_LI(ID3,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I)=ZAHYO_ORG
CALLRECALC(ID3,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP/(2.0*EPS1)
DIFF1(I)=-TMP
710 CONTINUE
c799 STOP
END
801 SUBROUTINERECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
INCLUDE'PARAMETER'
INTEGERI,I3
REAL*8 ZAHYO(MAX_ATOM3)
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
INTEGERJ,J3
INTEGERIDXJ
INTEGERM
REAL*8TMP_R
REAL*8TMP_CUTOFF
REAL*8SIN
REAL*8SQRT
I3=I*3-2
DO810M=1,LIST(0,I)
IDXJ=LIST(M,I)
J=LIST(M+MAX_LIST_N,I)
IDXJ=LIST(M,I)
J3=J*3-2
TMP_R=
# SQRT ((ZAHYO(J3)-ZAHYO(I3))**2+
# (ZAHYO(J3+1)-ZAHYO(I3+1))**2+
# (ZAHYO(J3+2)-ZAHYO(I3+2))**2 )
KYORI(IDXJ)=TMP_R
TMP_CUTOFF=1.0D0
IF(TMP_R.GE.(PRM_CR-PRM_CD) )THEN
TMP_CUTOFF=0.5D0*
# (1.0D0-SIN(PI*(TMP_R-PRM_CR)/(2.0D0*PRM_CD)))
ENDIF
IF(TMP_R.GE.(PRM_CR+PRM_CD) )THEN
TMP_CUTOFF=0.0D0
ENDIF
CUTOFF(IDXJ)=TMP_CUTOFF
810 CONTINUE
RETURN
c899 STOP
END
901 REAL*8FUNCTIONTERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
INCLUDE'PARAMETER'
INTEGERI
REAL*8 ZAHYO(MAX_ATOM3)
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
REAL*8 TERSOFF_I
INTEGERL
INTEGERJ
TERSOFF_LI=TERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
DO910L=1,LIST(0,I)
J=LIST(L+MAX_LIST_N,I)
TERSOFF_LI=TERSOFF_LI+
# TERSOFF_I(J,ZAHYO,LIST,KYORI,CUTOFF)
910 CONTINUE
RETURN
c999 STOP
END
1001SUBROUTINEMAKE_LIST2(N,LIST,LIST2)
INCLUDE'PARAMETER'
INTEGERN
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
INTEGERLIST2(0:MAX_LIST2_N,MAX_ATOM)
INTEGERI,J,K,L,M
DO1010I=1,N
DO1005J=1,LIST(0,I)
L=LIST(J+MAX_LIST_N,I)
DO1003K=1,LIST(0,L)
M=LIST(K+MAX_LIST_N,L)
IF(M.EQ.I)THEN
GOTO1003
ENDIF
IF(M.EQ.L)THEN
GOTO1003
ENDIF
IF(LIST2(0,I).EQ.MAX_LIST2_N)THEN
WRITE(*,*)'MAX_LIST2_Nisoverflow'
&,i,LIST2(0,i),MAX_LIST2_N,MAX_ATOM
STOP
ENDIF
LIST2(0,I)=LIST2(0,I)+1
LIST2(LIST2(0,I),I)=M
1003 CONTINUE
1005 CONTINUE
1010CONTINUE
RETURN
c1099STOP
END
1101SUBROUTINECALC_DIFF2
#(N,ZAHYO,LIST,LIST2,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
INCLUDE'PARAMETER'
INTEGERN,N3
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
INTEGERLIST2(0:MAX_LIST2_N,MAX_ATOM)
REAL*8 ZAHYO(MAX_ATOM3)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
INTEGERI,J,K
INTEGERL
N3=N*3
C DO1110I=1,MAX_ATOM3
C DIFF2(0,I)=0
C1110CONTINUE
DO1150I=1,N
C 同一原子の座標系
CALLCALC_DIFF2_1
# (I,ZAHYO,LIST,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
C 隣合う原子の座標系
DO1145L=1,LIST(0,I)
J=LIST(L+MAX_LIST_N,I)
CALLCALC_DIFF2_2
# (I,J,ZAHYO,LIST,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
1145 CONTINUE
C 第二隣接の座標系
DO1146L=1,LIST2(0,I)
J=LIST2(L,I)
CALLCALC_DIFF2_3
# (I,J,ZAHYO,LIST,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
1146 CONTINUE
1150CONTINUE
RETURN
c1199STOP
END
1201SUBROUTINECALC_DIFF2_1
# (I,ZAHYO,LIST,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
INCLUDE'PARAMETER'
INTEGERI
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
REAL*8 ZAHYO(MAX_ATOM3)
REAL*8 KYORI(MAX_DATA_IDX)
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 TERSOFF_LI
INTEGERI3,J3
INTEGERII,JJ
INTEGERI3II,J3JJ
REAL*8X,Y
REAL*8TMP
REAL*8DATA
I3=I*3-3
J3=I3
DO1230II=1,3
DO1220JJ=1,3
I3II=I3+II
J3JJ=J3+JJ
IF(I3II.NE.J3JJ)THEN
X=ZAHYO(I3II)
Y=ZAHYO(J3JJ)
ZAHYO(I3II)=X+EPS2
ZAHYO(J3JJ)=Y+EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-EPS2
ZAHYO(J3JJ)=Y-EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP+
# TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X+EPS2
ZAHYO(J3JJ)=Y-EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
-# TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-EPS2
ZAHYO(J3JJ)=Y+EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
-# TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
DATA=TMP/(8.0D0*EPS2*EPS2)
ZAHYO(I3II)=X
ZAHYO(J3JJ)=Y
ELSE
X=ZAHYO(I3II)
ZAHYO(I3II)=X+2.0D0*EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-2.0D0*EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP+
# TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X+EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
-# TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
-# TERSOFF_LI(I,ZAHYO,LIST,KYORI,CUTOFF)
DATA=TMP/(3.0D0*EPS2*EPS2)
ZAHYO(I3II)=X
ENDIF
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLADD_DIFF2(DATA,I3II ,J3JJ,DIFF2,DIFF2_DATA)
1220 CONTINUE
1230CONTINUE
RETURN
c1299STOP
1301SUBROUTINEADD_DIFF2(DATA,I3II,J3JJ,DIFF2,DIFF2_DATA)
INCLUDE'PARAMETER'
INTEGERI3II,J3JJ
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DATA
IF(DIFF2(0,I3II).EQ.MAX_DIFF2_LIST)THEN
WRITE(*,*)'MAX_DIFF2_LISTisoverflow.'
STOP
ENDIF
DIFF2(0,I3II)=DIFF2(0,I3II)+1
DIFF2(DIFF2(0,I3II),I3II)=J3JJ
DIFF2_DATA(DIFF2(0,I3II),I3II)=DATA
C WRITE(*,*)I3II,J3JJ,DATA
RETURN
c1349STOP
END
1351SUBROUTINEADD2_DIFF2(DATA,I3II,J3JJ,DIFF2,DIFF2_DATA)
INCLUDE'PARAMETER'
INTEGERI3II,J3JJ
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DATA
INTEGERI
DO1360I=1,DIFF2(0,I3II)
IF(DIFF2(I,I3II).EQ.J3JJ)THEN
DIFF2_DATA(I,I3II)=DIFF2_DATA(I,I3II)+DATA
GOTO1380
ENDIF
1360CONTINUE
WRITE(*,*)'errorADD2_DIFF2'
STOP
C WRITE(*,*)I3II,J3JJ,DATA
1380RETURN
c1399STOP
END
1401SUBROUTINECALC_DIFF2_2
# (I,J,ZAHYO,LIST,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
INCLUDE'PARAMETER'
INTEGERI,J
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
INTEGERLIST_IJ(0:MAX_LIST_IJ)
REAL*8 ZAHYO(MAX_ATOM3)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 TERSOFF_I
REAL*8 TERSOFF_LIJ
INTEGERI3,J3
INTEGERII,JJ
INTEGERI3II,J3JJ
REAL*8X,Y
REAL*8TMP
REAL*8DATA
CALLMAKE_LIST_IJ(I,J,LIST_IJ,LIST)
I3=I*3-3
J3=J*3-3
DO1430II=1,3
DO1420JJ=1,3
I3II=I3+II
J3JJ=J3+JJ
X=ZAHYO(I3II)
Y=ZAHYO(J3JJ)
ZAHYO(I3II)=X+EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
# +TERSOFF_I(J,ZAHYO,LIST,KYORI,CUTOFF)
# +TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-EPS2
ZAHYO(J3JJ)=Y-EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
# +TERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
# +TERSOFF_I(J,ZAHYO,LIST,KYORI,CUTOFF)
# +TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X+EPS2
ZAHYO(J3JJ)=Y-EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
# -TERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
# -TERSOFF_I(J,ZAHYO,LIST,KYORI,CUTOFF)
# -TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-EPS2
ZAHYO(J3JJ)=Y+EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
# -TERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
# -TERSOFF_I(J,ZAHYO,LIST,KYORI,CUTOFF)
# -TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
DATA=TMP/(8.0D0*EPS2*EPS2)
ZAHYO(I3II)=X
ZAHYO(J3JJ)=Y
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
CALLADD_DIFF2(DATA,I3II,J3JJ,DIFF2,DIFF2_DATA)
1420 CONTINUE
1430CONTINUE
RETURN
c1499STOP
END
1501SUBROUTINEMAKE_LIST_IJ(I,J,LIST_IJ,LIST)
INCLUDE'PARAMETER'
INTEGERI,J
INTEGERLIST_IJ(0:MAX_LIST_IJ)
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
INTEGERL,M
LIST_IJ(0)=0
DO1540L=1,LIST(0,I)
DO1530M=1,LIST(0,J)
IF(LIST(L+MAX_LIST_N,I)
# .NE.LIST(M+MAX_LIST_N,J))THEN
GOTO1530
ENDIF
IF(LIST_IJ(0).EQ.MAX_LIST_IJ) THEN
WRITE(*,*)'MAX_LIST_IJisoverflow.'
ENDIF
LIST_IJ(0)=LIST_IJ(0)+1
LIST_IJ(LIST_IJ(0))=LIST(L+MAX_LIST_N,I)
1530 CONTINUE
1540CONTINUE
RETURN
c1599STOP
END
1601REAL*8FUNCTIONTERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
INCLUDE'PARAMETER'
INTEGERLIST_IJ(0:MAX_LIST_IJ)
REAL*8 ZAHYO(MAX_ATOM3)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
REAL*8 TERSOFF_I
INTEGERI,L
TERSOFF_LIJ =0.0D0
DO1610L=1,LIST_IJ(0)
I=LIST_IJ(L)
TERSOFF_LIJ=TERSOFF_LIJ+
# TERSOFF_I(I,ZAHYO,LIST,KYORI,CUTOFF)
1610 CONTINUE
RETURN
c1699 STOP
END
1701SUBROUTINECALC_DIFF2_3
# (I,J,ZAHYO,LIST,KYORI,CUTOFF,DIFF2,DIFF2_DATA)
INCLUDE'PARAMETER'
INTEGERI,J
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
INTEGERLIST_IJ(0:MAX_LIST_IJ)
REAL*8 ZAHYO(MAX_ATOM3)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 TERSOFF_LIJ
INTEGERI3,J3
INTEGERII,JJ
INTEGERI3II,J3JJ
REAL*8X,Y
REAL*8TMP
REAL*8DATA
CALLMAKE_LIST_IJ(I,J,LIST_IJ,LIST)
I3=I*3-3
J3=J*3-3
DO1730II=1,3
DO1720JJ=1,3
I3II=I3+II
J3JJ=J3+JJ
X=ZAHYO(I3II)
Y=ZAHYO(J3JJ)
ZAHYO(I3II)=X+EPS2
ZAHYO(J3JJ)=Y+EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-EPS2
ZAHYO(J3JJ)=Y-EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
# +TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X+EPS2
ZAHYO(J3JJ)=Y-EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
# -TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X-EPS2
ZAHYO(J3JJ)=Y+EPS2
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
TMP=TMP
# -TERSOFF_LIJ(LIST_IJ,ZAHYO,LIST,KYORI,CUTOFF)
ZAHYO(I3II)=X
ZAHYO(J3JJ)=Y
CALLRECALC(I,ZAHYO,LIST,KYORI,CUTOFF)
CALLRECALC(J,ZAHYO,LIST,KYORI,CUTOFF)
CALLADD_DIFF2(DATA,I3II ,J3JJ,DIFF2,DIFF2_DATA)
1720 CONTINUE
1730CONTINUE
RETURN
c1799STOP
END
1801SUBROUTINECONJGRAD(N,ZAHYO,DIFF1,DIFF2,DIFF2_DATA)
INCLUDE'PARAMETER'
INTEGERN
REAL*8ZAHYO(MAX_ATOM3)
REAL*8DIFF1(MAX_ATOM3)
INTEGERDIFF2(0:MAX_DIFF2_LIST,MAX_ATOM3)
REAL*8 DIFF2_DATA(MAX_DIFF2_LIST,MAX_ATOM3)
INTEGERN3
REAL*8P(MAX_ATOM3)
REAL*8R(MAX_ATOM3)
REAL*8AP(MAX_ATOM3)
REAL*8X(MAX_ATOM3)
REAL*8NORM_R2,NORM_R2_
REAL*8PAP
REAL*8A
REAL*8C
INTEGERI,J,K
INTEGERL
N3=N*3
NORM_R2=0.0D0
DO1810I=1,N3
X(I)=0.0D0
R(I)=DIFF1(I)
P(I)=DIFF1(I)
NORM_R2=NORM_R2+R(I)*R(I)
1810CONTINUE
C WRITE(*,*)0,NORM_R2
DO1860K=1,N3
IF(NORM_R2 .LT.1.0D-6)THEN
GOTO1865
ENDIF
DO1820I=1,N3
AP(I)=0.0D0
DO1815L=1,DIFF2(0,I)
J=DIFF2(L,I)
AP(I)=AP(I)+P(J)*DIFF2_DATA(L,I)
1815 CONTINUE
1820 CONTINUE
PAP=0.0D0
DO1830I=1,N3
PAP=PAP+P(I)*AP(I)
1830 CONTINUE
A=NORM_R2/PAP
C WRITE(*,*)'A=',A
NORM_R2_=0.0D0
DO1840I=1,N3
X(I)=X(I)+A*P(I)
R(I)=R(I)-A*AP(I)
NORM_R2_=NORM_R2_+R(I)*R(I)
1840 CONTINUE
C WRITE(*,*)K,NORM_R2_
C=NORM_R2_/NORM_R2
IF(C.GT.1.0D0)THEN
GOTO1865
NORM_R2=NORM_R2_
DO1850I=1,N3
P(I)=R(I)+C*P(I)
1850 CONTINUE
1860CONTINUE
1865DO1870I=1,N3
ZAHYO(I)=ZAHYO(I)+X(I)
1870CONTINUE
RETURN
c1899STOP
END
1901SUBROUTINEMAKE_LIST_VDW(N,ZAHYO,LIST_VDW,KYORI,CUTOFF,IDX
#,LIST,LIST2,AN)
INCLUDE'PARAMETER'
INTEGERN,N3
REAL*8 ZAHYO(MAX_ATOM3)
CHARACTER*8 AN(MAX_ATOM)
INTEGERLIST_VDW(0:MAX_LIST_VN2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
INTEGERIDX
INTEGERLIST(0:MAX_LIST_N2,MAX_ATOM)
INTEGERLIST2(0:MAX_LIST2_N,MAX_ATOM)
INTEGERI,K,J
INTEGERE,F,G
INTEGERL,M
INTEGERI3
INTEGERTX,TY,TZ
INTEGERTI,TJ,TI3,TJ3
REAL*8TMP_R,TMP_CUTOFF
C REAL*8SIN
C INTEGERIDX
N3=N*3
C リスト配列の初期化
DO1905I=1,MAX_ATOM
LIST_VDW(0,I)=0
1905 CONTINUE
C 系の座標の最大値-最小値を検索する
C 一度、大まかに原子をグルーピングする。
C 最大値最小値から、メッシュの個数を決める
C 全ての原子をMESHに グルーピングする。
DO1925TI=1,N
DO1924TJ=1,N
IF(TI.LE.TJ)THEN
GOTO1924
ENDIF
TI3=TI*3-2
TJ3=TJ*3-2
C
C 原子間距離の計算
C
TMP_R=
# ((ZAHYO(TJ3 )-ZAHYO(TI3 ))**2+
# (ZAHYO(TJ3+1)-ZAHYO(TI3+1))**2+
# (ZAHYO(TJ3+2)-ZAHYO(TI3+2))**2 )**0.5d0
C
C 範囲外のとき飛ばす。
IF(TMP_R.GT.(PRM_VCR2+PRM_VCD2))THEN
GOTO1924
ENDIF
C 同一分子内は 無視
IF(AN(TI).EQ. AN(TJ))THEN
ENDIF
CWRITE(*,*)AN(TI),AN(TJ)
C LISTの配列があふれる時STOP
IF((LIST_VDW(0,TI).EQ.MAX_LIST_VN).OR.
# (LIST_VDW(0,TJ).EQ.MAX_LIST_VN))THEN
WRITE(*,*)'LIST_VNisoverflow:MAX_LIST_VN.'
STOP
ENDIF
IDX=IDX+1
KYORI(IDX)=TMP_R
LIST_VDW(0,TI)=LIST_VDW(0,TI)+1
LIST_VDW(LIST_VDW(0,TI),TI )=IDX
LIST_VDW(LIST_VDW(0,TI)+MAX_LIST_VN,TI)=TJ
LIST_VDW(0,TJ)=LIST_VDW(0,TJ)+1
LIST_VDW(LIST_VDW(0,TJ),TJ )=IDX
LIST_VDW(LIST_VDW(0,TJ)+MAX_LIST_VN,TJ)=TI
C
C カットオフ関数を計算
C
TMP_CUTOFF=0.0D0
IF(TMP_R.GE.(PRM_VCR-PRM_VCD))THEN
TMP_CUTOFF=0.5D0*
# (1.0D0+SIN(PI*(TMP_R-PRM_VCR)/(2.0D0*PRM_VCD)))
ENDIF
IF(TMP_R.GE.(PRM_VCR+PRM_VCD))THEN
TMP_CUTOFF=1.0D0
ENDIF
IF(TMP_R.GE.(PRM_VCR2-PRM_VCD2))THEN
TMP_CUTOFF=0.5D0*
# (1.0D0-SIN(PI*(TMP_R-PRM_VCR2)/(2.0D0*PRM_VCD2)))
ENDIF
IF(TMP_R.GE.(PRM_VCR2+PRM_VCD2))THEN
TMP_CUTOFF=0.0D0
ENDIF
CUTOFF(IDX) =TMP_CUTOFF
1924 CONTINUE
1925 CONTINUE
C DO1952I=1,N
C write(*,*)"--",I
C DO1951J=1,LIST_VDW(0,I)
C WRITE(*,*)LIST_VDW(J+MAX_LIST_VN,I)
C1951 CONTINUE
C1952 CONTINUE
1969 RETURN
c1999 STOP
END
C
C原子Iに関したVanDerWaalsポテンシャル
C
2001REAL*8FUNCTIONVDW_I(I,ZAHYO,LIST_VDW,KYORI,CUTOFF)
INCLUDE'PARAMETER'
INTEGERI
REAL*8 ZAHYO(MAX_ATOM3)
INTEGERLIST_VDW(0:MAX_LIST_VN2,MAX_ATOM)
REAL*8 KYORI(MAX_DATA_IDX)
REAL*8 CUTOFF(MAX_DATA_IDX)
REAL*8 R,U
REAL*8VDW_FUNC
INTEGERL,IDXJ
U=0.0D0
DO2010L=1,LIST_VDW(0,I)
IDXJ=LIST_VDW(L,I)
R=KYORI(IDXJ)
U=U+CUTOFF(IDXJ)*VDW_FUNC(R)
2010CONTINUE
VDW_I=U
C write(*,*)"VDW",U,R