PROGRAM MXDTRICL C======================================================================== C## ## C## Program : MXDTRICL ## C## ## C## by Katsuyuki Kawamura (University of Tokyo) ## C## (Okayama University) ## C## (Hokkaido University) ## C## (Tokyo Institute of Technology) ## C## ## C## Configuration and Energy for Non-Cubic Systems ## C## (Oblique parallelepiped) ## C## with Pressure Control by stress tensor, ## C## and Quantum Correction for energy and pressure ## C## ## C## 2nd order interpolation from U and F tables ## C## ## C## First cubic version on Hitac 8800/8700 1980 ## C## First orthogonal (crystal) version 1983-Oct ## C## on CDC7600 at Manchester Univ. ## C## HITAC M-280/IAP version 1985-Sep-12 ## C## (Px, Py, Pz) pressure control version 1987-Feb-07 ## C## Pressure tensor and fractional coordinates 1987-Oct-29 ## C## Five element and input data format and history 1987-Nov-05 ## C## PC9800RA+NDP-FORTRAN-386 version 1989-Jan-26 ## C## Reviced for JCPE 1990-Apr-14 ## C## (XDORTO : DEFECT) 1990-Apr-21 ## C## 3-body interaction (H2O, Kumagai & Kats) 1991-Feb-02 ## C## Integrated version of MD and XD (MXD) 1991-May-22 ## C## Rearranged 1991-Oct-23 ## C## Seven comonents, rearranged 1992-Jan-23 ## C## Quatum corrections (Nakao & Kats) 1992-Mar-04 ## C## Ten comonents, rearranged 1992-Mar-31 ## C## Extended Andersen's pressure control 1992-Apr-07 ## C## (Katsuta & Kats) ## C## Metal (main group) potential 1992-Apr-18 ## C## Revised for JCPE version 1992-Aug-01 ## C## 2nd order interpolation from tables 1992-Sep-05 ## C## 2nd order interpolation of velocity 1992-Dec-12 ## C## Nose's thermostat 1992-Dec-14 ## C## Correction for trancation of VW-term 1993-Dec-10 ## C## Reviced 3-body term by Kuma 1994-Jan-30 ## C## L-J potential 1994-Jun-28 ## C## Nose's thermostat + quantum 1994-Sep-01 ## C## Charge - Dipole Interaction 1994-Sep-10 ## C## Improvement of Semi-classical MD 1995-Jun-15 ## C## FILE09.DAT format changed 1996-Jul-18 ## C## Model by Belonoshko & Dubrovinsky 1996-Sep-05 ## C## Shear 1997-Feb-18 ## C## Electric (N.SAWAGUCHI) & Gravity Field 1997-Jun-30 ## C## Apply constant strain rate 1997-Jun-30 ## C## Diatomic 3 chrge model 1997-Oct-10 ## C## 3-body j-i-k with j<>k 1999-Nov-16 ## C## 3-body sqrt(k1xk2) -> k1xk2 2000-May-01 ## C## POSISION-VELOCITY-ENERGY option 2000-Dec-16 ## C## Modify EWALD direct term 2001-Mar-24 ## C## 3-body j-i-k : modified 2001-Sep-11 ## c## File07.dat : format 2001-Dec-02 ## C## Polyatomic molecule 2002-Feb-23 ## C## Modify NETWOK analysis (c.n.=5) 2002-Sep-14 ## C## file07.dat (i10) and 3-body 2003-Jul-09 ## C## New multi-3-body 2003-Jul-28 ## ///// C## Separate file08.dat (file081.dat) 2005-Aug-11 ## c## file09v format change 2008-Nov-23 ## C## file08 format changed 2009-Feb-24 ## C## file09p and file09pv(pos) -> 5 figures 2009-Feb-25 ## C## Triatomic molecule (H2O, CO2, ...) 2010-May-26 ## C## Elwctric field at atoms 2010-Jul-09 ## C==============================================================================| C Format and parameters of 'FILE05.DAT' file | C------------------------------------------------------------------------------| C 1 MD.......I....:....I....:....I....:....I....:....I....:....I....:....I | C XD.......I... : | C 2 START :TITLE(60 CHARACTERS) : | C CONTINUE : (CONT.) : | C RESTART : : | C STOP : : | C 3 ECONOMY :IRECRD(1):IRECRD(2):IRECRD(3):IRECRD(4):IRECRD(5): : | C NORMAL : : : (50) : (M50,X5): (5) : : | C DETAIL : : : : : : : | C 4 NOACCUM :DTIME :FORMULA :(RCUTL) :(RCUTS) : : : | C ACCUM : : : : : : : | C 5 T NO-CNTL: : : : : : : | C T [BLANK]: : : : [No control on temperature]| C T SCALING:TMPGET :DELTMP :NTSTEP :(TDUMP) : : : | C T NOSE :TMPGET :DELTMP :STEMP : : : : | C 6 P NO-CNTL: : : : : : : | C P [BLANK]: : : : :[No control on pressure]| C P SCALING: SPRES(1):SPRES(2) :SPRES(3) :(PDUMP) : : : | C P ANDERSEN SPRES(1):SPRES(2) :SPRES(3) :VIRM(1) :VIRM(2) :VIRM(3) : | C P SHEAR : SPRES(1):SPRES(2) :SPRES(3) :VIRM(1) :VIRM(2) :VIRM(3) : | C : SPRES(4):SPRES(5) :SPRES(6) :VIRM(4) :VIRM(5) :VIRM(6) : | C : Pyz : Pzx : Pxy : : : : | C 7 V [BLANK]: : : :[Volume is changed with P-control]| C V CONST. : : : : [Volume is kept constant]| C V CELL :BOX(1) :BOX(2) :BOX(3) :BOX(4) :BOX(5) :BOX(6) : | C V DENSITY:DENSTY : : : : : : | C D CONST. :DENSTY : : : : : : | C 8 BUSING :MODE,MXN2: (ALPHA) : : : : : | C MORSE : : : : (Busing+Morse) : : | C MORSEQ : : : : : : : | C MORSE-PL : : : : (charge-dipole) : : | C MORSE-AT : : : : : : : | C BMH-EXP : : : : 3-body sqrt(k1xk2) : | C BMH-EXP* : : : : 3-body k1xk2 : | C BELONO : : : : (Belonoshko & Dubrovinsky) : | C TOSIFUMI : : : : : : : | C WOODCOCK : : : : : : : | C PAULING : : : : (Woodcock+Pauling f.) : | C METAL : : : : : : : | C STSUNE : : : : (Tsuneyuki et al.): : | C L-J : : : : : : : | C 81 N A NO. :ZI :WI :AI :BI :CI(VW) :DI() : | C - : : : : : [- not moved] | C * : : : : : [* dummy atoms]| C = : : : : : [= Morse only] | C 81e[BLANK] : : : : : : : | C 82 I J :DMIJ :BEIJ :RSIJ : Rswich : [Morse] : | C I J : D1ij : Be1ij : D2ij : Be2ij : Rswich : i3 : | C D3ij : Be3ij : r3ij : : [BMH-EXP]: | C J I J :FK3BP :ANG3BP :R3BLIM :R3BGD : [3-body] : | C J I K :FK3BP(1) :ANG3BP(1):R3BLIM(1):R3BGD(1) : [3-body(J<>K)]| C : : :R3BLIM(2):R3BGD(2) : : : | C 82e[BLANK] : : : : : : : | C : : : : : : : | C------------------------------------------------------------------------------| C 91 STRUCTURE: : : : 9:[Detail of final structure]| C 92 NETWORK :NFCION(1):NFCION(2): : 10:[Network structure analys.]| C 93 VELOCITY :NS09PV :PVMULT : : 11:[Record particle velocity]| C POSITION :NS09PV :PVMULT : : 11:[... ... position]| C ENERGY :NS09PV :PVMULT : : 11:[....... energy ]| C POSVELENE:NS09PV : : : 11:[..... pos,velo,ener]| C 94 QUANTUM : : : : 12:[Quantum correction]| C 95 PCF, RDF : ISTEP : Rend(A) : : 13:[Table of PCF and RCN]| C*96 DIPOLE : : : : 14:[E(dipole moment)]| C 97 CENTER : : : : 15:[Centring of atom cluster]| C 98 NO(MV=0) : : : : 16:[No correction for morment]| C 99 CRYSTAL : : : : 17:[MD of crystal structure]| C 9A BINARY : : : : 18:[Binary data for file09x.]| C 9B PRESSURE : NPRESS : : : 19:[Pressure tensor on file11]| C 9C ELEC.FIELD EFD1 : EFD2 : EFD3 : EFFEQ : 20:[Electric field]| C 9D GRAV.FIELD GFD1 : GFD2 : GFD3 : : 21: [Gravity field]| C 9E CONSTSHEAR VX-RY : VY-RZ : VX-RZ :(ps)-1 : 22:[Const.shear rat]| C 9F DIATOMIC : DINTRA :iatom2(1):iatom2(2): : 23:[Diatomic molec]| C 9g TRIATOMIC: Zmole31 : Dintra31:iatom3(1,1): (1,2): : icont : | C : Zmole32 : Dintra32:iatom3(2,1): (2,2): : : | C : : : : : 33:[Triatomic molecule]| C 9I MOLECULE :dMOLintra: Mstart : Mend : : 26:[Define molecule] | C 9L POLYATOM :dMOLintra:MOLstart : MOLend : :29:[Polyatomic molecule]| C 9n ........ : : : : : : : | C 9e [BLANK] : : : : : : : | C : : : : : : : | C MD.......I....:....I....:....I....:....I....:....I....:....I....:..: | C REPEAT 1 TO 9 : | C==============================================================================| C IRECRD NRECRD : C ----------------------------- ----------------------------- : C 1 Total number of steps Current step No. from 'START' : C 2 Interval of print PCF etc. Accumulation No. of PCF etc. : C (I2=N2 when 'ACCUM') : C 3 Interval of FILE07 recording Current step number : C (default: 50) in the current job : C 4 Interval of FILE09P recording Number of records in FILE09P : C (default: 50:MD. 5:XD) : C 5 Interval of FILE09V recording Number of records in FILE09V : C (default: 5) : C 6 Number of steps of current HIST Number OF HISTRY informations : C 7-8 Not used Not used : C 9 Interval of FILE09PV recording Number of steps in FILE09PV : C=======================================================================I C I/O number FLNAME Filename : C 5 - input from keyboad : C 15 ( 5) FILE05.DAT in : C 6, * - screen output out : C 16 ( 6) FILE06.DAT out : C 17 ( 7) FILE07.DAT in/out : C 18 ( 8) FILE08.DAT in/out : C 38 (18) FILE081.DAT in/out : C 19 ( 9) FILE09P.DAT in/out : C 10 (10) FILE10.DAT in : C 29 (11) FILE09V.DAT in/out : C 28 (12) FILE09PV.DAT out : C 27 (13) FILE11.DAT out : C 22 (19) TEMPO.DAT in/out(work) : C=======================================================================I C Variables in PARAMERER statement : C LNI : Maximum number of particles (ion or atom) in a basic cell : C LTB : Maximum table length of Coulomb energy and force : C LSR : Table length of short range interactions : C LEL : Maximum number of particle species : C LEE : Number of pairs of particle species : C LCT : Maximum number of steps : C LNV : Maxinum number of reciprocal lattice points in EWALD sum. : C LAA : Maximum number of atoms in a asymmetric unit (XD) : C LAT : Maximum number of atoms in a crystal unit cell (XD) : C=======================================================================I C P(3,LNI) : Fractional coordinates of atoms, 0=
>>>>', 5X,
* '< No. of steps >---< Temperature / K >---< Pressure ',
* '/ GPa >---< Date (yymmdd) >',6X,'I')
2002 FORMAT ('I',130('='),'I')
2221 FORMAT ('I ',I7,I5,I3,I7,5X, 99X, ' I')
2222 FORMAT ('I ',I7,I5,I3,I7,4X, I7,I5,I3,I7,4X, 74X, ' I')
2223 FORMAT ('I ',I7,I5,I3,I7,4X, I7,I5,I3,I7,4X, I7,I5,I3,I7,5X,
* 47X, ' I')
2224 FORMAT ('I ',I7,I5,I3,I7,4X, I7,I5,I3,I7,4X, I7,I5,I3,I7,4X,
* I7,I5,I3,I7, 26X, ' I')
2225 FORMAT ('I ',I7,I5,I3,I7,4X, I7,I5,I3,I7,4X, I7,I5,I3,I7,4X,
* I7,I5,I3,I7,4X, I7,I5,I3,I7,' I')
END
C
C
C ========
C================================================================ F07F08
SUBROUTINE F07F08 (INOEND)
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VALUES/ VAL(LVA), TVAL(LVA),TVALL(LVA),VALMAX(LVA),
* VAL0(LVA),SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* AVA(LVA,L50), NAV,NAVT
REAL *8 VAL,TVAL,TVALL,VALMAX,VAL0,SVAL,SVALL,VALMIN
COMMON /RADIAL/ NRDF(LTB,LEE),IPRDF(2)
INTEGER *4 NRDF
COMMON /GEOMET/ DTO(2),AVTHT(12),MBR(8,8,2),NRG(9,2),ITBR(121,12),
* RTO(2),SVTHT(12),NBR(8,8,2),MEB(9,2), NTT(121,12),
* NTO(2),NVTHT(12),ANGL(3,12),TTAB(LST), NTBL
COMMON /ACOORD/ NPT, BOXO(6),NIU(LAA),P0C(3,LAT),NSYM,ISYM(LNI),
* NPTP,NBOX(3),JON(LAT),PPC(3,LAT),MATM,
* RS(3,3,96),PPS(3,LAT),IHEX
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
C
COMMON /WORK01/ V10(3,LNI)
REAL *8 V10
COMMON /TIMDAT/ KKTIME(7,2)
C
CHARACTER *10 RUNO18,RUNO19
CHARACTER *4 TITLE0(15), BIN
CHARACTER *1 DEFECT
integer *4 iform7
INTEGER *4 IYEAR,IMONTH,IDAY, IHOUR,IMINUT,ISECND,I100TH
C
IF (INOEND.EQ.1) GO TO 501
C --------------------------------------------- Read from FILE07.DAT
C system description, coordinates and velocities
iform7 = 0
OPEN (17, FILE=FLNAME(7), STATUS='OLD',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
7 READ (17,7007) TITLE0, NJOB, BIN,
* NTION, NCOMPO, (NRECRD(I),I=1,9)
C
IF (NTION.GT.LNI) THEN
WRITE (*,*) 'Error: No. of ions (', NTION, ') is too large',
* ' (LNI=', LNI, ') !!!'
STOP
END IF
IF (NCOMPO.GT.LEL) THEN
WRITE (*,*) 'Error: No. of ion species (',NCOMPO,') is ',
* 'too large (LEL=',LEL,') !!!'
STOP
END IF
RUNOPT(18) = ' '
IF (BIN.EQ.'BIN ') RUNOPT(18) = 'BINARY '
C
READ (17,7017) (ATOM(I),I=1,NCOMPO)
READ (17,7018) (NION(I),I=1,NCOMPO)
READ (17,7018) (IONS(1,I),I=1,NCOMPO)
READ (17,7018) (IONS(2,I),I=1,NCOMPO)
READ (17,7070) TEMP, DELTMP,TMPGET, (SPRES(I),I=1,3),
* DTIME, RUNOPT(51), BOX,
* DENSTY, RUNOPT(52), VBOX
IF (RUNOPT(51).EQ.'THERMOSTAT') READ (17,7080) STEMP, VSTEMP
IF (RUNOPT(52).EQ.'H-TENSOR ') THEN
DO 100 I = 1, 3
READ (17,7080) (H(I,J),J=1,3)
100 CONTINUE
END IF
C
if (iform7.eq.0) then
WRITE (*,1177) TITLE0, TITLE
1177 FORMAT (5X,14('='),' Titles in FILE07.DAT and FILE05.DAT are ',
* 14('=') / '=====[F7]: ',15A4,' ===== ' /
* '=====[F5]: ',15A4,' ===== ' )
end if
C
IF (NTION.GT.LNI) WRITE (*,*) 'The number of atoms :',NTION,
* ' is greater than LNI:',LNI
NTIOND = 0
DO 110 I = 1, NTION
IOND(I) = 1
if (iform7.eq.0 ) then
READ (17,7700,err=7878) (P(J,I),J=1,3),
* DEFECT, (V10(J,I),J=1,3), (P0(J,I),J=1,3)
else
READ (17,7702,err=7878) (P(J,I),J=1,3),
* DEFECT, (V10(J,I),J=1,3), (P0(J,I),J=1,3)
end if
if (abs(V10(1,i)-5.0)+abs(V10(2,i)-5.0)+
* abs(V10(3,i)-5.0) .gt. 3.0 ) then
if (iform7.eq.1) then
write (6,*) i,'-th atom is strange'
stop
end if
iform7 = 1
rewind 17
go to 7
end if
IF (DEFECT.NE.' ') THEN
write (6,*) i,defect
IOND(I) = 0
NTIOND = NTIOND + 1
V10(1,I) = 0.0D0
V10(2,I) = 0.0D0
V10(3,I) = 0.0D0
END IF
DO 105 J = 1, 3
V(J,I) = (V10(J,I)-5.0D0) * 0.1D0
105 CONTINUE
110 CONTINUE
IF (NTIOND.GT.0) WRITE (*,7979) NTIOND
7979 FORMAT (1X,I6,' DEFECTS WERE DETECTED ')
IF (NRECRD(6).GT.0) THEN
READ (17,7800,END=180,ERR=180) ((IHISTR(J,I),J=1,4),
* I=1,NRECRD(6))
GO TO 190
180 NRECRD(6) = 0
190 END IF
IRECRD(6) = 0
CLOSE (17)
if (iform7.eq.0) write (6,*) 'Format of file07.dat will be ',
* 'converted.'
go to 201
c
7878 write (6,*) 'File07.dat : error at the line ',i+9
stop
C
201 IF (RUNOPT(2).EQ.'RESTART ') THEN
RUNOPT(2) = 'START '
NRECRD(6) = 0
DO 210 I = 1,NTION
DO 210 J = 1, 3
P(J,I) = P0(J,I)
210 CONTINUE
END IF
C
NBOX(1) = 1
NBOX(2) = 1
NBOX(3) = 1
IF (RUNOPT(17).EQ.'CRYSTAL ') CALL FILE10
C
IF (TITLE(1).NE.'BENC' .OR.
* TITLE(2).NE. 'HMAR' ) THEN
C file09p.dat : COORDINATES AT EACH 5 STEP
OPEN (19, FILE=FLNAME(9), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
C file09v.dat : VALUES AT EACH 5 STEP
OPEN (29, FILE=FLNAME(11), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
END IF
C
IF (RUNOPT(2).EQ.'CONTINUE '.OR.RUNOPT(2).EQ.'CONTINUE ') THEN
NJOB(2) = NJOB(2) + 1
C ----------------------------------- Read from FILE08.DAT
C PCF, properties, etc.
OPEN (18, FILE=FLNAME(8), STATUS='OLD',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
REWIND 18
READ (18,8001) NCUT0,NRCUT(1),NRECRD(2),NAV,NAVT,NTBL,
* MXCUT,NPAIR
DO 301 J = 1, LEE
DO 301 N = 1, LTB
NRDF(N,J) = 0
301 CONTINUE
DO 311 I = NCUT0, NRCUT(1)
READ (18,8001) (NRDF(I,J),J=1,NPAIR)
311 CONTINUE
DO 321 I = 1, LVA
READ (18,8003) TVAL(I),SVAL(I),SVALL(I),VAL0(I)
321 CONTINUE
c DO 331 I = 1, NAV
c READ (18,8003) (AVA(J,I),J=1,LVA)
c 331 CONTINUE
READ (18,8003) (AU(I),I=1,NTION)
DO 341 I = 1, 12
READ (18,8003) (ANGL(J,I),J=1,3)
341 CONTINUE
DO 351 K = 1, 2
DO 351 J = 1, 8
READ (18,8001) (MBR(I,J,K),I=1,8)
351 CONTINUE
DO 361 J = 1, 2
READ (18,8001) (NRG(I,J),I=1,9)
361 CONTINUE
DO 371 I = 1, 121
READ (18,8005) (ITBR(I,J),J=1,12)
371 CONTINUE
IF (RUNOPT(17).EQ.'CRYSTAL ') THEN
READ (18,8004) ((PPC(J,N),J=1,3),
* (PPS(J,N),J=1,3),N=1,NPT)
END IF
CLOSE (18)
c
OPEN (38, FILE=FLNAME(18), STATUS='OLD',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
REWIND 38
DO 331 I = 1, NAV
READ (38,8003) (AVA(J,I),J=1,LVA)
331 CONTINUE
close (38)
C
CALL FILE09
ELSE
NJOB(1) = NJOB(1) + 1
NJOB(2) = 1
NRECRD(4) = 0
NRECRD(5) = 0
IF (TITLE(1).NE.'BENC' .OR.
* TITLE(2).NE. 'HMAR' ) THEN
REWIND 29
REWIND 19
END IF
END IF
RETURN
C
C ========================================= Output file07 and file08
501 NRECRD(6) = NRECRD(6) + 1
CALL KCLOCK (IYEAR,IMONTH,IDAY,IHOUR,IMINUT,ISECND,I100TH)
IHISTR(1,NRECRD(6)) = IRECRD(6)
IHISTR(2,NRECRD(6)) = INT(TMPGET)
IHISTR(3,NRECRD(6)) = INT((SPRES(1)+SPRES(2)+SPRES(3))/3.0)
IHISTR(4,NRECRD(6)) = IYEAR*10000 + IMONTH*100 + IDAY
IRECRD(6) = 0
IF (NRECRD(6).GT.1) THEN
KHIST = NRECRD(6) - 1
IF (IHISTR(2,NRECRD(6)).EQ.IHISTR(2,KHIST).AND.
* IHISTR(3,NRECRD(6)).EQ.IHISTR(3,KHIST)) THEN
IHISTR(1,KHIST)=IHISTR(1,NRECRD(6))+IHISTR(1,KHIST)
IHISTR(4,KHIST)=IHISTR(4,NRECRD(6))
NRECRD(6) = KHIST
END IF
END IF
IF (TITLE(1).EQ.'BENC' .AND.
* TITLE(2).EQ. 'HMAR' ) GO TO 699
C
RUNO18 = ' '
RUNO19 = 'H-TENSOR '
IF (RUNOPT(5).EQ.'T NOSE ') RUNO18 = 'THERMOSTAT'
C
C ---------------------------------------------- Write on FILE07.DAT
C system description, coordinates and velocities
C
OPEN (17, FILE=FLNAME(7), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
REWIND 17
BIN = ' '
IF (RUNOPT(18).EQ.'BINARY ') BIN = 'BIN '
WRITE (17,7007) TITLE, NJOB, BIN,
* NTION, NCOMPO, (NRECRD(I),I=1,9)
WRITE (17,7017) (ATOM(I),I=1,NCOMPO)
WRITE (17,7018) (NION(I),I=1,NCOMPO)
WRITE (17,7018) (IONS(1,I),I=1,NCOMPO)
WRITE (17,7018) (IONS(2,I),I=1,NCOMPO)
WRITE (17,7070) TEMP, DELTMP,TMPGET, (SPRES(I),I=1,3),
* DTIME, RUNO18, BOX,
* DENSTY, RUNO19, VBOX
IF (RUNO18.EQ.'THERMOSTAT') WRITE (17,7080) STEMP,VSTEMP
DO 503 I = 1, 3
WRITE (17,7080) (H(I,J),J=1,3)
503 CONTINUE
do 508 io = 1, ncompo
DO 507 I = ions(1,io), ions(2,io)
DO 505 J = 1, 3
V10(J,I) = V(J,I) * 10.0D0 + 5.0D0
505 CONTINUE
DEFECT = ' '
IF (IOND(I).EQ.0) DEFECT = '*'
WRITE (17,7702) (P(J,I),J=1,3),DEFECT,(V10(J,I),J=1,3),
* (P0(J,I),J=1,3), io
507 CONTINUE
508 continue
WRITE (17,7800) ((IHISTR(J,I),J=1,4),I=1,NRECRD(6))
ENDFILE (17)
REWIND 17
CLOSE (17)
C
C -------------------------------------------- Write on FILE08.DAT
C PCF, properties, etc.
DO 512 N = 1, NRCUT(1)
DO 511 J = 1, LEE
IF (NRDF(N,J).GT.0) GO TO 513
511 CONTINUE
512 CONTINUE
513 NCUT0 = N - 1
NPAIR = NCOMPO * (NCOMPO+1) / 2
OPEN (18, FILE=FLNAME(8), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
REWIND 18
WRITE (18,8001) NCUT0,NRCUT(1),NRECRD(2),NAV,NAVT,NTBL,MXCUT,
* NPAIR
DO 611 I = NCUT0, NRCUT(1)
WRITE (18,8001) (NRDF(I,J),J=1,NPAIR)
611 CONTINUE
DO 621 I = 1, LVA
WRITE (18,8003) TVAL(I),SVAL(I),SVALL(I),VAL0(I)
621 CONTINUE
c DO 631 I = 1, NAV
c WRITE (18,8003) (AVA(J,I),J=1,LVA)
c 631 CONTINUE
WRITE (18,8003) (AU(I),I=1,NTION)
DO 641 I = 1, 12
WRITE (18,8003) (ANGL(J,I),J=1,3)
641 CONTINUE
DO 651 K = 1, 2
DO 651 J = 1, 8
WRITE (18,8001) (MBR(I,J,K),I=1,8)
651 CONTINUE
DO 661 J = 1, 2
WRITE (18,8001) (NRG(I,J),I=1,9)
661 CONTINUE
DO 671 J = 1, 121
WRITE (18,8005) (ITBR(J,I),I=1,12)
671 CONTINUE
IF (RUNOPT(17).EQ.'CRYSTAL ') THEN
WRITE (18,8004) ((PPC(J,N),J=1,3),
* (PPS(J,N),J=1,3),N=1,NPT)
END IF
C
ENDFILE (18)
REWIND 18
CLOSE (18)
c
OPEN (38, FILE=FLNAME(18), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
REWIND 38
DO 631 I = 1, NAV
WRITE (38,8003) (AVA(J,I),J=1,LVA)
631 CONTINUE
ENDFILE (38)
REWIND 38
CLOSE (38)
C
699 WRITE (*,4001) IRECRD(1)
4001 FORMAT (15('='),' Files were updated ',13('='),
* ' End=',I6,2X,15('='))
WRITE (*,1178) TITLE
1178 FORMAT ('<<<===== ',15A4,' ====>>>')
RETURN
C
C -------------------------------------------- Formats of file07.dat
7007 FORMAT (15A4,2I5, 1X,A4 / I7,I3, 9I10)
7017 FORMAT (10(2X,A4) )
7018 FORMAT (10I6 )
7070 FORMAT (F10.2,F10.4,F10.2, 3F10.5 /
* E10.3, A10, 6F10.6 /
* F10.6, A10, 6F10.6 )
7080 FORMAT (10X,3F20.10)
7700 FORMAT (3F9.7, A1, 3F8.6, 1X, 3F9.6)
7701 FORMAT (3F9.7, A1, 3F8.6, 1X, 3F9.6, 1x,i2)
7702 FORMAT (3F10.8, A1, 3F9.7, 1X, 3F10.6, 1x,i2)
7800 FORMAT (3(I10,I5,I4,1X,I6))
C -------------------------------------------- Formats of file08.dat
8001 FORMAT (10I10)
8003 FORMAT (1P5E16.9)
8004 FORMAT (0P3F12.6,4X,3F12.6)
8005 FORMAT (12I8)
END
C
C
C ========
C================================================================ FILE09
SUBROUTINE FILE09
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VALUES/ VAL(LVA), TVAL(LVA),TVALL(LVA),VALMAX(LVA),
* VAL0(LVA),SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* AVA(LVA,L50), NAV,NAVT
REAL *8 VAL,TVAL,TVALL,VALMAX,VAL0,SVAL,SVALL,VALMIN
COMMON /ACOORD/ NPT, BOXO(6),NIU(LAA),P0C(3,LAT),NSYM,ISYM(LNI),
* NPTP,NBOX(3),JON(LAT),PPC(3,LAT),MATM,
* RS(3,3,96),PPS(3,LAT),IHEX
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
C
COMMON /WORK02/ IP(3,LNI), PP(3,LNI)
C
REAL *8 HH(3,3), VALVAL(LVA)
C
IF (TITLE(1).EQ.'BENC' .AND.
* TITLE(2).EQ. 'HMAR' ) RETURN
C --------------------------------------- Work file for continuation
OPEN (22, FILE = FLNAME(19), STATUS = 'UNKNOWN',
* ACCESS = 'SEQUENTIAL', FORM = 'FORMATTED' )
C
C -------------------------------------------- FILE09V.DAT
1991 FORMAT (F10.3,7F10.5 / 8F10.3 /
* F10.6, F10.4, 3F10.6,3F10.6 / 10F9.3 / 10F9.3 )
REWIND 29
REWIND 22
DO 410 K = 1, NRECRD(5)
READ (29,1991) (VALVAL(I),I=1,LVA)
WRITE (22,1991) (VALVAL(I),I=1,LVA)
410 CONTINUE
ENDFILE 22
REWIND 29
REWIND 22
DO 420 K = 1, NRECRD(5)
READ (22,1991) VALVAL
WRITE (29,1991) VALVAL
420 CONTINUE
C
C -------------------------------------------------- FILE09P.DAT
IF (RUNOPT(18).EQ.'BINARY ') THEN
CLOSE (22)
OPEN (22, FILE = FLNAME(19), STATUS = 'UNKNOWN',
* ACCESS = 'SEQUENTIAL', FORM = 'UNFORMATTED' )
END IF
MMMMM = NTION
IF (RUNOPT(17).EQ.'CRYSTAL ') MMMMM = NPTP
REWIND 19
REWIND 22
IF (RUNOPT(18).EQ.'BINARY ') THEN
DO 440 K = 1, NRECRD(4)
READ (19) L, HH
READ (19) ((PP(J,I),J=1,3),I=1,MMMMM)
WRITE (22) L, HH
WRITE (22) ((PP(J,I),J=1,3),I=1,MMMMM)
440 CONTINUE
REWIND 19
REWIND 22
DO 450 K = 1, NRECRD(4)
READ (22) L, HH
READ (22) ((PP(J,I),J=1,3),I=1,MMMMM)
WRITE (19) L, HH
WRITE (19) ((PP(J,I),J=1,3),I=1,MMMMM)
450 CONTINUE
ELSE
DO 460 K = 1, NRECRD(4)
READ (19,9002) L, HH
READ (19,9001) ((IP(J,I),J=1,3),I=1,MMMMM)
WRITE (22,9002) L, HH
WRITE (22,9001) ((IP(J,I),J=1,3),I=1,MMMMM)
460 CONTINUE
REWIND 19
REWIND 22
DO 470 K = 1, NRECRD(4)
READ (22,9002) L, HH
READ (22,9001) ((IP(J,I),J=1,3),I=1,MMMMM)
WRITE (19,9002) L, HH
WRITE (19,9001) ((IP(J,I),J=1,3),I=1,MMMMM)
470 CONTINUE
END IF
C
CLOSE (22)
RETURN
C ----------------------------------------- Formats of file09a.dat's
9001 FORMAT (18I5)
9002 FORMAT (I5,5X, 9F7.3)
END
C
C
C ========
C================================================================ FILE10
SUBROUTINE FILE10
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /ACOORD/ NPT, BOXO(6),NIU(LAA),P0C(3,LAT),NSYM,ISYM(LNI),
* NPTP,NBOX(3),JON(LAT),PPC(3,LAT),MATM,
* RS(3,3,96),PPS(3,LAT),IHEX
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
C
CHARACTER *4 HEX
C
C ------------------------------ Input file of xtal geometry
OPEN (10,FILE=FLNAME(10),STATUS='OLD',
* ACCESS='SEQUENTIAL',FORM='FORMATTED')
REWIND 10
READ (10,5010) BOXO,
* NBOX,NPT,NPTP,NSYM,HEX,MATM
READ (10,5012) (ATMXTL(J),J=1,MATM)
READ (10,5014) (NIU(J),J=1,MATM)
READ (10,5020) (JON(N),(P0C(J,N),J=1,3),N=1,NPTP)
READ (10,5030) (((RS(J,I,N),J=1,3),I=1,3),N=1,NSYM)
READ (10,5040) (ISYM(N),N=1,NTION)
REWIND 10
CLOSE (10)
IHEX = 0
IF (HEX.EQ.'HEX ') IHEX = 1
RETURN
5010 FORMAT (3F10.7,3F10.8 / 6I5,5X,A4,I6 )
5012 FORMAT ( 18A4 )
5014 FORMAT ( 18I4 )
5020 FORMAT (I5,3F10.7)
5030 FORMAT (9F6.1)
5040 FORMAT (12I6)
END
C
C
C ========
C================================================================ INITIA
SUBROUTINE INITIA (INOEND)
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C -------------------------------------------- Initial reading, etc.
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /VECTOR/ FNV(LNV), UNV(LNV), PNV(3,3,LNV),
* VEC(3,LNV), ZIA(LEM), ALPHA, UCSELF,UCSeLFI(LEM),
* MODE, NVN, NVEC(3,LNV)
REAL *8 FNV,UNV,PNV,VEC,ZIA,ALPHA,UCSELF,ucselfi
COMMON /VALUES/ VAL(LVA), TVAL(LVA),TVALL(LVA),VALMAX(LVA),
* VAL0(LVA),SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* AVA(LVA,L50), NAV,NAVT
REAL *8 VAL,TVAL,TVALL,VALMAX,VAL0,SVAL,SVALL,VALMIN
COMMON /RADIAL/ NRDF(LTB,LEE),IPRDF(2)
INTEGER *4 NRDF
COMMON /GEOMET/ DTO(2),AVTHT(12),MBR(8,8,2),NRG(9,2),ITBR(121,12),
* RTO(2),SVTHT(12),NBR(8,8,2),MEB(9,2), NTT(121,12),
* NTO(2),NVTHT(12),ANGL(3,12),TTAB(LST), NTBL
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
common /STRCTU/ lentab
COMMON /OUTERF/ EFD(3),EFREQ, GFD(3), STRT(3), MEFD
REAL *8 EFD, EFREQ, GFD, STRT
COMMON /MOLECU/ ZMOLE2(2), DMOLE2(4,LNI), DINTRA2(2),
* ZMOLE3(2), DMOLE3(4,LNI), DINTRA3(2),
* dMOLintra,
* NDMOLE2, IDMOLE2(3,LNI), IATOM2(2), MOLstart,
* ndmole3, idmole3(4,lni), iatom3(2,2),
* NMOLE, IMOLE(38,LNI), MMOLE(LNI), MOLend
real *8 zmole2,zmole3,dmole2,dmole3,dintra2,dintra3
COMMON /WORK01/ VV(3,LNI), DUM(3,LNI)
COMMON /WORK02/ IPV(3,LNI),IDUMMY(3,LNI)
C
REAL *8 BOXA(6), FA(3), param1,param2,param3,param4,param5
CHARACTER *4 AAX, ATY, THS1,THS2, RUNOP1
CHARACTER *10 RUNRUN, DUMMY
ATMNET(1) = ' '
ATMNET(2) = ' '
DO 10 I = 1, 53
RUNOPT(I) = ' '
10 CONTINUE
NRECRD(9) = 0
dMOLintra = 0.0
MOLstart = 0
MOLend = 0
do i=1, 2
dintra2(i) = 0.0
iatom2(i) = 0
zmole2(i) = 0.0
end do
C
C --------------------------------------- Data input from FILE05.DAT
C
IP0 = 0
INOEND = 0
30 READ (15,1001,END=888) RUNOPT(1)
RUNOP1 = RUNOPT(1)
IF (RUNOP1.EQ.'MDX.') THEN
RUNOPT(1) = 'MD.......:'
RUNOP1 = 'MD..'
IP0 = 1
END IF
IF (RUNOP1.EQ.'MD..') THEN
RUNOPT(1) = 'MD........'
RUNOPT(17) = 'AMORPHOUS '
END IF
IF (RUNOP1.EQ.'XD..') THEN
RUNOPT(1) = 'XD........'
RUNOPT(17) = 'CRYSTAL '
END IF
IF (RUNOP1.NE.'MD..' .AND.
* RUNOP1.NE.'XD..' ) GO TO 30
READ (15,1001,END=888) RUNOPT(2),TITLE
IF (RUNOPT(2).EQ.' ' .OR.
* RUNOPT(2).EQ.'STOP ' .OR.
* RUNOPT(2).EQ.'END ' ) GO TO 888
IF (RUNOPT(2).EQ.'CONT. ')
* RUNOPT(2) = 'CONTINUE '
GO TO 50
C
888 INOEND = -1
RETURN
C
C -------------------------------- Read file07.dat, file08.dat, etc.
50 CALL F07F08 (INOEND)
C -------------------------------------- Input file of xtal geometry
CALL TITLET (1,0)
C ------------------------------------------- Economy, normal detail
READ (15,1000) RUNOPT(3), AREC1, AREC2, AREC3, AREC4, AREC5
IRECRD(1) = INT(AREC1)
IRECRD(2) = INT(AREC2)
IRECRD(3) = INT(AREC3)
IRECRD(4) = INT(AREC4)
IRECRD(5) = INT(AREC5)
IF (IRECRD(1).GT.LCT) THEN
WRITE (6,*) 'The number of steps:',IRECRD(1),
* 'is too large (LCT=',LCT,')'
WRITE (6,*) 'Please chage all the LCT parameters'
STOP
END IF
IF (IRECRD(1).LT.IRECRD(2)) IRECRD(2) = IRECRD(1)
IF (MOD(IRECRD(1),IRECRD(2)).NE.0) IRECRD(2) = IRECRD(1)
IF (IRECRD(3).LE.0) IRECRD(3) = 50
IF (IRECRD(2).LT.IRECRD(3)) IRECRD(3) = IRECRD(2)
IF (IRECRD(4).LE.0) THEN
IF (RUNOP1.EQ.'MD..') IRECRD(4) = IRECRD(3)
IF (RUNOP1.EQ.'XD..') IRECRD(4) = 5
END IF
IF (IRECRD(5).LE.0) IRECRD(5) = 5
C ------------------------------------------------- Accume, noaccume
READ (15,1000) RUNOPT(4), DDT, FORMUL, RCUT(1), RCUT(2)
C ------------------------------------------------------ Temperatute
READ (15,1000) RUNRUN, TARGT, DELT, STEMP0, TDUMP
IF (RUNRUN.EQ.'T ') RUNOPT(5) = 'T NO-CNTL '
IF (RUNRUN.EQ.'T NO ') RUNOPT(5) = 'T NO-CNTL '
IF (RUNRUN.EQ.'T NO-CNTL ') RUNOPT(5) = 'T NO-CNTL '
IF (RUNRUN.EQ.'T SCALING ') THEN
RUNOPT(5) = 'T SCALING '
NTSTEP = STEMP0
END IF
IF (RUNRUN.EQ.'T NOSE ') RUNOPT(5) = 'T NOSE '
C --------------------------------------
IF (RUNOPT(5) .NE.'T NOSE ' .OR.
* RUNOPT(2) .NE.'CONTINUE ' .OR.
* RUNOPT(51).NE.'THERMOSTAT' ) THEN
STEMP = STEMP0
VSTEMP = 0.0
END IF
IF (NTSTEP.LE.0) NTSTEP = 1
DELTMP = DELT
TMPGET = TARGT
IF (TDUMP.LT.0.001) TDUMP = 0.5
C --------------------------------------------------------- Pressure
READ (15,1000) RUNRUN, (SPRES(I),I=1,3), VIRM(1),VIRM(2),VIRM(3)
pdump = 1.0
IF (RUNRUN.EQ.'P ') RUNOPT(6) = 'P NO-CNTL '
IF (RUNRUN.EQ.'P NO ') RUNOPT(6) = 'P NO-CNTL '
IF (RUNRUN.EQ.'P NO-CNTL ') RUNOPT(6) = 'P NO-CNTL '
IF (RUNRUN.EQ.'P SCALING ') THEN
RUNOPT(6) = 'P SCALING '
SPRES(4) = 0.0
SPRES(5) = 0.0
SPRES(6) = 0.0
pdump = virm(1)
if (pdump.lt.0.01) pdump = 1.0
end if
IF (RUNRUN.EQ.'P ANDERSEN') THEN
RUNOPT(6) = 'P ANDERSEN'
IF (ABS(VBOX(2)).LT.1.0E-9.AND.
* ABS(VBOX(3)).LT.1.0E-9 ) THEN
VBOX(1) = 0.0
VBOX(2) = 0.0
VBOX(3) = 0.0
END IF
END IF
C --------------------------------------------
IF (RUNOPT(6).NE.'P ANDERSEN'.AND.
* ABS(VBOX(2)).GT.1.0E-9.AND.
* ABS(VBOX(3)).GT.1.0E-9 ) THEN
VBOX(1) = 0.0
VBOX(2) = 0.0
VBOX(3) = 0.0
END IF
IF (RUNRUN.EQ.'P SHEAR ') THEN
RUNOPT(6) = 'P SHEAR '
READ (15,1000) DUMMY, (SPRES(I),I=4,6),
* (VIRM(I),I=4,6)
END IF
C ----------------------------------------------------------- Volume
READ (15,1000) RUNRUN, BOXA
IF (RUNRUN.EQ.' ') RUNOPT(7) = 'V FREE '
IF (RUNRUN.EQ.'V CONST. ') RUNOPT(7) = 'V CONST. '
IF (RUNRUN.EQ.'V CONTROL ') RUNOPT(7) = 'V CONST. '
IF (RUNRUN.EQ.'D CONST. ') RUNOPT(7) = 'D CONST. '
IF (RUNRUN.EQ.'D CONTROL ') RUNOPT(7) = 'D CONST. '
C --------------------------------------- Change cell size
IF (RUNRUN.EQ.'V CELL ') THEN
RUNOPT(7) = 'V CELL '
DO 400 J = 1, 3
FA(J) = BOXA(J) / BOX(J)
BOX(J) = BOXA(J)
400 CONTINUE
BOX(4) = BOXA(4)
BOX(5) = BOXA(5)
BOX(6) = BOXA(6)
C ----------------------------------------- Change density
ELSE IF (RUNRUN.EQ.'V DENSITY ') THEN
RUNOPT(7) = 'V DENSITY '
FA(1) = (DENSTY/BOXA(1))**(1.0/3.0)
FA(2) = FA(1)
FA(3) = FA(1)
DO 440 I = 1, 3
BOX(I) = BOX(I) * FA(I)
440 CONTINUE
END IF
C
C -------------------------------------------------- Potential model
READ (15,1000) RUNOPT(8), AMODE, ALPHA
MODE = INT(AMODE)
IF (RUNOPT(8).NE.' ' .AND.
* RUNOPT(8).NE.'BUSING ' .AND.
* RUNOPT(8).NE.'MORSE ' .AND.
* RUNOPT(8).NE.'MORSEQ ' .AND.
* RUNOPT(8).NE.'MORSE-PL ' .AND.
* RUNOPT(8).NE.'MORSE-AT ' .AND.
* RUNOPT(8).NE.'BMH-EXP ' .AND.
* RUNOPT(8).NE.'BMH-EXP* ' .AND.
* runopt(8).ne.'BMH-EXPQ ' .and.
* RUNOPT(8).NE.'BELONO ' .AND.
* RUNOPT(8).NE.'TOSIFUMI ' .AND.
* RUNOPT(8).NE.'WOODCOCK ' .AND.
* RUNOPT(8).NE.'PAULING ' .AND.
* RUNOPT(8).NE.'STSUNE ' .AND.
* RUNOPT(8).NE.'L-J ' .AND.
* RUNOPT(8).NE.'METAL ' ) THEN
WRITE (*,*) 'Interatomic potential model ',
* RUNOPT(8),' is not recognized'
STOP
END IF
C
ZSUM = 0.0
DO 110 I = 1, LEM
ATOM(I) = ' '
ZIO(I) = 0.0
WIO(I) = 0.0
AIO(I) = 0.0
BIO(I) = 0.0
CIO(I) = 0.0
DIO(I) = 0.0
NION(I) = 0
IION(I) = 0
110 CONTINUE
NCOMPO = 0
C --------------------------------------------- Read atom parameters
DO 220 J = 1, LEL+1
READ (15,1300,END=230) I,ATY,AAX,ANJ,ZJ,WJ,AJ,BJ,CJ,DJ
IF (I.LE.0.OR.AAX.EQ.' ') GO TO 230
ATOM(I) = AAX
ZIO(I) = ZJ
WIO(I) = WJ
AIO(I) = AJ
BIO(I) = BJ
CIO(I) = CJ
DIO(I) = DJ
NION(I) = INT(ANJ)
IF (I.NE.1) ZSUM = ZSUM + ZJ * ANJ
IF (ATY.EQ.'-') IION(I) = -1
IF (ATY.EQ.'*') IION(I) = -999
IF (ATY.EQ.'=') IION(I) = 1
NCOMPO = NCOMPO + 1
220 CONTINUE
230 ZI1 = - ZSUM / REAL(NION(1))
IF (ABS(ZI1-ZIO(1)).GT.0.00001) THEN
WRITE (*,*) 'Warnning on total charge neutralization! ',
* ZIO(1),ZI1
C ZIO(1) = ZI1
END IF
IO1 = NCOMPO + 1
DO 240 IO = IO1, LEL
IF (NION(IO).GT.0) NCOMPO = IO
240 CONTINUE
write (6,*) 'Number of components is ',NCOMPO
C ------------------------------------------------------------------
DTMO = DTIME
IF (RUNOPT(2).EQ.'START ') THEN
IF (DDT.GT.0.0001) DTIME = DDT * 1.0E-15
IF (DTIME.LT.1.0E-18) DTIME = 2.0E-15
IF (RUNOPT(17).EQ.'AMORPHOUS '.AND.IP0.EQ.0) THEN
DO 330 I = 1,NTION
DO 330 J = 1, 3
P0(J,I) = P(J,I)
330 CONTINUE
END IF
NAVT = 0
NAV = 0
DO 350 I = 1, LVA
TVAL(I) = 0.0
SVAL(I) = 0.0
VAL0(I) = 0.0
350 CONTINUE
MXCUT = 99999
NRECRD(1) = 0
NRECRD(2) = 0
C VBOX(1) = 1.0
END IF
C
CALL PREPAR (FORMUL)
C
C ---------------------------------------- Configuration and heading
C
NREM = IRECRD(1) - NRECRD(1)
NSTEP1 = NRECRD(1) + 1
THS1 = 'th'
IF (MOD(NSTEP1,10).EQ.1) THS1 = 'st'
IF (MOD(NSTEP1,10).EQ.2) THS1 = 'nd'
IF (MOD(NSTEP1,10).EQ.3) THS1 = 'rd'
THS2 = 'th'
IF (MOD(IRECRD(1),10).EQ.1) THS2 = 'st'
IF (MOD(IRECRD(1),10).EQ.2) THS2 = 'nd'
IF (MOD(IRECRD(1),10).EQ.3) THS2 = 'rd'
WRITE (16, 2000) RUNOPT(2),NREM,NSTEP1,THS1,IRECRD(1),THS2,DTIME,
* IRECRD(2),
* RUNOPT(5),TEMP,DELTMP,NTSTEP,TMPGET,RUNOPT(4),
* NRECRD(2),NRECRD(4)
IF (RUNOPT(5).EQ.'T NOSE ') WRITE (16,2010) STEMP
C
IF (RUNOPT(6).EQ.'P SCALING ') THEN
WRITE (16,2020) RUNOPT(6),(SPRES(I),I=1,3)
ELSE IF (RUNOPT(6).EQ.'P ANDERSEN') THEN
WRITE (16,2027) RUNOPT(6), (SPRES(I),I=1,3),
* (VIRM(LL),LL=1,3)
ELSE IF (RUNOPT(6).EQ.'P SHEAR ') THEN
WRITE (16,2028) RUNOPT(6), (SPRES(I),I=1,3),
* (SPRES(I),I=4,6)
ELSE
WRITE (16,2022) RUNOPT(6)
END IF
C
CALL TABLER (1)
C
C ------------------------------------------ Read RUNOPT(9),...,(22)
lentab = lst
IPRDF(1) = 2
IPRDF(2) = 9999
520 READ (15,1000) RUNRUN,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,PARAM6
IF (RUNRUN.NE.' ') THEN
IF (RUNRUN.EQ.'STRUCTURE ') then ! STRUCTURE [9]
RUNOPT(9) = 'STRUCTURE '
lentab = param1
if (lentab.lt.1) lentab = lst
if (lentab.gt.LST) lentab = lst
end if
IF (RUNRUN.EQ.'NETWORK ') THEN ! NETWORK [10]
RUNOPT(10) = 'NETWORK '
NATX = 0
IO = PARAM1
IF (IO.GT.0.AND.IO.LE.LEE) THEN
NATX = NATX + 1
ATMNET(NATX) = ATOM(IO)
END IF
IO = PARAM2
IF (IO.GT.0.AND.IO.LE.LEE) THEN
NATX = NATX + 1
ATMNET(NATX) = ATOM(IO)
END IF
write (6,*) 'Network forming cation(s) is(are)',
* (i,atmnet(i),i=1,natx)
END IF
C
IF (RUNRUN.EQ.'VELOCITY ') THEN ! VELOCITY [11]
RUNOPT(11) = 'VELOCITY '
IRECRD(9) = PARAM1
PVMULT = 50000.0
IF (PARAM2.GT.0.0) PVMULT = PARAM2
IF (IRECRD(9).LE.0) IRECRD(9) = 1
END IF
IF (RUNRUN.EQ.'POSITION ') THEN ! POSITION [11]
RUNOPT(11) = 'POSITION '
IRECRD(9) = PARAM1
PVMULT = 90000.0
IF (PARAM2.GT.0.0) PVMULT = PARAM2
IF (IRECRD(9).LE.1) IRECRD(9) = 1
END IF
IF (RUNRUN.EQ.'ENERGY ') THEN
RUNOPT(11) = 'ENERGY '
IRECRD(9) = PARAM1
PVMULT = 1.0E12
IF (PARAM2.GT.0) PVMULT = PARAM2
IF (IRECRD(9).LE.1) IRECRD(9) = 1
END IF
IF (RUNRUN.EQ.'POSVELENE ') THEN
RUNOPT(11) = 'POSVELENE '
IRECRD(9) = PARAM1
PVMULT = 1.0E12
C IF (PARAM2.GT.0) PVMULT = PARAM2
IF (IRECRD(9).LE.1) IRECRD(9) = 1
END IF
IF (RUNRUN.EQ.'QUANTUM ') THEN ! QUANTUM [12]
RUNOPT(12) = 'QUANTUM '
CALL QCTABL
END IF
IF (RUNRUN.EQ.'PCF '.OR. ! PCF or RDF.[13]
* RUNRUN.EQ.'RDF ') THEN
RUNOPT(13) = 'PCF '
IF (PARAM1.GT.0.999) IPRDF(1) = PARAM1
IF (PARAM2.GT.0.5 .AND. PARAM2.LT.20.0)
* IPRDF(2) = PARAM2*100
END IF
IF (RUNRUN.EQ.'DIPOLE ') THEN ! DIPOLE [14]
RUNOPT(14) = 'DIPOLE '
END IF
IF (RUNRUN.EQ.'CENTER ') THEN ! CENTER [15]
RUNOPT(15) = 'CENTER '
END IF
IF (RUNRUN.EQ.'NO(MV=0) ') THEN ! NO(MV=0) [16]
RUNOPT(16) = 'NO(MV=0) '
END IF
IF (RUNRUN.EQ.'CRYSTAL ') THEN ! CRYSTAL [17]
RUNOPT(17) = 'CRYSTAL '
END IF
IF (RUNRUN.EQ.'BINARY ') THEN ! BINARY [18]
RUNOPT(18) = 'BINARY '
IF (RUNOPT(2).EQ.'START ') THEN
CLOSE (19)
OPEN (19, FILE=FLNAME(9), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='UNFORMATTED' )
END IF
END IF
IF (RUNRUN.EQ.'PRESSURE ') THEN ! PRESSURE [19]
RUNOPT(19) = 'PRESSURE '
OPEN (27, FILE=FLNAME(13), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
REWIND 27
END IF
IF (RUNRUN.EQ.'ELEC.FIELD') THEN ! ELEC.FIELD [20]
RUNOPT(20) = 'ELEC.FIELD'
MEFD = INT(PARAM1) ! Mode of elec.field
EFD(1) = DBLE(PARAM2) *1.00D5 ! [EFD]==[V/m]
EFD(2) = DBLE(PARAM3) *1.00D5 ! 1 CV/m = 1 J/m
EFD(3) = DBLE(PARAM4) *1.00D5 ! = 10^5 erg/cm
EFREQ = DBLE(PARAM5) ! Hz
c write(6,*) MEFD, EFREQ
c write(6,*) EFD(1),EFD(2),EFD(3)
END IF
if (runrun.eq.'GRAV.FIELD') then ! GRAV.FIELD [21]
runopt(21) = 'GRAV.FIELD'
gfd(1) = param1
gfd(2) = param2
gfd(3) = param3
end if
if (runrun.eq.'CONSTSHEAR') then ! CONSTSHEAR [22]
runopt(22) = 'CONSTSHEAR'
C ----- Shear rate / ps ( dvx/dry )
C ( dvy/drz )
C ( dvx/drz )
STRT(1) = param1
STRT(2) = param2
STRT(3) = param3
IF (RUNOPT(6).EQ.'P SCALING '.OR.
* RUNOPT(6).EQ.'P ANDERSEN'.OR.
* RUNOPT(6).EQ.'P SHEAR ' )then
write (6,*) 'Error ',runopt(6),runopt(22)
stop
end if
end if
if (runrun.eq.'DIATOMIC ') then ! Diatomic molecule ========
runopt(23) = 'DIATOMIC '
write (6,*) param1,param2,param3
DINTRA2(1) = param2
IATOM2(1) = param3
zmole2(1) = param1-zio(iatom2(1))*2.0
MOLstart = param3
MOLend = param3
if (param6.gt.0.0001) then
READ (15,1000) RUNRUN,PARAM1,PARAM2,PARAM3,PARAM4,
* PARAM5,PARAM6
DINTRA2(2) = param2
IATOM2(2) = param3
zmole2(2) = param1-zio(iatom2(2))*2.0
MOLstart = param3
MOLend = param3
end if
CALL DIATOM
write (16,7011) atom(MOLstart), zmole2(1),
* zio(iatom2(1))*2+zmole2(1)
7011 format ('I Diatomic molecule : ',A2, '2 : ',
* 'Charge at molecular center is ',F8.4,
* ', molecular charge is',f8.4,32x, 'I')
if (iatom2(2).gt.0) then
write (16,7012) atom(MOLstart),zmole2(2),
* zio(iatom2(2))*2+zmole2(2)
7012 format ('I : ',A2, '2 : ',
* 'Charge at molecular center is ',F8.4,
* ', molecular charge is',f8.4,32x, 'I')
end if
end if
if (runrun.eq.'TRIATOMIC ') then ! Triatomic molecule =======
runopt(33) = 'TRIATOMIC ' ! 1st 3 atom mol. ex. H2O
DINTRA3(1) = param2 ! ex. O-H
IATOM3(1,1) = param3 ! Center atom O
IATOM3(1,2) = param4 ! H
if (param6.gt.0.0001) then ! 2nd 3 atom mol. ex. CO2
READ (15,1000) RUNRUN,PARAM1,PARAM2,PARAM3,
* PARAM4,PARAM5,PARAM6
DINTRA3(2) = param2 ! C-O
IATOM3(2,1) = param3 ! Center atom C
IATOM3(2,2) = param4 ! O
end if
call triatom
end if
if (runrun.eq.'MOLECULE ') then
runopt(26) = 'MOLECULE '
dMOLintra = param1
MOLstart = param2
MOLend = param3
call MOLECULE
end if
if (runrun.eq.'POLYATOMS ') then
runopt(29) = 'POLYATOMS '
dMOLintra = param1
MOLstart = param2
MOLend = param3
call MOLECULE
end if
GOTO 520
END IF
WRITE (16,2030) (I,RUNOPT(I),I=1,40)
C ---------------------------------------------------- Check P and V
CALL CHECKP (DTMO)
C ------------------------------------------------------ file09p.dat
IF (RUNOPT(2).EQ.'START ') THEN
IF (TITLE(1).NE.'BENC' .OR.
* TITLE(2).NE. 'HMAR' ) THEN
IF (RUNOPT(17).EQ.'AMORPHOUS ') THEN
NRECRD(4) = 1
IF (RUNOPT(18).EQ.'BINARY ') THEN
WRITE (19) NRECRD(4), 0, ((H(J,I),J=1,3),I=1,3)
WRITE (19) ((P(J,I),J=1,3),I=1,NTION)
ELSE
DO 450 I = 1, NTION
DO 450 J = 1, 3
IPV(J,I) = P(J,I) * 90000.0
450 CONTINUE
DUMMY = ' '
WRITE (19,9002) NRECRD(4), 0,
* ((H(J,I),J=1,3),I=1,3)
WRITE (19,9001) ((IPV(J,I),J=1,3),I=1,NTION)
END IF
END IF
END IF
END IF
C ----------------------------------------------------- file09PV.dat
IF (RUNOPT(11).NE.' ') THEN
IF (RUNOPT(18).EQ.'BINARY ') THEN
OPEN (28, FILE=FLNAME(12), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='UNFORMATTED' )
ELSE
OPEN (28, FILE=FLNAME(12), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
END IF
REWIND 28
NRECRD(9) = 1
IF (RUNOPT(11).EQ.'VELOCITY ') THEN
IF (RUNOPT(18).EQ.'BINARY ') THEN
DO 550 I = 1, NTION
DO 550 J = 1, 3
VV(J,I) = V(J,I) / DTIME
550 CONTINUE
WRITE (28) NRECRD(1),IRECRD(9)
WRITE (28) ((VV(J,I),J=1,3),I=1,NTION)
ELSE
DO 560 I = 1, NTION
DO 560 J = 1, 3
IPV(J,I)= V(J,I)*PVMULT*1E-15/DTIME +5000.0
560 CONTINUE
WRITE (28,9002) NRECRD(1),IRECRD(9)
WRITE (28,9001) ((IPV(J,I),J=1,3),I=1,NTION)
END IF
END IF
IF (RUNOPT(11).EQ.'POSITION ') THEN
IF (RUNOPT(18).EQ.'BINARY ') THEN
WRITE (28) NRECRD(1),IRECRD(9), H
WRITE (28) ((P(J,I),J=1,3),I=1,NTION)
ELSE
DO 580 I = 1, NTION
DO 580 J = 1, 3
IPV(J,I) = P(J,I) * PVMULT
580 CONTINUE
WRITE (28,9002) NRECRD(1),IRECRD(9), H
WRITE (28,9001) ((IPV(J,I),J=1,3),I=1,NTION)
END IF
END IF
9001 FORMAT (18I5)
9002 FORMAT (I7,i3, 9F7.3)
END IF
C ------------------------------------------------------------------
IF (NREM.LE.0) GO TO 2222
CALL TITLET (0, 1)
RETURN
C
2222 WRITE (*,2233) RUNOPT(2)
2233 FORMAT ('>>>>> The number of steps to be calculated is less',
* ' than one >>>>>' /
* '>>>>> Mode=', A9, ' Please increase the number ',
* 'of steps >>>>>' )
STOP
C
1000 FORMAT (A10,6F10.5)
1001 FORMAT (A10,15A4)
1300 FORMAT (I1,A1,A2, F6.0,6F10.0)
2000 FORMAT ('I [ ',A10,' ] ',I7,' steps-run from',I7,'-',A2,
* ' to ',I7,'-',A2,' step with time step of',
* 1PE9.2,' sec. RDF''s at every', I7,' step I' /
* 'I [ ',A10,' ] Temperature=',0PF7.1,' K changed ',
* 'with a rate of',F6.1,' K per ', I3, ' steps until',
* F7.1,' K (',A8,' : ',I5,' : ',I4,') I' )
2010 FORMAT ('I',18X,'"Mass" of Nose''s thermostat is ',E12.4,
* ' g.cm2',63X,'I' )
2022 FORMAT ('I [ ',A10,' ] MD basic cell is fixed at the present ',
* 'size and shape. ', 57X, 'I')
2020 FORMAT ('I [ ',A10,' ] Pressure is controlled at Px=',F8.4,
* ' Py=',F8.4,
* ' Pz=',F8.4,
* ' GPa using forced scaling of cell ',
* 'dimensions.',5X,'I')
2027 FORMAT ('I [ ',A10,' ] Pressure is controlled at ',3F9.4,
* ' GPa by Andersen''s mass ',3(1X,G9.2E3),
* ' g I')
2028 FORMAT ('I [ ',A10,' ] Pressure is controlled at Px=',F9.4,
* ' Py=',F9.4,
* ' Pz=',F9.4,
* ' GPa using forced scaling of cell ',
* 'dimensions.',2X,'I'/
* 'I ',10X,30X, 'Pyz=',F8.4,' Pzx=',F8.4,' Pxy=',F8.4,
* ' GPa',44X, 'I' )
2030 format ('I',130('-'),'I' /
* 'I [Options] ',8(I3,':',A10),' I' /
* 'I ',8(I3,':',A10),' I' /
* 'I ',8(I3,':',A10),' I' /
* 'I ',8(I3,':',A10),' I' /
* 'I ',8(I3,':',A10),' I' )
END
C
C
C ==========
C============================================================== MOLECULE
SUBROUTINE MOLECULE
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C ======================================recognize diatomic molecules
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
c
COMMON /MOLECU/ ZMOLE2(2), DMOLE2(4,LNI), DINTRA2(2),
* ZMOLE3(2), DMOLE3(4,LNI), DINTRA3(2),
* dMOLintra,
* NDMOLE2, IDMOLE2(3,LNI), IATOM2(2), MOLstart,
* ndmole3, idmole3(4,lni), iatom3(2,2),
* NMOLE, IMOLE(38,LNI), MMOLE(LNI), MOLend
real *8 zmole2,zmole3,dmole2,dmole3,dintra2,dintra3
c
real *8 pix,piy,piz, pjx,pjy,pjz, rx,ry,rz, dx,dy,dz,
* pjx0,pjy0,pjz0, rij2
integer mi(lni), ndistr(38)
c
cut2 = dMOLintra**2
do 10 I = 1, ntion
mi(i) = 0
10 continue
do 20 n = 1, 38
ndistr(n) = 0
20 continue
nnn = 1 ! No. of molecules
imole(1,nnn) = ions(1,MOLstart)
mi(ions(1,MOLstart)) = 1
mmole(nnn) = 1 ! No. of atoms in the molecule
C------------------------------------------- calc distance between atoms
do 590 io = MOLstart, MOLend
do 510 i = ions(1,io), ions(2,io)
if (mi(i).gt.0) go to 510
c
do 500 n = 1, nnn
do 400 k = 1, mmole(n)
j=imole(k,n)
if (i.eq.j) go to 510
RX = P(1,i) - P(1,j)
RY = P(2,i) - P(2,j)
RZ = P(3,i) - P(3,j)
if (RX.lt.-0.5) RX = RX + 1.0
if (RX.gt. 0.5) RX = RX - 1.0
if (RY.lt.-0.5) RY = RY + 1.0
if (RY.gt. 0.5) RY = RY - 1.0
if (RZ.lt.-0.5) RZ = RZ + 1.0
if (RZ.gt. 0.5) RZ = RZ - 1.0
c --------- delete these if-statements for triclinic
c IF (ABS(RX).GT.0.5) RX = RX - SIGN(1.0D0,RX)
c IF (ABS(RY).GT.0.5) RY = RY - SIGN(1.0D0,RY)
c IF (ABS(RZ).GT.0.5) RZ = RZ - SIGN(1.0D0,RZ)
DX = H(1,1)*RX + H(1,2)*RY + H(1,3)*RZ
DY = H(2,1)*RX + H(2,2)*RY + H(2,3)*RZ
DZ = H(3,1)*RX + H(3,2)*RY + H(3,3)*RZ
c DX = RX * BOX(1)
c DY = RY * BOX(2)
c DZ = RZ * BOX(3)
RIJ2 = DX*DX + DY*DY + DZ*DZ
IF (RIJ2.gt.CUT2) GO TO 400
mmole(n) = mmole(n) + 1
IMOLE(mmole(n),n) = i
mi(i) = 1
go to 510
400 CONTINUE
c
500 continue
nnn=nnn+1
imole(1,nnn) = i
mi(i)=1
mmole(nnn) = 1
510 CONTINUE
590 continue
c
c write (6,*) (mmole(n),n=1,nnn)
c
do 660 n2=2, nnn
mm2=mmole(n2)
if (mm2.le.0) go to 660
do 650 n1 = 1, n2-1
mm1=mmole(n1)
mm2=mmole(n2)
if (mm1.le.0) go to 650
do 630 m1=1, mm1
do 640 m2=1, mm2
i=imole(m1,n1)
j=imole(m2,n2)
RX = P(1,i) - P(1,j)
RY = P(2,i) - P(2,j)
RZ = P(3,i) - P(3,j)
if (RX.lt.-0.5) RX = RX + 1.0
if (RX.gt. 0.5) RX = RX - 1.0
if (RY.lt.-0.5) RY = RY + 1.0
if (RY.gt. 0.5) RY = RY - 1.0
if (RZ.lt.-0.5) RZ = RZ + 1.0
if (RZ.gt. 0.5) RZ = RZ - 1.0
c --------- delete these if-statements for triclinic
c IF (ABS(RX).GT.0.5) RX = RX - SIGN(1.0D0,RX)
c IF (ABS(RY).GT.0.5) RY = RY - SIGN(1.0D0,RY)
c IF (ABS(RZ).GT.0.5) RZ = RZ - SIGN(1.0D0,RZ)
DX = H(1,1)*RX + H(1,2)*RY + H(1,3)*RZ
DY = H(2,1)*RX + H(2,2)*RY + H(2,3)*RZ
DZ = H(3,1)*RX + H(3,2)*RY + H(3,3)*RZ
c DX = RX * BOX(1)
c DY = RY * BOX(2)
c DZ = RZ * BOX(3)
RIJ2 = DX*DX + DY*DY + DZ*DZ
IF (RIJ2.le.CUT2) then
mmm1=mmole(n1)
do m=1, mm2
imole(mmm1+m,n1)=imole(m,n2)
mmole(n1)=mmm1+mm2
mmole(n2)=0
end do
go to 660
end if
640 continue
630 continue
650 continue
660 continue
c
c
nmole=0
do n=1, nnn
na = mmole(n)
if (na.gt.38) na=38
if (na.gt.0) then
ndistr(na)=ndistr(na)+1
nmole=nmole+1
mmole(nmole)=mmole(n)
do i=1, mmole(n)
imole(i,nmole)=imole(i,n)
end do
end if
end do
c write (6,*) (mmole(n),n=1,nmole)
c
write (6,1001) nmole
1001 format (' Total number of molecules is',I5)
write (6,1002) (n,n=1,38), (ndistr(n),n=1,38)
1002 format ('N.A',19I4 / 3X,19I4 / 'N.M',19I4 / 3x,19I4)
RETURN
END
C
C
C ========
C================================================================ DIATOM
SUBROUTINE DIATOM
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C ======================================recognize diatomic molecules
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
c
COMMON /MOLECU/ ZMOLE2(2), DMOLE2(4,LNI), DINTRA2(2),
* ZMOLE3(2), DMOLE3(4,LNI), DINTRA3(2),
* dMOLintra,
* NDMOLE2, IDMOLE2(3,LNI), IATOM2(2), MOLstart,
* ndmole3, idmole3(4,lni), iatom3(2,2),
* NMOLE, IMOLE(38,LNI), MMOLE(LNI), MOLend
real *8 zmole2,zmole3,dmole2,dmole3,dintra2,dintra3
real *8 pix,piy,piz, pjx,pjy,pjz, rx,ry,rz, dx,dy,dz,
* pjx0,pjy0,pjz0, rij2
c
C---------------------------------------------calc distance of atoms
nnn = 0
do 900 iii = 1, 2
cut2 = dintra2(iii)**2
io = iatom2(iii)
if (io.le.0 .or. io.gt.ncompo) go to 900
i1 = ions(1,io)
i2 = ions(2,io)
DO 810 I=i1, i2-1
pix = p(1,i)
piy = p(2,i)
piz = p(3,i)
do 800 J=i+1,i2
pjx0 = p(1,j)
pjy0 = p(2,j)
pjz0 = p(3,j)
if (pjx0.lt.pix) pjx0 = pjx0 + 1.0
if (pjy0.lt.piy) pjy0 = pjy0 + 1.0
if (pjz0.lt.piz) pjz0 = pjz0 + 1.0
DO 250 K = 1, 8
pjx = pjx0 - transx(k)
pjy = pjy0 - transy(k)
pjz = pjz0 - transz(k)
RX = PIX - PjX
RY = PIY - PjY
RZ = PIZ - PjZ
c - - - - - delete these if-statements for triclinic
C IF (ABS(RX).GT.0.5) RX = RX - SIGN(1.0D0,RX)
C IF (ABS(RY).GT.0.5) RY = RY - SIGN(1.0D0,RY)
C IF (ABS(RZ).GT.0.5) RZ = RZ - SIGN(1.0D0,RZ)
DX = H(1,1)*RX + H(1,2)*RY + H(1,3)*RZ
DY = H(2,1)*RX + H(2,2)*RY + H(2,3)*RZ
DZ = H(3,1)*RX + H(3,2)*RY + H(3,3)*RZ
c DX = RX * BOX(1)
c DY = RY * BOX(2)
c DZ = RZ * BOX(3)
RIJ2 = DX*DX + DY*DY + DZ*DZ
IF (RIJ2.LE.CUT2) GO TO 255
250 CONTINUE
go to 800
C ----------------------------------Kumiawase of diatomic
255 nnn = nnn +1
IDMOLE2(1,nnn) = I
IDMOLE2(2,nnn) = J
idmole2(3,nnn) = iii
DMOLE2(1,nnn) = DX
DMOLE2(2,nnn) = Dy
DMOLE2(3,nnn) = DZ
DMOLE2(4,nnn) = SQRT(RIJ2)
C -----------------------------------P of center of mass
Pix=(Pix+Pjx)/2.
Piy=(Piy+Pjy)/2.
Piz=(Piz+Pjz)/2.
if (pix.lt.0.0) pix = pix + 1.0
if (pix.gt.1.0) pix = pix - 1.0
if (piy.lt.0.0) piy = piy + 1.0
if (piy.gt.1.0) piy = piy - 1.0
if (piz.lt.0.0) piz = piz + 1.0
if (piz.gt.1.0) piz = piz - 1.0
p(1,ntion+nnn) = pix
p(2,ntion+nnn) = piy
p(3,ntion+nnn) = piz
C
C WRITE(*,*) nnn,IDMOLE2(1,nnn),IDMOLE2(2,nnn),
C * pix,piy,piz
C
800 CONTINUE
810 continue
900 CONTINUE
ndmole2 = nnn
RETURN
END
C
C
C ==========
C=============================================================== TRIATOM
SUBROUTINE TRIATOM
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C ======================================recognize diatomic molecules
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
c
COMMON /MOLECU/ ZMOLE2(2), DMOLE2(4,LNI), DINTRA2(2),
* ZMOLE3(2), DMOLE3(4,LNI), DINTRA3(2),
* dMOLintra,
* NDMOLE2, IDMOLE2(3,LNI), IATOM2(2), MOLstart,
* ndmole3, idmole3(4,lni), iatom3(2,2),
* NMOLE, IMOLE(38,LNI), MMOLE(LNI), MOLend
real *8 zmole2,zmole3,dmole2,dmole3,dintra2,dintra3
real *8 pix,piy,piz, pjx,pjy,pjz, rx,ry,rz, dx,dy,dz,
* pjx0,pjy0,pjz0, rij2
c
C---------------------------------------------calc distance of atoms
nnn = 0
do 900 iii = 1, 2
cut2 = dintra3(iii)**2
io = iatom3(iii,1)
jo = iatom3(iii,2)
if (io.le.0 .or. io.gt.ncompo) go to 900
if (jo.le.0 .or. jo.gt.ncompo) go to 900
i1 = ions(1,io)
i2 = ions(2,io)
j1 = ions(1,jo)
j2 = ions(2,jo)
DO 810 I=i1, i2
pix = p(1,i)
piy = p(2,i)
piz = p(3,i)
k1=0
k2=0
mmm = 0
do 800 J=j1, j2
pjx0 = p(1,j)
pjy0 = p(2,j)
pjz0 = p(3,j)
if (pjx0.lt.pix) pjx0 = pjx0 + 1.0
if (pjy0.lt.piy) pjy0 = pjy0 + 1.0
if (pjz0.lt.piz) pjz0 = pjz0 + 1.0
DO 250 K = 1, 8
pjx = pjx0 - transx(k)
pjy = pjy0 - transy(k)
pjz = pjz0 - transz(k)
RX = PIX - PjX
RY = PIY - PjY
RZ = PIZ - PjZ
c - - - - - delete these if-statements for triclinic
C IF (ABS(RX).GT.0.5) RX = RX - SIGN(1.0D0,RX)
C IF (ABS(RY).GT.0.5) RY = RY - SIGN(1.0D0,RY)
C IF (ABS(RZ).GT.0.5) RZ = RZ - SIGN(1.0D0,RZ)
DX = H(1,1)*RX + H(1,2)*RY + H(1,3)*RZ
DY = H(2,1)*RX + H(2,2)*RY + H(2,3)*RZ
DZ = H(3,1)*RX + H(3,2)*RY + H(3,3)*RZ
c DX = RX * BOX(1)
c DY = RY * BOX(2)
c DZ = RZ * BOX(3)
RIJ2 = DX*DX + DY*DY + DZ*DZ
IF (RIJ2.GT.CUT2) GO TO 250
if (mmm.eq.0) then
k1=j
dx1=dx
dy1=dy
dz1=dz
pjx1=pjx
pjy1=pjy
pjz1=pjz
dij1=sqrt(rij2)
end if
if (mmm.eq.1) then
k2=j
dx2=dx
dy2=dy
dz2=dz
pjx2=pjx
pjy2=pjy
pjz2=pjz
dij2=sqrt(rij2)
end if
if (mmm.ge.2) then
write (6,*) 'Broken structure > 2'
write (6,*) i, pix, piy, piz
write (6,*) k1, pjx1,pjy1,pjz1,dij1
write (6,*) k2, pjx2,pjy2,pjz2,dij2
write (6,*) j, pjx, pjy, pjz, sqrt(rij2)
do l=1,8
write (6,*) transx(l),transy(l),transz(l)
end do
stop
end if
mmm = mmm + 1
250 CONTINUE
800 continue
if (mmm.ne.2) then
write (6,*) 'Broken structure < 2'
stop
end if
C ----------------------------Atoms in Triatomic molecule
nnn = nnn +1
IDMOLE3(1,nnn) = I
IDMOLE3(2,nnn) = K1
IDMOLE3(3,nnn) = K2
IDmole3(4,nnn) = iii
C -----------------------------------P of center of mass
Pix=(Pix+Pjx1+pjx2)/3.0
Piy=(Piy+Pjy1+pjy2)/3.0
Piz=(Piz+Pjz1+pjz2)/3.0
C
C WRITE(*,*) nnn,IDMOLE2(1,nnn),IDMOLE2(2,nnn),
C * pix,piy,piz
C
810 continue
900 CONTINUE
ndmole3 = nnn
write (6,*) 'ndmole3=',ndmole3
RETURN
END
C
C
C ========
C================================================================ PREPAR
SUBROUTINE PREPAR (FORMUL)
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C ----------------------------------- Preparing some variables, etc.
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
COMMON /VALUES/ VAL(LVA), TVAL(LVA),TVALL(LVA),VALMAX(LVA),
* VAL0(LVA),SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* AVA(LVA,L50), NAV,NAVT
REAL *8 VAL,TVAL,TVALL,VALMAX,VAL0,SVAL,SVALL,VALMIN
COMMON /RADIAL/ NRDF(LTB,LEE),IPRDF(2)
INTEGER *4 NRDF
COMMON /GEOMET/ DTO(2),AVTHT(12),MBR(8,8,2),NRG(9,2),ITBR(121,12),
* RTO(2),SVTHT(12),NBR(8,8,2),MEB(9,2), NTT(121,12),
* NTO(2),NVTHT(12),ANGL(3,12),TTAB(LST), NTBL
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
NELM = 0
TWEGHT = 0.0D0
DO 260 IO = 1, LEL
IONS(1,IO) = NELM + 1
NELM = NELM + NION(IO)
IONS(2,IO) = NELM
NIOND(IO) = 0
DO 250 J = IONS(1,IO), IONS(2,IO)
IF (IOND(J).NE.0) NIOND(IO) = NIOND(IO) + 1
250 CONTINUE
TWEGHT = TWEGHT + WIO(IO) * REAL(NIOND(IO))
260 CONTINUE
NFORML = NION(2)
IF (NFORML.EQ.0) NFORML = NION(3)
IF (FORMUL.GT.0.0) NFORML = NION(1) / FORMUL
FJMOL = ANA / 1.0D10 / REAL(NFORML)
IF (NELM.GT.NTION) GO TO 4444
IF (NELM.LT.NTION) WRITE (*,1004) NELM,NTION
NTION = NELM
C
DO 500 I = 1, LVA
VALMAX (I) = -9.9D19
VALMIN (I) = 9.9D19
500 CONTINUE
C
TPRE = TEMP
RETURN
C
4444 WRITE (*,4455)
4455 FORMAT ('***** THE NUMBER OF PARTICLES IN FILE05 IS MORE THAN ',
* 'THAT IN FILE07 *****')
STOP
C
1004 FORMAT ('******* Warnning ***** NTION(new)=',I5,' (old)=',
* I5,7('*'))
1111 FORMAT (15A4)
END
C
C
C ========
C================================================================ CHECKP
SUBROUTINE CHECKP (DTMO)
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C ----------------------------------- Preparing some variables, etc.
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
REAL *8 DL,FV,TT,RL,CENTER
C
C ----------------------- Check and correct velocity and momentum
FV = 1.0D0
TT = TEMP
IF (TT.LT.0.001) TT = 0.001
IF ((TMPGET-TEMP)*DELTMP.LT.0.0D0) TEMP = TMPGET
FV = SQRT(TEMP/TT) * (DTIME/DTMO)
DO 370 J = 1, 3
DL = 0.0D0
DO 330 IO = 1, NCOMPO
RL = 0.0D0
IF (NION(IO).GT.0) THEN
I1 = IONS(1,IO)
I2 = IONS(2,IO)
DO 310 I = I1, I2
IF (IOND(I).NE.0) RL = RL + V(J,I)
310 CONTINUE
END IF
DL = DL + RL * WIO(IO)
330 CONTINUE
DL = DL / TWEGHT
IF (RUNOPT(16).EQ.'NO(MV=0) ') THEN
DL = 0.0D0
END IF
DO 350 I = 1, NTION
IF (P(J,I).LT.0.0D0) P(J,I) = P(J,I) + 1.0D0
IF (P(J,I).GE.1.0D0) P(J,I) = P(J,I) - 1.0D0
IF (IOND(I).NE.0) V(J,I) = (V(J,I) - DL) * FV
IF (IOND(I).EQ.0) V(J,I) = 0.0
IF (P(J,I)-P0(J,I).GT. 0.5) P0(J,I) = P0(J,I) + 1.0
IF (P(J,I)-P0(J,I).LT.-0.5) P0(J,I) = P0(J,I) - 1.0
350 CONTINUE
IF (RUNOPT(15).EQ.'CENTER ') THEN
CENTER = 0.0D0
DO 360 I = 1, NTION
CENTER = CENTER + P(J,I)
360 CONTINUE
CENTER = CENTER / NTION - 0.5D0
DO 362 I = 1, NTION
P(J,I) = P(J,I) - CENTER
P0(J,I) = P0(J,I) - CENTER
362 CONTINUE
END IF
370 CONTINUE
C
RETURN
END
C
C
C ========
C================================================================ TABLER
SUBROUTINE TABLER (IPR)
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C --------------------------------------------- Heading of MD output
C Preparing tables for force and energy calculations
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VECTOR/ FNV(LNV), UNV(LNV), PNV(3,3,LNV),
* VEC(3,LNV), ZIA(LEM), ALPHA, UCSELF,UCSeLFI(LEM),
* MODE, NVN, NVEC(3,LNV)
REAL *8 FNV,UNV,PNV,VEC,ZIA,ALPHA,UCSELF,ucselfi
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
CHARACTER *63 LOGO1(18), LOGO2(18), LOGO3(12)
DATA LOGO1 /
*' ******* ************************** ',
*' **** *********** ******** ',
*' ***** ********* ******** ',
*' ****** ********** ********* ',
*' ******* *********** *********',
*' **** *** ************ *********',
*' *** *** *** ********* *********',
*' *** *** *** ********* Oblique *********',
*' *** *** *** ********* *********',
*' *** *** *** ********* ******** ',
*' *** ******* ********* ******* ',
*' **** ***** ********* ******* ',
*' ***** *** ********* ******* ',
*' ***** * ********* ******* ',
*' ******* ********* ****** ',
*' ******** *********** ****** ',
*'*********** ************************ R',
*' '/
DATA LOGO2 /
*'************ ************************* ',
*' ********* ************ ******* ',
*' ******** *********** ******* ',
*' ******* *** ******** ******** ',
*' ****** *** ******** ******** ',
*' ****** *** ******** ********',
*' ****** *** ******** ********',
*' ******** ******** Oblique ********',
*' ****** ******** ********',
*' ******** ******** ******* ',
*' *** ****** ******** ******* ',
*' *** ****** ******** ******* ',
*' *** ****** ******** ******* ',
*' *** ****** ******** ****** ',
*' **** ****** ******** ****** ',
*' ****** ******* ********** ****** ',
*'********** *************************** R',
*' '/
DATA LOGO3 /
*'Ms-Fortran-PowerStation Ver.4.0 Version ',
*'386DX+FPU/486DX/Pentium + NDP-FORTRAN/xxx Version ',
*'LUNA-88K (88100+88200) + f77 Version ',
*'Transputer (T805) + Parallel fortran (3L) Version ',
*'HP 9000 Series (PA-RISC) + f77 Version ',
*'IBM-AIX-FORT Version ',
*'F77 on Sony NEWS-WS Version ',
*'FTN compilar on DN10000 Version ',
*'Hitachi Super Computer (S820-80) Version ',
*'F77 on CRAY Super Computer Version ',
*'DEC Fortran for Windows NT Version ',
*' Version '/
c if (FLNAME(3).eq.'Ms-Fortran ') logo3(1) = logo3(1)
if (FLNAME(3).eq.'NDP-FORTRAN386') logo3(1) = logo3(2)
IF (FLNAME(3).EQ.'LUNA88K ') LOGO3(1) = LOGO3(3)
IF (FLNAME(3).EQ.'PARALLEL-F77 ') LOGO3(1) = LOGO3(4)
IF (FLNAME(3).EQ.'HP-9000 ') LOGO3(1) = LOGO3(5)
if (FLNAME(3).eq.'IBM-AIX-FORT ') logo3(1) = logo3(6)
if (FLNAME(3).eq.'NEWS-F77 ') logo3(1) = logo3(7)
if (FLNAME(3).eq.'DN10000 ') logo3(1) = logo3(8)
if (FLNAME(3).eq.'S820-80 ') logo3(1) = logo3(9)
if (FLNAME(3).eq.'CRAY-F77 ') logo3(1) = logo3(10)
if (FLNAME(3).eq.'DEC Fortran ') logo3(1) = logo3(11)
if (FLNAME(3).eq.'Dummy ') logo3(1) = logo3(12)
C
IF (RUNOPT(17).EQ.'CRYSTAL ') THEN
DO 10 I = 1, 18
LOGO1(I) = LOGO2(I)
10 CONTINUE
END IF
C
IDX = 0
IF (RUNOPT(52).EQ.'H-TENSOR ') IDX =1
CALL TMATRX (IDX)
C
IF (RUNOPT(8).NE.'METAL ') CALL COULMB
C
C -------------------------------------------------------- LOGO mark
IF (IPR.EQ.1) THEN
WRITE (16,5000) (REAL(NION(I))/REAL(NFORML),ATOM(I),I=1,LEM)
WRITE (16,5001) BOX(1),BOX(4),
* BOX(2),BOX(5), LOGO1(1),
* BOX(3),BOX(6), LOGO1(2), LOGO1(3),
* DENSTY, VOL, LOGO1(4), LOGO1(5)
WRITE (16,5002) MODE,NVN, RCUT(2), LOGO1(6),
* RUNOPT(8),ALPHA, RCUT(1), LOGO1(7),
* LOGO1(8), LOGO1(9)
5000 FORMAT('I--', 128('-'), 'I' /
* 'I Formula = ',10(F6.3,A2,1X), 26X,' I' /
* 'I--', 126('-'), '--I' )
5001 FORMAT('I Basic cell : A=',F10.5,' A cos(alpha)=',F9.5,
* 10X,'I ',63X, ' I'/
* 'I B=',F10.5,' A cos(beta )=',F9.5,
* 10X,'I ',A63, ' I'/
* 'I C=',F10.5,' A cos(gamma)=',F9.5,
* 10X,'I ',A63, ' I'/
* 'I--',60('-'),'I ', A63, ' I' /
* 'I Density : ',F12.7,' g/cm3 Cell.Vol :',
* F12.5, 3X,'I ',A63, ' I' /
* 'I--',60('-'),'I ',A63, ' I' )
5002 FORMAT('I P-model : Mode=',I3,' (N(Nv)=',I4,') ',
* 'Rcut(S)=',F7.3,' A I ',A63,' I' /
* 'I ',A8,' : Alpha(EWALD)=',F6.3,' A-1 ',
* 'Rcut(L)=',F7.3,' A I ',A63,' I' /
* 'I--',60('-'),'I ', A63,' I' /
* 'I Atom No Z W A B',
* 7X,'C D I ',A63,' I' )
C
DO 110 I = 1, 8
WRITE (16,5005) I, ATOM(I), NION(I), ZIO(I), WIO(I),
* AIO(I), BIO(I), CIO(I), DIO(I),
* LOGO1(I+9)
5005 FORMAT('I', I3, 2X, A3, I6, F8.3, F7.2, F8.4, 3F8.3,
* ' I ',A63,' I' )
110 CONTINUE
I = 9
WRITE (16,5006) I, ATOM(I), NION(I), ZIO(I), WIO(I),
* AIO(I), BIO(I), CIO(I), DIO(I),
* LOGO3(1),FLNAME(2)
I = 10
WRITE (16,5006) I, ATOM(I), NION(I), ZIO(I), WIO(I),
* AIO(I), BIO(I), CIO(I), DIO(I),
* ' ', ' '
5006 FORMAT('I', I3, 2X, A3, I6, F8.3, F7.2, F8.4, 3F8.3,
* ' I ',A52,A11,' I' )
END IF
C
C ------------------------------------------------------ Short range
IF (RUNOPT(8).EQ.'METAL ') CALL METALP (IPR)
IF (IPR.EQ.1) THEN
r3limax = 0.0
IF (RUNOPT(8).EQ.' ') CALL BUSING
IF (RUNOPT(8).EQ.'BUSING ') CALL BUSING
IF (RUNOPT(8).EQ.'STSUNE ') CALL BUSING
IF (RUNOPT(8).EQ.'MORSE ') CALL MORSEP
if (runopt(8).eq.'MORSEQ ') CALL MORSEQ
IF (RUNOPT(8).EQ.'MORSE-PL ') CALL MORSEP
IF (RUNOPT(8).EQ.'MORSE-AT ') CALL MORSEP
IF (RUNOPT(8).EQ.'BMH-EXP ') CALL BMHEXP
IF (RUNOPT(8).EQ.'BMH-EXP* ') CALL BMHEXP
if (runopt(8).eq.'BMH-EXPQ ') call BMHEXPQ
IF (RUNOPT(8).EQ.'BELONO ') CALL MORSEP
IF (RUNOPT(8).EQ.'TOSIFUMI ') CALL TOSIFU
IF (RUNOPT(8).EQ.'WOODCOCK ') CALL ANGELP
IF (RUNOPT(8).EQ.'PAULING ') CALL ANGELP
IF (RUNOPT(8).EQ.'L-J ') CALL LJMODL
C
IF (RUNOPT(3).EQ.'DETAIL ') THEN
DO 200 I = 70, 300, 10
RIJ = I * 0.01
WRITE (16,6666) RIJ, E0(I)*1E8,
* (E1(I,J)*1E8,J=1,NPAIR)
200 CONTINUE
WRITE (16,6666)
DO 210 I = 70, 300, 10
RIJ = I * 0.01
WRITE (16,6666) RIJ,F0(I),(F1(I,J),J=1,NPAIR)
210 CONTINUE
6666 FORMAT (2X,F5.2,1X,F10.6,1X,10F11.7)
END IF
END IF
C
ECORR = 0.0
VCORR = 0.0
IF (RUNOPT(8).EQ.' ' .OR. RUNOPT(8).EQ.'BUSING ' .OR.
* RUNOPT(8).EQ.'STSUNE ' .OR. RUNOPT(8).EQ.'MORSE ' .OR.
* RUNOPT(8).EQ.'MORSE-PL ' .OR. RUNOPT(8).EQ.'MORSE-AT ' .OR.
* RUNOPT(8).EQ.'BMH-EXP ' .OR. RUNOPT(8).EQ.'BELONO ' .OR.
* RUNOPT(8).EQ.'BMH-EXP* ' .OR. runopt(8).eq.'MORSEQ ' .or.
* runopt(8).eq.'BMH-EXPQ ' .or.
* RUNOPT(8).EQ.'TOSIFUMI ' .OR. RUNOPT(8).EQ.'WOODCOCK ' .OR.
* RUNOPT(8).EQ.'PAULING ' .OR. RUNOPT(8).EQ.'L-J ') THEN
CALL VWCORR
END IF
RETURN
END
C
C
C ========
C================================================================ TMATRX
SUBROUTINE TMATRX (IDX)
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
C
REAL *8 SINA(3), COSA(3), DET, GG, BOXIJ
C
C
C -- (0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,1,0),(1,0,1),(0,1,1),(1,1,1)
C
N = 0
DO 10 I = 0, 1
DO 10 J = 0, 1
DO 10 K = 0, 1
N = N + 1
TRANSX(N) = DBLE(I)
TRANSY(N) = DBLE(J)
TRANSZ(N) = DBLE(K)
10 CONTINUE
C
IF (IDX.NE.0) THEN
DO 50 I = 1, 3
BOX(I) = SQRT(H(1,I)**2 + H(2,I)**2 + H(3,I)**2)
50 CONTINUE
DO 68 I = 1, 3
K1 = 2
K2 = 3
IF (I.EQ.2) THEN
K1 = 1
K2 = 3
ELSE IF (I.EQ.3) THEN
K1 = 1
K2 = 2
END IF
BOXIJ= H(1,K1)*H(1,K2)+H(2,K1)*H(2,K2)+H(3,K1)*H(3,K2)
COSA(I) = BOXIJ / (BOX(K1)*BOX(K2))
BOX(I+3) = COSA(I)
SINA(I) = SQRT(1.0D0 - COSA(I)**2)
68 CONTINUE
GO TO 150
END IF
C
C ---------------------------- cos and sin of alpha, beta, and gamma
DO 120 I = 1, 3
COSA(I) = BOX(I+3)
IF (BOX(I+3).GT.1.0) THEN
COSA(I) = COS(BOX(I+3)*PI/180.0D0)
BOX(I+3) = COSA(I)
END IF
SINA(I) = SQRT(1.0D0 - COSA(I)**2)
120 CONTINUE
C
C ------------------ Transformation matrix from crystal to Cartesian
C
H(1,3) = 0.0D0
H(2,3) = 0.0D0
H(3,3) = BOX(3)
H(1,2) = 0.0D0
H(2,2) = BOX(2)*SINA(1)
H(3,2) = BOX(2)*COSA(1)
H(3,1) = BOX(1)*COSA(2)
cc H(2,1) = BOX(1)*COSA(3)*SINA(1)
cc H(1,1) = BOX(1)*SQRT(1.0D0-COSA(2)**2-(COSA(3)*SINA(1))**2)
H(2,1) = -BOX(1)*(COSA(1)*COSA(2)-COSA(3))/SINA(1)
H(1,1) = BOX(1)*SQRT(1-COSA(1)**2-COSA(2)**2-COSA(3)**2+
* 2*COSA(1)*COSA(2)*COSA(3))/SINA(1)
VOL = H(3,1)*(H(1,2)*H(2,3) - H(2,2)*H(1,3)) -
* H(2,1)*(H(1,2)*H(3,3) - H(3,2)*H(1,3)) +
* H(1,1)*(H(2,2)*H(3,3) - H(3,2)*H(2,3))
IF (VOL.LE.0.0D0) THEN
H(1,1) = - H(1,1)
H(2,1) = - H(2,1)
H(3,1) = - H(3,1)
VOL = - VOL
END IF
C
C WRITE (*,*) H(1,1), H(2,1), H(3,1)
C WRITE (*,*) H(1,2), H(2,2), H(3,2)
C WRITE (*,*) H(1,3), H(2,3), H(3,3)
C WRITE (*,*) VOL
C
C ------------------ Transformation matrix from Cartesian to crystal
C
150 CALL INVERS (H, DET, HINV)
C
C WRITE (*,*) HINV(1,1), HINV(2,1), HINV(3,1)
C WRITE (*,*) HINV(1,2), HINV(2,2), HINV(3,2)
C WRITE (*,*) HINV(1,3), HINV(2,3), HINV(3,3)
C
VOL = DET
DENSTY = TWEGHT / (ANA * VOL * 1.0D-24)
C
C ---------------------------------------------------- Metric tensor
DO 180 I = 1, 3
DO 180 J = 1, 3
GG = 0.0 D0
DO 170 K = 1, 3
GG = GG + H(K,J) * H(K,I)
170 CONTINUE
G(J,I) = GG
180 CONTINUE
CALL INVERS (G, DET, GINV)
C
C --------------------------------------- Reciprocal cell parameters
RBOX(1) = BOX(2)*BOX(3)*SINA(1) / VOL
RBOX(2) = BOX(1)*BOX(3)*SINA(2) / VOL
RBOX(3) = BOX(1)*BOX(2)*SINA(3) / VOL
RBOX(4) = (COSA(2)*COSA(3)-COSA(1)) / (SINA(2)*SINA(3))
RBOX(5) = (COSA(1)*COSA(3)-COSA(2)) / (SINA(1)*SINA(3))
RBOX(6) = (COSA(1)*COSA(2)-COSA(3)) / (SINA(1)*SINA(2))
C ---------------------------------------
IF (RCUT(1).LT.0.01) RCUT(1) = 15.0
IF (RCUT(1).GT.1.0/RBOX(1)/2) RCUT(1) = 1.0/RBOX(1)/2
IF (RCUT(1).GT.1.0/RBOX(2)/2) RCUT(1) = 1.0/RBOX(2)/2
IF (RCUT(1).GT.1.0/RBOX(3)/2) RCUT(1) = 1.0/RBOX(3)/2
NRCUT(1) = INT(RCUT(1)*100.0 + 2.5)
C IF (NRCUT(1).LT.LSR) NRCUT(1) = LSR
IF (MXCUT.GT.NRCUT(1)) MXCUT = NRCUT(1)
IF (RCUT(2).LT.0.01) RCUT(2) = 7.5
IF (RCUT(2).GT.RCUT(1)) RCUT(2) = RCUT(1)
IF (RCUT(2).GT.(LSR-1)*0.01) RCUT(2) = (LSR-1)*0.01
NRCUT(2) = INT(RCUT(2)*100.0 +3.01)
RETURN
END
C
C
C ========
C================================================================ INVERS
SUBROUTINE INVERS (X, DET, XINV)
C -------------------------------------------- Given 3 by 3 matrix X
C Store determinant at D and inverse at Xinv
C
REAL *8 DET, X(3,3), XINV(3,3)
C
DET = X(1,1)*X(2,2)*X(3,3) + X(1,2)*X(2,3)*X(3,1) +
* X(1,3)*X(2,1)*X(3,2) - X(1,3)*X(2,2)*X(3,1) -
* X(1,2)*X(2,1)*X(3,3) - X(1,1)*X(2,3)*X(3,2)
IF (DET.EQ.0.0D0) GO TO 10
XINV(1,1) = (X(2,2)*X(3,3) - X(3,2)*X(2,3)) / DET
XINV(1,2) = (X(3,2)*X(1,3) - X(1,2)*X(3,3)) / DET
XINV(1,3) = (X(1,2)*X(2,3) - X(2,2)*X(1,3)) / DET
XINV(2,1) = (X(2,3)*X(3,1) - X(3,3)*X(2,1)) / DET
XINV(2,2) = (X(3,3)*X(1,1) - X(1,3)*X(3,1)) / DET
XINV(2,3) = (X(1,3)*X(2,1) - X(2,3)*X(1,1)) / DET
XINV(3,1) = (X(2,1)*X(3,2) - X(3,1)*X(2,2)) / DET
XINV(3,2) = (X(3,1)*X(1,2) - X(1,1)*X(3,2)) / DET
XINV(3,3) = (X(1,1)*X(2,2) - X(2,1)*X(1,2)) / DET
RETURN
C --------------------------------------------- TEST FOR SINGULARITY
10 IF (DET.EQ.0) WRITE (*,6180)
6180 FORMAT(5X,'*** The matrix is singular ***')
RETURN
END
C
C
C ========
C================================================================ PTOXYZ
C
SUBROUTINE PTOXYZ (I)
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
C
REAL *8 PX,PY,PZ
C
C -------------------------------- TRANSFORMATION OF ION COORDINATES
C FROM CRYSTAL TO CARTESIAN (X,Y,Z)
C
PX = P(1,I)
PY = P(2,I)
PZ = P(3,I)
Q(1,I) = H(1,1)*PX + H(1,2)*PY + H(1,3)*PZ
Q(2,I) = H(2,1)*PX + H(2,2)*PY + H(2,3)*PZ
Q(3,I) = H(3,1)*PX + H(3,2)*PY + H(3,3)*PZ
C
PX = P0(1,I)
PY = P0(2,I)
PZ = P0(3,I)
Q0(1,I) = H(1,1)*PX + H(1,2)*PY + H(1,3)*PZ
Q0(2,I) = H(2,1)*PX + H(2,2)*PY + H(2,3)*PZ
Q0(3,I) = H(3,1)*PX + H(3,2)*PY + H(3,3)*PZ
RETURN
END
C
C
C ========
C================================================================ XYZTOP
C
SUBROUTINE XYZTOP
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
C
REAL *8 QX,QY,QZ
C
C -------------------------------- TRANSFORMATION OF ION COORDINATES
C FROM CARTESIAN (X,Y,Z) TO CRYSTAL
C
DO 100 I = 1, NTION
QX = Q(1,I)
QY = Q(2,I)
QZ = Q(3,I)
P(1,I) = HINV(1,1)*QX + HINV(1,2)*QY + HINV(1,3)*QZ
P(2,I) = HINV(2,1)*QX + HINV(2,2)*QY + HINV(2,3)*QZ
P(3,I) = HINV(3,1)*QX + HINV(3,2)*QY + HINV(3,3)*QZ
C
QX = Q0(1,I)
QY = Q0(2,I)
QZ = Q0(3,I)
P0(1,I) = HINV(1,1)*QX + HINV(1,2)*QY + HINV(1,3)*QZ
P0(2,I) = HINV(2,1)*QX + HINV(2,2)*QY + HINV(2,3)*QZ
P0(3,I) = HINV(3,1)*QX + HINV(3,2)*QY + HINV(3,3)*QZ
100 CONTINUE
RETURN
END
C
C
C ========
C================================================================ COULMB
SUBROUTINE COULMB
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C ------------------------------------ Table for Coulomb interaction
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VECTOR/ FNV(LNV), UNV(LNV), PNV(3,3,LNV),
* VEC(3,LNV), ZIA(LEM), ALPHA, UCSELF,UCSeLFI(LEM),
* MODE, NVN, NVEC(3,LNV)
REAL *8 FNV,UNV,PNV,VEC,ZIA,ALPHA,UCSELF,ucselfi
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
COMMON /CARTES/ H(3,3),HINV(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
* ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,RBOX,Q,TRANSX,TRANSY,TRANSZ
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
REAL *8 XN,FCT, AL2PI,RIJ,ARIJ,PIAL2,VN2,EXPVN,
* YN,UCT, ELC2,ASP,ERFC,PAAV2,alphal,
* ZN,PCT, Z, X0,X1,X2,X3, Y1,Y2,Y3,Y4
INTEGER *4 MXNV(6)
C MODE 1 2 3 4 5 6
C MAXIMUM of NV**2 7 15 23 28 31 39
DATA MXNV / 7, 11, 15, 23, 28, 31 /
C No. of NVs 40 85 125 230 309 369 510
C
C ----- FUNCTION ERFC(X) : VERSION 5662 IN "COMPUTER APPROXIMATIONS"
DATA X0,X1,X2,X3 / 10.00464, 8.426553, 3.460259, .5623536 /
DATA Y0,Y1,Y2,Y3,Y4/ 10.00464, 19.71558, 15.70229, 6.090749, 1.0/
C
ELC2 = ELC**2
DO 10 I = 10, NRCUT(1)+1
E0(I) = 0.0
F0(I) = 0.0
10 CONTINUE
NVN = 0
UCSELF = 0.0D0
DO 30 IO = 1, LEL
ZIA(IO) = 0.0
30 CONTINUE
az = 0.0
do 40 io = 1, ncompo
az = az + abs(zio(io))
40 continue
IF (MODE.LE.-998. .or. az.lt.0.00001) RETURN
C --------------------------------------- Gaussian (alpha) parameter
MAXNV2 = ABS(MODE)
IF (MAXNV2.LE.6) THEN
IF (MAXNV2.LE.0) MAXNV2 = 1
MAXNV2 = REAL(MXNV(MAXNV2))
END IF
ABC2 = MAXNV2 /(RCUT(1)*2.0)**2 * 1.0001
AB = SQRT(ABC2)
IF (ALPHA.LT.0.001) THEN
ALPHAL = MAXNV2 * 0.064D0 + 3.714D0 +
* RCUT(1) * 2.0 * 0.027D0
ALPHA = ALPHAL / (RCUT(1)*2.0D0)
END IF
C ------------------------------------------------------ Coulomb [1]
AL2PI = 2.0D0 * ALPHA / SQRT(PI)
DO 125 I = 10, NRCUT(1)+3
RIJ = REAL(I) * 0.01D0
ARIJ = 1.0D0 / RIJ
C --- FUNCTION ERFC(X) : VERSION 5662
C --- in "COMPUTER APPROXIMATIONS"
Z = ABS(ALPHA * RIJ)
ERFC = EXP(-Z*Z) *
* (X0+Z*(X1+Z*(X2+Z*X3))) /
* (Y0+Z*(Y1+Z*(Y2+Z*(Y3+Z*Y4))))
ERFC = SIGN(ERFC,Z)
IF (Z.LT.0.0D0) ERFC = 2.0D0 + ERFC
E0(I) = ERFC * (ARIJ*1.0D8) * ELC2
F0(I) = ( AL2PI * EXP(-(ALPHA*RIJ)**2) * RIJ + ERFC ) *
* (ARIJ*1.0D8)**2 * ELC2 * ARIJ
125 CONTINUE
C ------------------------------------------------------ Coulomb [2]
C Generate reciprocal vectors for EWALD summation
C Semi-sphere part only
FCT = 4.0 * ELC2 * 1.0D-8 / (VOL*1.0D-24)
UCT = 2.0 * ELC2 * 1.0D-16 / (PI * VOL*1.0D-24)
PCT = 2.0 * ELC2 * 1.0D-16 / (2.0D0 * PI * VOL*1.0D-24)
PIAL2 = PI**2 / ALPHA**2
IL = INT(BOX(1) * AB + 1.5)
JL = INT(BOX(2) * AB + 1.5)
KL = INT(BOX(3) * AB + 1.5)
IL2 = IL * 2 + 1
JL2 = JL * 2 + 1
KL2 = KL + 1
C
DO 270 II = 1, IL2
KX = IL + 1 - II
DO 260 JJ = 1, JL2
KY = JL + 1 - JJ
DO 250 KK = 1, KL2
KZ = KK - 1
IF (KZ.GT.0) GO TO 230
IF (KY.LT.0) GO TO 250
IF (KY.EQ.0 .AND. KX.LE.0) GO TO 250
230 XN = HINV(1,1)*KX + HINV(2,1)*KY + HINV(3,1)*KZ
YN = HINV(1,2)*KX + HINV(2,2)*KY + HINV(3,2)*KZ
ZN = HINV(1,3)*KX + HINV(2,3)*KY + HINV(3,3)*KZ
VN2 = XN**2 + YN**2 + ZN**2
IF (VN2.GT.ABC2) GO TO 250
NVN = NVN + 1
IF (NVN.GT.LNV) THEN
WRITE (*,9901) ABS(MODE)
9901 FORMAT (' ***** SET [MODE] LESS THAN ',I2,
* ' *****')
STOP
END IF
NVEC(1,NVN) = KX
NVEC(2,NVN) = KY
NVEC(3,NVN) = KZ
VEC(1,NVN) = XN
VEC(2,NVN) = YN
VEC(3,NVN) = ZN
XNN = HINV(1,1)*XN + HINV(1,2)*YN + HINV(1,3)*ZN
YNN = HINV(2,1)*XN + HINV(2,2)*YN + HINV(2,3)*ZN
ZNN = HINV(3,1)*XN + HINV(3,2)*YN + HINV(3,3)*ZN
EXPVN = EXP(- VN2 * PIAL2) / VN2
FNV(NVN) = FCT * EXPVN
UNV(NVN) = UCT * EXPVN
PAAV2 = 2.0D0 * (PIAL2 + 1.0D0/VN2)
PCTEXV = PCT * EXPVN
PNV(1,1,NVN)= PCTEXV* H(1,1)*(HINV(1,1)-PAAV2*XNN*XN)
PNV(2,1,NVN)= PCTEXV* H(1,2)*(HINV(1,2)-PAAV2*XNN*YN)
PNV(3,1,NVN)= PCTEXV* H(1,3)*(HINV(1,3)-PAAV2*XNN*ZN)
PNV(1,2,NVN)= PCTEXV* H(2,1)*(HINV(2,1)-PAAV2*YNN*XN)
PNV(2,2,NVN)= PCTEXV* H(2,2)*(HINV(2,2)-PAAV2*YNN*YN)
PNV(3,2,NVN)= PCTEXV* H(2,3)*(HINV(2,3)-PAAV2*YNN*ZN)
PNV(1,3,NVN)= PCTEXV* H(3,1)*(HINV(3,1)-PAAV2*ZNN*XN)
PNV(2,3,NVN)= PCTEXV* H(3,2)*(HINV(3,2)-PAAV2*ZNN*YN)
PNV(3,3,NVN)= PCTEXV* H(3,3)*(HINV(3,3)-PAAV2*ZNN*ZN)
250 CONTINUE
260 CONTINUE
270 CONTINUE
C ------------------------------------------------------ Coulomb [3]
ASP = - (ALPHA*1.0D8) * ELC2 / SQRT(PI)
DO 310 IO = 1, NCOMPO
UCSELF = UCSELF + DBLE(NION(IO))*ZIO(IO)**2*ASP
UCSeLFI(IO) = DBLE(NION(IO))*ZIO(IO)**2*ASP
ZIA(IO) = 2.0 * ZIO(IO)**2*ASP
310 CONTINUE
RETURN
END
C
C
C ========
C================================================================ VWCORR
SUBROUTINE VWCORR
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C --------- Correction of energy and pressur for Van der Waals terms
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VECTOR/ FNV(LNV), UNV(LNV), PNV(3,3,LNV),
* VEC(3,LNV), ZIA(LEM), ALPHA, UCSELF,UCSeLFI(LEM),
* MODE, NVN, NVEC(3,LNV)
REAL *8 FNV,UNV,PNV,VEC,ZIA,ALPHA,UCSELF,ucselfi
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
real *8 pi4, SATOMS
C
PI4 = 4.0D0 * PI
C BETA = CAL * 1.0D10 / ANA
C IF (RUNOPT(8).EQ.'TOSIFUMI ') BETA = 1.0D-19 * 1.0D7
ECORR = 0.0
VCORR = 0.0
N = 0
DO 230 I = 1, NCOMPO
DO 220 J = 1, I
N = N + 1
SATOMS = NION(I) * NION(J) / VOL * PI4
IF (I.EQ.J) SATOMS = SATOMS / 2.0D0
ECORR = ECORR - SATOMS*CIJ(N) / 3.0D0 / RCUT(1)**3
* - SATOMS*DIJ(N) / 5.0D0 / RCUT(1)**5
VCORR = VCORR - 6.0D0*SATOMS*CIJ(N) / 3.0D0 / RCUT(1)**3
* - 8.0D0*SATOMS*DIJ(N) / 5.0D0 / RCUT(1)**5
IF (RUNOPT(8).EQ.'MORSE-PL ') THEN
ECORR = ECORR - SATOMS*D4IJ(N) / RCUT(1)
* - SATOMS*D7IJ(N) / 4.0 / RCUT(1)**4
VCORR = VCORR - 4.0*SATOMS*D4IJ(N) / RCUT(1)
* - 7.0*SATOMS*D7IJ(N) / 4.0 / RCUT(1)**4
END IF
220 CONTINUE
230 CONTINUE
C WRITE (*,1000) ECORR*FJMOL,
C * VCORR / (3.0D0*VOL*1.0D-24)*1.0D-10
C1000 FORMAT (11X, 'Ecorr=',F7.3,'kJ/mol Pcorr=',F6.3,'GPa')
RETURN
END
C
C
C ========
C================================================================ BUSING
SUBROUTINE BUSING
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C ----------------------- IDA-GILBERT-BUSING type potential function
C BORN-MAYER-HUGGINS type
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
REAL *8 BETA,EX,RIJ,ARIJ,ARB
C
BETA = CAL * 1.0D10 / ANA
C
N = 0
DO 110 I = 1, NCOMPO
II = I
DO 100 J = 1, II
C N = (II-1)*LEL +J -(2*LEL-II)*(II-1)/2
N = N + 1
ZIJ(N) = ZIO(II)*ZIO(J)
AIJ(N) = ABS(AIO(II) + AIO(J))
BIJ(N) = ABS(BIO(II) + BIO(J))
CIJ(N) = CIO(II) * CIO(J) * BETA
DIJ(N) = DIO(II) * DIO(J) * BETA
D4IJ(N) = 0.0
D7IJ(N) = 0.0
IF (RUNOPT(8).EQ.'STSUNE ') THEN
IF (I.EQ.J .AND. ATOM(I).EQ.'SI ') CIJ(N) = 0.0
END IF
100 CONTINUE
110 CONTINUE
C
DO 150 I = 10, NRCUT(2)
RIJ = REAL(I) * 0.01
ARIJ = 1.0 / RIJ
DO 140 J = 1, LEE
E1(I,J) = 0.0
F1(I,J) = 0.0
IF (ABS(AIJ(J)).LT.1.0E-5) GO TO 140
EX = 0.0
IF (BIJ(J).GT.0.0001) THEN
ARB = (AIJ(J) - RIJ) / BIJ(J)
IF (ARB.GT.-128.0) EX = EXP(ARB)
END IF
E1(I,J) = BETA * BIJ(J)*EX
C * - CIJ(J)*ARIJ**6
F1(I,J) = BETA * EX * 1.0D8 * ARIJ
C F1(I,J) = (BETA * EX - 6.0*CIJ(J)*ARIJ**7) *
C * 1.0D8 * ARIJ
140 CONTINUE
150 CONTINUE
C
RETURN
END
C
C
C =======
C================================================================ MORSEP
SUBROUTINE MORSEP
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C ----------------------- IDA-GILBERT-BUSING type potential function
C BORN-MAYER-HUGGINS type
C plus MORSE function
C plus three body
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
COMMON /PMORSE/ DMIJ(LEF), BEIJ(LEF), RSIJ(LEF),
* DM1IJ(LEF), BE1IJ(LEF), DM2IJ(LEF),BE2IJ(LEF)
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
REAL *8 E1M,AM1, RIJ, ELC2,BETA,ARB, EPSIJ(LEF),
* F1M,AM2,ARIJ, EX,ZFORML(LEM), SEPij(LEF)
CHARACTER *40 FMT1, FMT2
C
ELC2 = ELC * ELC
BETA = CAL * 1.0D10 / ANA
C
N3BP = 0
DO 10 I = 1, l3p
I3BP(1,I) = 0
i3BP(2,I) = 0
i3bp(3,i) = 0
10 CONTINUE
NPAIR = NCOMPO * (NCOMPO+1) / 2
N = 0
DO 110 I = 1, NCOMPO
II = I
DO 100 J = 1, II
C N = (II-1)*LEL +J -(2*LEL-II)*(II-1)/2
N = N + 1
ZIJ(N) = ZIO(II) * ZIO(J)
AIJ(N) = ABS(AIO(II) + AIO(J))
BIJ(N) = ABS(BIO(II) + BIO(J))
CIJ(N) = CIO(II) * CIO(J) * BETA
DIJ(N) = 0.0
D4IJ(N) = (DIO(II)*ZIO(J)**2 + DIO(J)*ZIO(II)**2 ) / 2.0D0
* * ELC2 * 1.0D8
D7IJ(N) = 2.0D0 * ZIO(II)*ZIO(J) * DIO(II)*DIO(J)
* * ELC2 * 1.0D8
DMIJ(N) = 0.0
BEIJ(N) = 0.0
RSIJ(N) = 0.0
RSWTCH(N) = 0.0
EPSij(N) = 1.0
SEPij(N) = 1.0
100 CONTINUE
110 CONTINUE
C
IF (RUNOPT(8).EQ.'MORSE '.OR.RUNOPT(8).EQ.'MORSE-AT '.OR.
* RUNOPT(8).EQ.'MORSE-PL '.OR.RUNOPT(8).EQ.'BELONO ') THEN
120 READ (15,5555) IP,JP, KP, DIJP, BEIJP, RSIJP, R3BG
5555 FORMAT (3I2,4X,5F10.0)
IF (IP.NE.0.AND.MOD(IP,10).EQ.0) IP = IP / 10
IF (JP.NE.0.AND.MOD(JP,10).EQ.0) JP = JP / 10
IF (KP.NE.0.AND.MOD(KP,10).EQ.0) KP = KP / 10
IF (IP.GE.1.AND.IP.LE.NCOMPO .AND.
* JP.GE.1.AND.JP.LE.NCOMPO ) THEN
IF (KP.EQ.0) THEN
IF (JP.GT.IP) THEN
IJ = IP
IP = JP
JP = IJ
END IF
N = (IP - 1) * IP / 2 + JP
DMIJ(N) = DIJP
BEIJ(N) = BEIJP
RSIJ(N) = RSIJP
RSWTCH(N) = R3BG
ELSE IF (IP.EQ.KP) THEN
N3BP = N3BP +1
I3BP(1,N3BP) = iP
i3BP(2,N3BP) = jP
i3BP(3,N3BP) = KP
C ------------------------------------ F:kJ/mol
FK3BP(N3BP) = DIJP
ANG3BP(N3BP) = BEIJP
R3BLIM(1,N3BP) = RSIJP
R3BGRD(1,N3BP) = R3BG
IF (ANG3BP(N3BP) .LE.0.01) ANG3BP(N3BP) =90.0
IF (R3BLIM(1,N3BP).LE.0.01) R3BLIM(1,N3BP)= 1.2
IF (R3BGRD(1,N3BP).LE.0.01) R3BGRD(1,N3BP)=20.0
R3BLIM(2,N3BP) = R3BLIM(1,N3BP)
R3BGRD(2,N3BP) = R3BGRD(1,N3BP)
ELSE IF (IP.NE.KP) THEN
N3BP = N3BP +1
I3BP(1,N3BP) = iP
i3BP(2,N3BP) = jP
i3BP(3,N3BP) = KP
C ------------------------------------ F:kJ/mol
FK3BP(N3BP) = DIJP
ANG3BP(N3BP) = BEIJP
R3BLIM(1,N3BP) = RSIJP
R3BGRD(1,N3BP) = R3BG
IF (ANG3BP(N3BP).LE.0.01) ANG3BP(N3BP) =90.0
IF (R3BLIM(1,N3BP).LE.0.01) R3BLIM(1,N3BP)= 1.2
IF (R3BGRD(1,N3BP).LE.0.01) R3BGRD(1,N3BP)=20.0
READ (15,5566) R3BLIM2, R3BGRD2
5566 FORMAT (30X,2F10.0)
IF (R3BLIM2.LE.0.01) R3BLIM2 = R3BLIM(1,N3BP)
IF (R3BGRD2.LE.0.01) R3BGRD2 = R3BGRD(1,N3BP)
R3BLIM(2,N3BP) = R3BLIM2
R3BGRD(2,N3BP) = R3BGRD2
ELSE
STOP 'Something wrong in potetial param.'
END IF
GO TO 120
END IF
if (runopt(8).eq.'BELONO ') then
read (15,5577) zforml
5577 format (10f5.0)
N = 0
DO 131 I = 1, NCOMPO
II = I
DO 130 J = 1, II
N = N + 1
epsij(N) = ABS(zio(II)/zforml(II))*
* abs(zio(J) /zforml(J))
sepij(N) = SQRT(1.0 - epsij(N))
130 CONTINUE
131 CONTINUE
end if
LCOMPO = NCOMPO
IF (LCOMPO.GT.7) LCOMPO = 7
LPAIR = LCOMPO*(LCOMPO+1)/2
FMT1 = '( 3H I ,9X, 3(5X,A2,1H-,A2),90X,1HI )'
FMT2 = '( 3H I ,4X,A4,1X, 3F10.3, 90X,1HI )'
IF (NCOMPO.EQ.3) THEN
FMT1 = '( 3H I ,9X, 6(5X,A2,1H-,A2),60X,1HI )'
FMT2 = '( 3H I ,4X,A4,1X, 6F10.3, 60X,1HI )'
ELSE IF (NCOMPO.EQ.4) THEN
FMT1 = '( 3H I ,9X, 10(5X,A2,1H-,A2), 20X,1HI )'
FMT2 = '( 3H I ,4X,A4,1X, 10F10.3, 20X,1HI )'
ELSE IF (NCOMPO.EQ.5) THEN
FMT1 = '( 3H I ,7X, 15(3X,A2,1H-,A2), 2X,1HI )'
FMT2 = '( 3H I ,2X,A4,1X, 15F8.2, 2X,1HI )'
ELSE IF (NCOMPO.EQ.6) THEN
FMT1 = '( 3H I ,3X, 21(1X,A2,1H-,A2), 1HI )'
FMT2 = '( 3H I ,A3, 21F6.2, 1HI )'
ELSE IF (NCOMPO.GE.7) THEN
FMT1 = '( 3H I ,5X, 28(1X,A1,1H-,A1),12X,1HI )'
FMT2 = '( 3H I ,1X,A4,1X, 28F4.1, 12X,1HI )'
END IF
WRITE (16, 6661)
6661 FORMAT ('I ', 60(' '), 'I--', 63('-'), '--I' )
WRITE (16,FMT1) ((ATOM(I),ATOM(J),J=1,I),I=1,LCOMPO)
WRITE (16,FMT2) 'Dij ', (DMIJ(J),J=1,LPAIR)
WRITE (16,FMT2) 'BEij', (BEIJ(J),J=1,LPAIR)
WRITE (16,FMT2) 'RSij', (RSIJ(J),J=1,LPAIR)
WRITE (16,FMT2) 'Rsw ', (RSWTCH(J),J=1,LPAIR)
if (RUNOPT(8).EQ.'BELONO ') then
write (16,fmt2) 'EPij', (EPSij(J),J=1,LPAIR)
write (16,fmt2) 'SEij', (SEPij(J),J=1,LPAIR)
end if
if (N3BP.GT.0) THEN
WRITE (16,6666)
6666 FORMAT ('I ',60(' '),' ', 63(' '),' I' /
* 'I',5X,'3-body potential ATOM(J)--ATOM(I)',
* '--ATOM(J) FK3BP ANG3BP ',
* ' R3BLIM ',
* ' R3BGRD R3LIM ',15X, 'I')
DO 140 N = 1, N3BP
IF (I3BP(2,N)*i3BP(1,N).GT.0) THEN
R3LIM(1,n) = LOG(0.999999D0/0.000001)/R3BGRD(1,N)
* + R3BLIM(1,N)
r3lim(2,n) = r3lim(1,n)
if (r3limax.lt.r3lim(1,n)) r3limax=r3lim(1,n) /////
WRITE (16,6667) ATOM(i3BP(1,N)), i3BP(1,N),
* ATOM(I3BP(2,N)), I3BP(2,N),
* ATOM(i3BP(3,N)), i3BP(3,N),
* FK3BP(N),ANG3BP(N),i3bp(2,n),i3bp(1,n),
* R3BLIM(1,N), R3BGRD(1,N), R3LIM(1,n)
6667 FORMAT ('I',22X, 3X,A2,'(',I2,')--', A2,'(',
* I2,')--',A2,'(',I2,')', F15.8, F11.3,
* i6,'-',i2, 2F10.3, F12.4,16X, 'I')
if (i3BP(1,N).ne.i3BP(3,N)) then
R3LIM(2,n) = LOG(0.999999D0/0.000001) /
* R3BGRD(2,N) + R3BLIM(2,N)
if (r3limax.lt.r3lim(2,n)) r3limax=r3lim(2,n) /////
WRITE (16,6668) i3bp(2,n),i3bp(3,n),
* R3BLIM(2,N),
* R3BGRD(2,N), R3LIM(2,n)
6668 FORMAT ('I',73X, i6,'-',i2,
* 2F10.3, F12.4,16X, 'I')
end if
END IF
140 CONTINUE
END IF
end if
C
DO 250 I = 10, NRCUT(2)
RIJ = REAL(I) * 0.01
ARIJ = 1.0 / RIJ
DO 240 J = 1, NPAIR
E1(I,J) = 0.0
F1(I,J) = 0.0
E1M = 0.0
F1M = 0.0
IF (ABS(AIJ(J)).LT.1.0E-5) GO TO 220
EX = 0.0
IF (BIJ(J).GT.0.00001) THEN
ARB = (AIJ(J) - RIJ) / BIJ(J)
IF (ARB.GT.-128.0) EX = EXP(ARB)
END IF
E1(I,J) = BETA * BIJ(J)*EX*EPSij(J)
c * - CIJ(J)*ARIJ**6
C * - D4IJ(J)*ARIJ**4 - D7IJ(J)*ARIJ**7
F1(I,J) = BETA * EX*EPSij(J)
c * - 6.0*CIJ(J)*ARIJ**7
C * - 4.0*D4IJ(J)*ARIJ**5 - 7.0*D7IJ(J)*ARIJ**8
C * - 4.0*D4IJ(J)*ARIJ**5 - D4IJ(J)*ARIJ**4/4.43
220 IF (DMIJ(J).LT.0.01) GO TO 230
IF (RUNOPT(8).EQ.'MORSE ' .OR.
* RUNOPT(8).EQ.'MORSE-PL ' .OR.
* RUNOPT(8).EQ.'BELONO ' ) THEN
AM1 = EXP(-2.0*BEIJ(J)*(RIJ-RSIJ(J)))
AM2 = EXP(-1.0*BEIJ(J)*(RIJ-RSIJ(J)))
E1M = BETA*DMIJ(J) *(AM1-2.0*AM2) *SEPij(J)
F1M = BETA*BEIJ(J) * DMIJ(J) * 2.0*(AM1 - AM2)
* *SEPij(J)
END IF
IF (RUNOPT(8).EQ.'MORSE-AT ') THEN
AM2 = DMIJ(J) * EXP(-BEIJ(J)*RIJ)
E1M = - BETA * AM2
F1M = - BETA * BEIJ(J) * AM2
END IF
IF (RSWTCH(J).LT.1.0E-6) THEN
E1(I,J) = E1(I,J) + E1M
F1(I,J) = F1(I,J) + F1M
ELSE IF (RIJ.LE.RSWTCH(J)) THEN
E1(I,J) = E1M
F1(I,J) = F1M
END IF
230 F1(I,J) = F1(I,J)*1.0D8 * ARIJ
240 CONTINUE
250 CONTINUE
RETURN
END
C
C
C =======
C================================================================ MORSEP
SUBROUTINE MORSEQ
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C ----------------------- IDA-GILBERT-BUSING type potential function
C BORN-MAYER-HUGGINS type
C plus MORSE function
C plus three body
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
COMMON /PMORSE/ DMIJ(LEF), BEIJ(LEF), RSIJ(LEF),
* DM1IJ(LEF), BE1IJ(LEF), DM2IJ(LEF),BE2IJ(LEF)
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
REAL *8 E1M,AM1, RIJ, ELC2,BETA,ARB, EPSIJ(LEF),
* F1M,AM2,ARIJ, EX,ZFORML(LEM), SEPij(LEF)
CHARACTER *40 FMT1, FMT2
C
ELC2 = ELC * ELC
BETA = CAL * 1.0D10 / ANA
BETAJ = 1.0D10 / ANA
C
N3BP = 0
DO 10 I = 1, l3p
I3BP(1,I) = 0
i3BP(2,I) = 0
i3bp(3,i) = 0
10 CONTINUE
NPAIR = NCOMPO * (NCOMPO+1) / 2
N = 0
DO 110 I = 1, NCOMPO
II = I
DO 100 J = 1, II
N = N + 1
AIJ(N) = CIO(II) + CIO(J)
BIJ(N) = BIO(II) * BIO(J)
CIJ(N) = AIO(II)*AIO(J)*BETAJ
DIJ(N) = 0.0
D4IJ(N) = 0.0
D7IJ(N) = 0.0
DMIJ(N) = 0.0
BEIJ(N) = 0.0
RSIJ(N) = 0.0
RSWTCH(N) = 0.0
epsij(n) = 1.0
sepij(n) = 1.0
100 CONTINUE
110 CONTINUE
C
IF (RUNOPT(8).EQ.'MORSEQ ') THEN
120 READ (15,5555) IP,JP, KP, ijkl,
* DIJP, BEIJP, RSIJP, R3BG
5555 FORMAT (3I2,i2,2x,5F10.0)
IF (IP.NE.0.AND.MOD(IP,10).EQ.0) IP = IP / 10
IF (JP.NE.0.AND.MOD(JP,10).EQ.0) JP = JP / 10
IF (KP.NE.0.AND.MOD(KP,10).EQ.0) KP = KP / 10
IF (IP.GE.1.AND.IP.LE.NCOMPO .AND.
* JP.GE.1.AND.JP.LE.NCOMPO ) THEN
IF (KP.EQ.0) THEN
IF (JP.GT.IP) THEN
IJ = IP
IP = JP
JP = IJ
END IF
N = (IP - 1) * IP / 2 + JP
DMIJ(N) = DIJP
BEIJ(N) = BEIJP
RSIJ(N) = RSIJP
RSWTCH(N) = R3BG
ELSE IF (IP.EQ.KP) THEN
N3BP = N3BP +1
I3BP(1,N3BP) = IP
i3BP(2,N3BP) = JP
i3BP(3,N3BP) = KP
C -------------------------------------- F:kJ/mol
FK3BP(N3BP) = DIJP
ANG3BP(N3BP) = BEIJP
R3BLIM(1,N3BP) = RSIJP
R3BGRD(1,N3BP) = R3BG
IF (ANG3BP(N3BP).LE.0.01) ANG3BP(N3BP)= 90.0
IF (R3BLIM(1,N3BP).LE.0.01) R3BLIM(1,N3BP)= 1.2
IF (R3BGRD(1,N3BP).LE.0.01) R3BGRD(1,N3BP)=20.0
R3BLIM(2,N3BP) = R3BLIM(1,N3BP)
R3BGRD(2,N3BP) = R3BGRD(1,N3BP)
ELSE IF (IP.NE.KP) THEN
N3BP = N3BP +1
I3BP(1,N3BP) = IP
i3BP(2,N3BP) = jP
i3BP(3,N3BP) = KP
C ------------------------------------ F:kJ/mol
FK3BP(N3BP) = DIJP
ANG3BP(N3BP) = BEIJP
R3BLIM(1,N3BP) = RSIJP
R3BGRD(1,N3BP) = R3BG
IF (ANG3BP(N3BP).LE.0.01) ANG3BP(N3BP) =90.0
IF (R3BLIM(1,N3BP).LE.0.01) R3BLIM(1,N3BP)= 1.2
IF (R3BGRD(1,N3BP).LE.0.01) R3BGRD(1,N3BP)=20.0
READ (15,5566) R3BLIM2, R3BGRD2
5566 FORMAT (30X,2F10.0)
IF (R3BLIM2.LE.0.01) R3BLIM2 = R3BLIM(1,N3BP)
IF (R3BGRD2.LE.0.01) R3BGRD2 = R3BGRD(1,N3BP)
R3BLIM(2,N3BP) = R3BLIM2
R3BGRD(2,N3BP) = R3BGRD2
ELSE
STOP 'Something wrong in potetial param.'
END IF
GO TO 120
END IF
if (runopt(8).eq.'BELONO ') then
read (15,5577) zforml
5577 format (10f5.0)
N = 0
DO 131 I = 1, NCOMPO
II = I
DO 130 J = 1, II
N = N + 1
epsij(N) = ABS(zio(II)/zforml(II))*
* abs(zio(J) /zforml(J))
sepij(N) = SQRT(1.0 - epsij(N))
130 CONTINUE
131 CONTINUE
end if
LCOMPO = NCOMPO
IF (LCOMPO.GT.7) LCOMPO = 7
LPAIR = LCOMPO*(LCOMPO+1)/2
FMT1 = '( 2HI ,9X, 3(5X,A2,1H-,A2),90X,1HI ) '
FMT2 = '( 2HI ,2X,A6,1X, 3F10.3, 90X,1HI ) '
IF (NCOMPO.EQ.3) THEN
FMT1 = '( 2HI ,9X, 6(5X,A2,1H-,A2),60X,1HI ) '
FMT2 = '( 2HI ,2X,A6,1X, 6F10.3, 60X,1HI ) '
ELSE IF (NCOMPO.EQ.4) THEN
FMT1 = '( 2HI ,9X, 10(5X,A2,1H-,A2), 20X,1HI ) '
FMT2 = '( 2HI ,2X,A6,1X, 10F10.3, 20X,1HI ) '
ELSE IF (NCOMPO.EQ.5) THEN
FMT1 = '( 2HI ,7X, 15(3X,A2,1H-,A2), 2X,1HI ) '
FMT2 = '( 2HI ,1X,A5,1X, 15F8.2, 2X,1HI ) '
ELSE IF (NCOMPO.EQ.6) THEN
FMT1 = '( 2HI ,3X, 21(1X,A2,1H-,A2), 1HI ) '
FMT2 = '( 2HI ,A3, 21F6.2, 1HI ) '
ELSE IF (NCOMPO.EQ.7) THEN
FMT1 = '( 2HI ,5X, 28(1X,A1,1H-,A1),12X,1HI ) '
FMT2 = '( 2HI ,A5, 1X, 28F4.1, 12X,1HI ) '
END IF
WRITE (16, 6661)
6661 FORMAT ('I ', 60(' '), 'I--', 63('-'), '--I' )
WRITE (16,FMT1) ((ATOM(I),ATOM(J),J=1,I),I=1,LCOMPO)
WRITE (16,FMT2) 'Dij ', (DMIJ(J),J=1,LPAIR)
WRITE (16,FMT2) 'BEij ', (BEIJ(J),J=1,LPAIR)
WRITE (16,FMT2) 'RSij ', (RSIJ(J),J=1,LPAIR)
write (16,fmt2) 'Rswtch',(RSWTCH(J),j=1,lpair)
if (RUNOPT(8).EQ.'BELONO ') then
write (16,fmt2) 'EPij', (EPSij(J),J=1,LPAIR)
write (16,fmt2) 'SEij', (SEPij(J),J=1,LPAIR)
end if
if (N3BP.GT.0) THEN
WRITE (16,6666)
6666 FORMAT ('I ',60(' '),' ', 63(' '),' I' /
* 'I',5X,'3-body potential ATOM(J)--ATOM(I)',
* '--ATOM(J) FK3BP ANG3BP ',
* ' R3BLIM ',
* ' R3BGRD R3LIM ',15X, 'I')
DO 140 N = 1, N3BP
IF (I3BP(2,N)*i3BP(1,N).GT.0) THEN
R3LIM(1,n) = LOG(0.999999D0/0.000001)/R3BGRD(1,N)
* + R3BLIM(1,N)
r3lim(2,n) = r3lim(1,n)
if (r3limax.lt.r3lim(1,n)) r3limax=r3lim(1,n)
WRITE (16,6667) ATOM(i3BP(1,N)), i3BP(1,N),
* ATOM(I3BP(2,N)), I3BP(2,N),
* ATOM(i3BP(3,N)), i3BP(3,N),
* FK3BP(N),ANG3BP(N),i3bp(2,n),i3bp(1,n),
* R3BLIM(1,N), R3BGRD(1,N), R3LIM(1,n)
6667 FORMAT ('I',22X, 3X,A2,'(',I2,')--', A2,'(',
* I2,')--',A2,'(',I2,')', F15.8, F11.3,
* i6,'-',i2, 2F10.3, F12.4,16X, 'I')
if (i3BP(1,N).ne.i3BP(3,N)) then
R3LIM(2,n) = LOG(0.999999D0/0.000001) /
* R3BGRD(2,N) + R3BLIM(2,N)
if (r3limax.lt.r3lim(2,n)) r3limax=r3lim(2,n)
WRITE (16,6668) i3bp(2,n),i3bp(3,n),
* R3BLIM(2,N),
* R3BGRD(2,N), R3LIM(2,n)
6668 FORMAT ('I',73X, i6,'-',i2,
* 2F10.3, F12.4,16X, 'I')
end if
END IF
140 CONTINUE
END IF
END IF
C
DO 250 I = 10, NRCUT(2)
RIJ = REAL(I) * 0.01
ARIJ = 1.0 / RIJ
DO 240 J = 1, NPAIR
E1(I,J) = 0.0
F1(I,J) = 0.0
E1M = 0.0
F1M = 0.0
EX = BIJ(j)*EXP(-Aij(j)*Rij)
E1(I,J) = BETAj * EX
F1(I,J) = BETAj * AIJ(j)*EX
AM1 = EXP(-2.0*BEIJ(J)*(RIJ-RSIJ(J)))
AM2 = EXP(-1.0*BEIJ(J)*(RIJ-RSIJ(J)))
E1M= BETA*DMIJ(J) *(AM1 - 2.0*AM2) * SEPij(J)
F1M= BETA*BEIJ(J) *DMIJ(J) * (2.0*AM1 -
* 2.0*AM2) * SEPij(J)
IF (RIJ.GT.RSWTCH(j)) THEN
E1(I,J) = E1(I,J)
F1(I,J) = F1(I,J)
ELSE IF (RIJ.LE.RSWTCH(J)) THEN
E1(I,J) = E1M
F1(I,J) = F1M
END IF
230 F1(I,J) = F1(I,J)*1.0D8 * ARIJ
240 CONTINUE
250 CONTINUE
RETURN
END
C
C
C =======
C================================================================ BMHEXP
SUBROUTINE BMHEXP
PARAMETER (LNI=36789, LTB=3004, LEL=8, LEM=10, LCT=2000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=19876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4, LVA=24+LEM*2, L3P=7, LRG=LNI*3 )
C
C ----------------------- IDA-GILBERT-BUSING type potential function
C BORN-MAYER-HUGGINS type plus Expornential type function
C plus gauss type function
C plus three body
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST(3,3),SPRES(6),PPXYZ(7),
* FJMOL, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,TDUMP,PDUMP,
* STEMP,VSTEMP,PREST,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL, RCUT(2),
* VIRM(6),DENSTY, NFORML,NRCUT(2),MXCUT
REAL *8 BOX, VBOX, VOL, RCUT, VIRM,DENSTY
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /PARAMT/ AIO(LEM),BIO(LEM),CIO(LEM),DIO(LEM),ZIO(LEM),
* AIJ(LEF),BIJ(LEF),CIJ(LEF),DIJ(LEF),D4IJ(LEF),
* PLIJ(LEF),RSWTCH(LEF),ZIJ(LEF), D7IJ(LEF),
* ECORR,VCORR, WIO(LEM),TWEGHT, AKFI(LEM),
* ANG3bp(l3p), r3blim(2,l3p),
* FK3bp(l3p), r3bgrd(2,l3p), R3lim(2,l3p),r3limax,
* i3bp(3,l3p), N3BP
COMMON /TABLES/ F1(LSR,LEE),E1(LSR,LEE),F0(LTB),E0(LTB)
REAL *8 F1,E1,F0,E0
COMMON /PMORSE/ DMIJ(LEF), BEIJ(LEF), RSIJ(LEF),
* DM1IJ(LEF), BE1IJ(LEF), DM2IJ(LEF),BE2IJ(LEF)
COMMON /CHARAC/ TITLE(15),RUNOPT(53),ATOM(LEM),ATMNET(2),
* ATMXTL(LAA),FLNAME(19)
CHARACTER *4 TITLE,ATOM,ATMNET,ATMXTL
CHARACTER *10 RUNOPT
CHARACTER *16 FLNAME
COMMON /CONSTS/ PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
REAL *8 PI, ANA, AKB, AHP, EP0, ELC, CVL, CAL
C
REAL *8 E1M,AM1, RIJ, ELC2,BETA,ARB, EPSIJ(LEF),
* F1M,AM2,ARIJ, EX,ZFORML(LEM), SEPij(LEF)
real *8 am3, dm3ij(lef), be3ij(lef), r03ij(lef)
integer ipara(2,10), npara
real *4 apara(8,10)
CHARACTER *54 FMT1, FMT2
C
ELC2 = ELC * ELC
BETA = CAL * 1.0D10 / ANA
C
N3BP = 0
DO 10 I = 1, l3p
I3BP(1,I) = 0
i3BP(2,I) = 0
i3bp(3,i) = 0
10 CONTINUE
NPAIR = NCOMPO * (NCOMPO+1) / 2
N = 0
DO 110 I = 1, NCOMPO
II = I
DO 100 J = 1, II
C N = (II-1)*LEL +J -(2*LEL-II)*(II-1)/2
N = N + 1
ZIJ(N) = ZIO(II) * ZIO(J)
AIJ(N) = ABS(AIO(II) + AIO(J))
BIJ(N) = ABS(BIO(II) + BIO(J))
CIJ(N) = CIO(II) * CIO(J) * BETA
DIJ(N) = 0.0
D4IJ(N) = (DIO(II)*ZIO(J)**2 + DIO(J)*ZIO(II)**2 ) / 2.0D0
* * ELC2 * 1.0D8
D7IJ(N) = 2.0D0 * ZIO(II)*ZIO(J) * DIO(II)*DIO(J)
* * ELC2 * 1.0D8
DM1IJ(N) = 0.0
BE1IJ(N) = 0.0
DM2IJ(N) = 0.0
BE2IJ(N) = 0.0
DM3IJ(N) = 0.0
BE3IJ(N) = 0.0
r03ij(n) = 0.0
RSWTCH(N) = 0.0
EPSij(N) = 1.0
SEPij(N) = 1.0
100 CONTINUE
110 CONTINUE
C
npara = 0
120 READ (15,5555) IP,JP, KP, ijkl,
* D1, BE1, D2, BE2, RSIJP, ggg
5555 FORMAT (3I2,i2,2X,6F10.0)
5556 format (10x,3F10.0)
IF (IP.NE.0.AND.MOD(IP,10).EQ.0) IP = IP / 10
IF (JP.NE.0.AND.MOD(JP,10).EQ.0) JP = JP / 10
IF (KP.NE.0.AND.MOD(KP,10).EQ.0) KP = KP / 10
IF (IP.GE.1.AND.IP.LE.NCOMPO .AND.
* JP.GE.1.AND.JP.LE.NCOMPO ) THEN
IF (KP.EQ.0) THEN
IF (JP.GT.IP) THEN
IJ = IP
IP = JP
JP = IJ
END IF
N = (IP - 1) * IP / 2 + JP
if (ijkl.eq.1) then
AIJ(N) = 0.0
BIJ(N) = 0.0
CIJ(N) = 0.0
DIJ(N) = 0.0
D4IJ(N) = 0.0
D7IJ(N) = 0.0
end if
DM1IJ(N) = D1
BE1IJ(N) = BE1
DM2IJ(N) = D2
BE2IJ(N) = BE2
RSWTCH(N) = RSIJP
if (ggg.gt.0.0) then
read (15,5556) dm3ij(n),be3ij(n),r03ij(n)
end if
npara = npara + 1
ipara(1,npara) = ip
ipara(2,npara) = jp
apara(1,npara) = d1
apara(2,npara) = be1
apara(3,npara) = d2
apara(4,npara) = be2
apara(5,npara) = dm3ij(n)
apara(6,npara) = be3ij(n)
apara(7,npara) = r03ij(n)
apara(8,npara) = rsijp
ELSE IF (IP.EQ.KP) THEN
N3BP = N3BP +1
I3BP(1,N3BP) = iP
i3BP(2,N3BP) = jP
i3BP(3,N3BP) = KP
C ------------------------------------ F:kJ/mol
FK3BP(N3BP) = D1
ANG3BP(N3BP) = BE1
R3BLIM(1,N3BP) = D2
R3BGRD(1,N3BP) = BE2
IF (ANG3BP(N3BP).LE.0.01) ANG3BP(N3BP) =90.0
IF (R3BLIM(1,N3BP).LE.0.01) R3BLIM(1,N3BP)= 1.2
IF (R3BGRD(1,N3BP).LE.0.01) R3BGRD(1,N3BP)=20.0
R3BLIM(2,N3BP) = R3BLIM(1,N3BP)
R3BGRD(2,N3BP) = R3BGRD(1,N3BP)
ELSE IF (IP.NE.KP) THEN
N3BP = N3BP +1
I3BP(1,N3BP) = iP
i3BP(2,N3BP) = jP
i3BP(3,N3BP) = KP
C ------------------------------------ F:kJ/mol
FK3BP(N3BP) = D1
ANG3BP(N3BP) = BE1
R3BLIM(1,N3BP) = D2
R3BGRD(1,N3BP) = BE2
IF (ANG3BP(N3BP).LE.0.01) ANG3BP(N3BP) =90.0
IF (R3BLIM(1,N3BP).LE.0.01) R3BLIM(1,N3BP)= 1.2
IF (R3BGRD(1,N3BP).LE.0.01) R3BGRD(1,N3BP)=20.0
READ (15,5566) R3BLIM2, R3BGRD2
5566 FORMAT (30X,2F10.0)
IF (R3BLIM2.LE.0.01) R3BLIM2 = R3BLIM(1,N3BP)
IF (R3BGRD2.LE.0.01) R3BGRD2 = R3BGRD(1,N3BP)
R3BLIM(2,N3BP) = R3BLIM2
R3BGRD(2,N3BP) = R3BGRD2
ELSE
STOP 'Something wrong in potetial param.'
END IF
GO TO 120
END IF
C
write (16,6661)
6661 format ('I ', 60(' '), 'I--', 63('-'), '--I' /
* 'I',19x,'DM1ij BE1ij DM2ij BE2ij',
* 7x,'DM3ij BE3ij R03ij Rswch',29x,'I' )
if (npara.gt.0) then
do 130 i = 1, npara
WRITE (16, 6663) ATOM(Ipara(1,i)),ipara(1,i),
* ATOM(ipara(2,i)),ipara(2,i), (apara(j,i),j=1,8)
6663 format ('I',2x,A2,'(',i1,')-',A2,'(',i1,') ',
* 3(F12.3, F10.3),F10.3,F10.3, 29x, 'I')
130 continue
end if
C
if (N3BP.GT.0) THEN
WRITE (16,6666)
6666 FORMAT ('I ',60(' '),' ', 63(' '),' I' /
* 'I',5X,'3-body potential ATOM(J)--ATOM(I)',
* '--ATOM(J) FK3BP ANG3BP R3BLIM ',
* ' R3BGRD R3LIM ',24X, 'I')
DO 140 N = 1, N3BP
IF (I3BP(2,N)*i3BP(1,N).GT.0) THEN
R3LIM(1,n) = LOG(0.999999D0/0.000001)/R3BGRD(1,N)
* + R3BLIM(1,N)
if (runopt(8).eq.'BMH-EXP* ') then
R3LIM(1,n) = LOG(0.9999D0/0.0001D0) /
* R3BGRD(1,N) + R3BLIM(1,N)
end if
r3lim(2,n) = r3lim(1,n)
if (r3limax.lt.r3lim(1,n)) r3limax=r3lim(1,n) /////
WRITE (16,6667) ATOM(i3BP(1,N)), i3BP(1,N),
* ATOM(I3BP(2,N)), I3BP(2,N),
* ATOM(i3BP(3,N)), i3BP(3,N),
* FK3BP(N),ANG3BP(N),i3bp(2,n),i3bp(1,n),
* R3BLIM(1,N), R3BGRD(1,N), R3LIM(1,n)
6667 FORMAT ( 'I',22X, 3X,A2,'(',I2,')--', A2,'(',
* I2,')--',A2,'(',I2,')', F15.8, F11.3,
* i6,'-',i2,2F10.3, F12.4,16X, 'I')
if (i3BP(1,N).ne.i3BP(3,N)) then
R3LIM(2,n) = LOG(0.999999D0/0.000001) /
* R3BGRD(2,N) + R3BLIM(2,N)
if (runopt(8).eq.'BMH-EXP* ') then
R3LIM(2,n) = LOG(0.9999D0/0.0001D0) /
* R3BGRD(2,N) + R3BLIM(2,N)
end if
if (r3limax.lt.r3lim(2,n)) r3limax=r3lim(2,n) /////
WRITE (16,6668) i3bp(2,n),i3bp(3,n),
* R3BLIM(2,N),
* R3BGRD(2,N), R3LIM(2,n)
6668 FORMAT ( 'I',73X, i6,'-',i2,
* 2F10.3, F12.4,16X, 'I')
end if
END IF
140 CONTINUE
END IF
C
DO 250 I = 10, NRCUT(2)
RIJ = REAL(I) * 0.01
ARIJ = 1.0 / RIJ
DO 240 J = 1, NPAIR
E1(I,J) = 0.0
F1(I,J) = 0.0
E1M = 0.0
F1M = 0.0
IF (ABS(AIJ(J)).LT.1.0E-5) GO TO 220
EX = 0.0
IF (BIJ(J).GT.0.00001) THEN
ARB = (AIJ(J) - RIJ) / BIJ(J)
IF (ARB.GT.-128.0) EX = EXP(ARB)
END IF
E1(I,J) = BETA * BIJ(J)*EX*EPSij(J)
c * - CIJ(J)*ARIJ**6
C * - D4IJ(J)*ARIJ**4 - D7IJ(J)*ARIJ**7
F1(I,J) = BETA * EX*EPSij(J)
c * - 6.0*CIJ(J)*ARIJ**7
C * - 4.0*D4IJ(J)*ARIJ**5 - 7.0*D7IJ(J)*ARIJ**8
C * - 4.0*D4IJ(J)*ARIJ**5 - D4IJ(J)*ARIJ**4/4.43
C
220 AM1 = DM1IJ(J)*EXP(-BE1IJ(J)*RIJ)
AM2 = DM2IJ(J)*EXP(-BE2IJ(J)*RIJ)
am3 = dm3ij(j)*exp(-be3ij(j)*(rij-r03ij(j))**2)
E1M = BETA * (AM1 + AM2 + am3)
F1M = BETA * (BE1IJ(J)*AM1 + BE2IJ(J)*AM2 +
* 2.0*be3ij(j)*(rij-r03ij(j))*am3)
E1(I,J) = E1(I,J) + E1M ! BMH+VdW+EXP
F1(I,J) = F1(I,J) + F1M
IF (RSWTCH(J).LT.1.0E-6) THEN
E1(I,J) = E1(I,J) + E1M
F1(I,J) = F1(I,J) + F1M
ELSE IF (RIJ.LE.RSWTCH(J)) THEN
E1(I,J) = E1M ! R