PROGRAM MXDORTHO C=============================================================================== C## ## C## Program : MXDORTHO ## C## ## C## by Katsuyuki Kawamura (Univeristy of Tokyo) ## C## (Okayama University) ## C## (Hokkaido University) ## C## (Tokyo Institute of Technology) ## C## ## C## Configuration and Energy for Cubic and Non-Cubic Systems ## C## (Rectangular parallelepiped) ## C## with Pressure Control by stress tensor, ## C## and Quantum Correction for energy and pressure ## C## ## C## 2nd order interpolation of energy and force from tables ## C## ## C## First cubic version on Hitac 8800/8700(Univ.Tokyo) 1980 ## C## First orthogonal (crystal) version using CDC7600 1983-Oct ## C## at Manchester University ## C## HITAC M-280/IAP version (Univ. Tokyo) 1985-Sep-12 ## C## PX, PY, PZ pressure control version 1987-Feb-07 ## C## Pressure tensor and fractional coordinates 1987-Oct-29 ## C## Five elements and input data format and history 1987-Nov-05 ## C## PC9800RA+NDP-FORTRAN-386 (MS-DOS 32-bit PC) version 1989-Jan-26 ## C## Reviced for distribution by JCPE 1990-Apr-14 ## C## (XDORTO : DEFECT) 1990-Apr-21 ## C## 3-body interaction (H2O, Kumagai & Kats, 1994) 1991-Feb-02 ## C## Integrated version of MDorto and XDorto (MXD) 1991-May-22 ## C## Rearranged 1991-Oct-23 ## C## Seven comonents (atoms), rearranged 1992-Jan-23 ## C## Quantum corrections (Nakao & Kats) 1992-Mar-04 ## C## Ten comonents (atoms), rearranged 1992-Mar-31 ## C## Extended Andersen's pressure control (Katsuta & Kats) 1992-Apr-07 ## C## Metal (main group, Na) potential 1992-Apr-18 ## C## Revised for JCPE version 1992-Aug-01 ## C## 2nd order interpolation from E and F tables 1992-Sep-05 ## C## 2nd order interpolation of particle velocity 1992-Dec-12 ## C## Nose's thermostat 1992-Dec-14 ## C## Correction for trancation of V.der W-term 1993-Dec-10 ## C## Reviced 3-body by Kumagai 1994-Jan-30 ## C## L-J potential 1994-Jun-28 ## C## Nose's thermostat + quantum corrections 1994-Sep-01 ## C## Improvement of Semi-classical MD 1995-Jun-15 ## C## FILE09.DAT format changed 1995-Jul-18 ## C## IP model by Belonoshko & Dubrovinsky 1996-Sep-05 ## C## Electric Field (N.SAWAGUCHI) & Gravity Field 1997-Jun-30 ## C## Diatomic 3 chrge model (N2 and O2) 1997-Oct-20 ## C## 'ENERGY' and 'CUBE' options 1998-Aug-24 ## C## 'CONVEC' option for convection motion 1999-Feb-09 ## C## 'P ANDERS-C' for cubic Andersen 1999-Aug-23 ## C## Pair type potential model (PAIR-P) 1999-Sep-27 ## C## 3-body potential: j-i-k with j<>k 1999-Nov-16 ## C## 'EXCLUSION' : column and so on 2000-Apr-15 ## C## 3-body term sqrt(k1xk2) -> k1xk2 2000-May-01 ## C## Gradual Cell change with time 2000-May-28 ## C## POSISION-VELOCITY-ENERGY option (file09pv.dat) 2000-Dec-16 ## C## Soft repulsive wall in a basic cell 2001-Mar-07 ## C## Modify EWALD direct term 2001-Mar-24 ## C## 3-body term j-i-k : modified 2001-Sep-11 ## C## File07.dat : format (figures) 2001-Dec-02 ## C## Polyatomic molecules 2002-Feb-23 ## C## Modify NETWORK 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## Extended diatomic molecule (ion) 2004-Sep-26 ## C## Separate file08.dat (file081.dat) 2005-Aug-11 ## C## CUBE-F option (forced CUBE cell) 2005-Nov-07 ## C## temperature gradient in a basic cel 2007-Jul-31 ## C## file09v format change 2008-Nov-22 ## C## One- or Two- Dimensional change of cell 2008-Dec-15 ## C## file08 format changed (RDF) 2009-Feb-24 ## C## file09p and file09pv(pos) -> 5 figures 2009-Feb-25 ## C## Triatomic molecule (H2O, CO2, ...) 2010-May-26 ## C## Electric field at atoms 2010-Jul-09 ## C## Index of molecule at each atom in file07.dat 2010-Nov-21 ## C## Create or revice molecule table 2010-Nov-22 ## C## Limi of the number of 3-body term: 7 -> 17 2010-Dec-11 ## C==============================================================================| C Format and parameters of 'FILE05.DAT' file | C------------------------------------------------------------------------------| C 1 MD.......I....:....I....:....I....:....I....:....I....:....I....:....I | C XD.......I... : | C MDX......|... : | 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 SCALE-A: TMPGET : DELTMP : NTSTEP : TDUMP : [Scale each atom]| C T NOSE : TMPGET : DELTMP : STEMP : : : : | C T GRAD : : : STEMP : TDUMP : [Temperature grad]| 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 ANDERS-C SPRES(1): : :VIRM(1) : : : | 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 V CHANGE : ICAXIS : BTAGET : BCNGR(A par step) : : : | C 8 BUSING :MODE,MXN2: (ALPHA) : : : : : | C MORSE : : : : : : : | C MORSEQ : : : : : : : | C MORSE-AT : : : : : : : | C BMH-EXP : 3-body sqrt(k1xk2) : : : : | C BMH-EXP* : 3 body k1xk2 : : : : | C BELONO : : : : : : : | C TOSIFUMI : : : : : : : | C WOODCOCK : : : : : : : | C PAULING : : : : : : : | C METAL : : : : : : : | C PAIR-P : : : : : : : | C STSUNE : : : : : : : | C L-J : : : : : : : | C 81 N A NO. : ZI : WI : AI : BI : CI(VW) : DI() : | C - : : : : : [- not moved ]| C x,y,z : : : :[x,y,z fix x or y or z coordinate]| C * : : : : : [* dummy atoms ]| C = : : : : : [= Morse only ]| C / : : : : : [/ no T-control]| 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 82 I J : AIJ : BIJ : CIJ : : (eV) : [Pair-U]| C 82 I J : AIJ : BIJ : CIJ : (kJ/mol) : [Pair-P]| C 82e[BLANK] : : : : : : : | C : : : : : : : | C------------------------------------------------------------------------------| C 91 STRUCTURE: : : : [9:Show distance etc.] | C 92 NETWORK :NFCION(1):NFCION(2): : [10:Network structure analy.]| C : : : : NFCION(1) should be 2. | C : : : : NFCION(2) should be 0 or 3. | C 93 VELOCITY :NS09PV :PVMULT : : [11:Record particle velocity]| C POSITION :NS09PV :PVMULT : : : [....... position]| C ENERGY :NS09PV :PVMULT : : : [....... energy ]| C POSVELENE:NS09PV : : : : [..... pos,velo,ener]| C FORCE :NS09PV : : : : [.......... force]| C 94 QUANTUM : : : : : [12:Quantum correction]| C 95 PCF, RDF : ISTEP : Rend(A) : : :[13:Format of PCF table]| C*96 DIPOLE : : : : : [14:E(dipole moment)]| C 97 CENTER : : : : [15:Centering of atom cluster]| C CENTERING: iaxcen : : : : : : | C 98 NO(MV=0) : : : : [16:No correction for morment]| C AM(MV=0) : Iamv : Namv : : [Moment correction for Iamv]| C : : : if Namv>0 then oly Namv atoms used | 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 :Zmole2(1):DINTRA2(1):iatom2(1): : : icont : | C :Zmole2(2):Dintra2(2):iatom2(2): : : : | C : [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 9h CUBE : : : : : [24:to Cubic cell] | C CUBE-F : : : : : [forced CUBE] | C 9i CONVEC : FCONVC : : : : [25:Convection] | C 9j MOLECULE :dMOLintra:MOLstart : MOLend : : [26:Define molecule]| C 9k EXCLUSION: : : : : [27:Exclusion] | C COLUMN : iaex : Rexcl(radius) F : : (R>0 out) : | C SLAB : iaex : Rexcl(Thickness/2) F : : (R<0 in ) : | C CUBE : Rexcl(edge/2) : Fexcl : : : : | C SPHERE : Rexcl(radius) : Fexcl : : : : | C HONEYCOMB: iaex : Rexcl(radius) : Fexcl : : : : | C 9l WALL : A : B : : :[28:Soft repulsive wall]| C 9m POLYATOMS:dMOLintra:MOLstart : MOLend : :[29:Polyatomic molecule]| C 9n REMOVE : RMZL : RMZH : RMVZ : [30:Remove atom condition]| C 9o T GRAD : IAXTGR : T000 : T050 : [31:Temperature gradient] | C 9p CELL CHAN: 0./1. : 0./1. : 0./1. : [32:1- or 2- dimensional NPT]| C 9q MOLTABLE : : : : :[33:make molecule table]| C 9r ........ : : : : : : : | C 9s [BLANK] : : : : : : : | C 9 : : : : : : : | C MD.......I....:....I....:....I....:....I....:....I....:....I....:....:....| C REPEAT 1 TO 9 | C==============================================================================| C IRECRD(1-9) NRECRD(1-9) : 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 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,5X, 73X, ' I')
2223 FORMAT ('I ',I7,I5,I3,I7,4X, I7,I5,I3,I7,5X, I6,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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C
COMMON /COUNTS/ NJOB(2),IRECRD(9),NRECRD(9),IHISTR(4,111),PVMULT
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VALUES/ TVAL(LVA),TVALL(LVA),VALMAX(LVA),VAL0(LVA),
* SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* VAL(LVA),AVA(LVA,L50), NAV,NAVT
REAL *8 TVAL,SVAL,TVALL,SVALL,VALMAX,VALMIN,VAL0,VAL,AVA
COMMON /RADIAL/ NRDF(LTB,LEE),IPRDF(2)
INTEGER *4 NRDF
COMMON /GEOMET/ DTO(2),AVTHT(12),MBR(8,8,2),NRG(13,2),
* RTO(2),SVTHT(12),NBR(8,8,2),MEB(13,2),
* NTO(2),NVTHT(12),ANGL(3,12),TTAB(LST),NTT(121,12),
* ANCN(7,2),NTBL, ITBR(121,12)
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
REAL *8 H(3,3)
CHARACTER *10 RUNO18, RUNO19
CHARACTER *4 TITLE0(15), BIN
CHARACTER *1 DEFECT, ANS
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,
* 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 (6X,14('='),' Titles in FILE07.DAT and FILE05.DAT are ',
* 14('=') / '=====[F7]: ',15A4,' =====' /
* '=====[F5]: ',15A4,' =====' )
end if
C
CT ------------------- delete this block-if in case of oblique system
IF (BOX(4)**2+BOX(5)**2+BOX(6)**2.GT.1.E-6) THEN
WRITE (*,*) 'Error: The box shape is not suitable for ',
* 'MXDORTO !!!'
WRITE (*,1131) BOX(4),BOX(5),BOX(6)
1131 FORMAT (' BOX(4 to 6) are ',3F12.7)
WRITE (*,*) 'Is it posibble to change BOX(4), BOX(5), and',
* ' BOX(6) as zero ? (y/n)'
READ (5,1141) ANS
1141 FORMAT (A1)
IF (ANS.EQ.'n' .OR. ANS.EQ.'N') STOP
BOX(4) = 0.0
BOX(5) = 0.0
BOX(6) = 0.0
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),
* iioo, ixmole(i)
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
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.'
c
go to 201
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
C -------------------------------------- Input file of xtal geometry
NBOX(1) = 1
NBOX(2) = 1
NBOX(3) = 1
IF (RUNOPT(17).EQ.'CRYSTAL ') CALL FILE10
c ------------------------------------------------------------------
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,13)
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 = ' '
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,
* DTIME, RUNO18, BOX,
* DENSTY, RUNO19, VBOX
IF (RUNO18.EQ.'THERMOSTAT') WRITE (17,7080) STEMP,VSTEMP
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, ixmole(i)
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,13)
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 ',12('='),
* ' End=',I7,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,1x,i6)
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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VALUES/ TVAL(LVA),TVALL(LVA),VALMAX(LVA),VAL0(LVA),
* SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* VAL(LVA),AVA(LVA,L50), NAV,NAVT
REAL *8 TVAL,SVAL,TVALL,SVALL,VALMAX,VALMIN,VAL0,VAL,AVA
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 /WORK02/ IP(3,LNI), PP(3,LNI)
C
REAL *8 H(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.7 / 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, H
READ (19) ((PP(J,I),J=1,3),I=1,MMMMM)
WRITE (22) L, H
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, H
READ (22) ((PP(J,I),J=1,3),I=1,MMMMM)
WRITE (19) L, H
WRITE (19) ((PP(J,I),J=1,3),I=1,MMMMM)
450 CONTINUE
ELSE
DO 460 K = 1, NRECRD(4)
READ (19,9002) L, H
READ (19,9001) ((IP(J,I),J=1,3),I=1,MMMMM)
WRITE (22,9002) L, H
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, H
READ (22,9001) ((IP(J,I),J=1,3),I=1,MMMMM)
WRITE (19,9002) L, H
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 (I7,3X, 9F7.3)
END
C
C
C ========
C================================================================ FILE10
SUBROUTINE FILE10
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /BOXCNG/ BTAGET, BCNGR, ICAXIS
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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 /VECTOR/ FNV(LNV),UNV(LNV),PNV(6,LNV),ZIA(LEM),UCSELF,
* ALPHA, UCSELFI(LEM), MODE, NVN, NVEC(3,LNV)
REAL *8 FNV,UNV,PNV,ZIA,UCSELF,ALPHA,UCSELFI
COMMON /VALUES/ TVAL(LVA),TVALL(LVA),VALMAX(LVA),VAL0(LVA),
* SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* VAL(LVA),AVA(LVA,L50), NAV,NAVT
REAL *8 TVAL,SVAL,TVALL,SVALL,VALMAX,VALMIN,VAL0,VAL,AVA
COMMON /RADIAL/ NRDF(LTB,LEE),IPRDF(2)
INTEGER *4 NRDF
COMMON /GEOMET/ DTO(2),AVTHT(12),MBR(8,8,2),NRG(13,2),
* RTO(2),SVTHT(12),NBR(8,8,2),MEB(13,2),
* NTO(2),NVTHT(12),ANGL(3,12),TTAB(LST),NTT(121,12),
* ANCN(7,2),NTBL, ITBR(121,12)
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), fconvc, MEFD
REAL *8 EFD, EFREQ, GFD
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 /EXCLUS/ REXCL, Fexcl, iaex, iextype
common /WALLP/ WALLa, WALLb
common /REMOVE/ RMZL,RMZH,RMVZ
COMMON /WORK01/ VV(3,LNI),DUM(3,LNI)
COMMON /WORK02/ IPV(3,LNI),IDUMMY(3,LNI)
C
REAL *8 BOXA(6), FA(3)
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.' .or.
* RUNOP1.EQ.'mdx.') THEN
RUNOPT(1) = 'MD........'
RUNOP1 = 'MD..'
IP0 = 1
END IF
IF (RUNOP1.EQ.'MD..' .or.
* RUNOP1.EQ.'md..') THEN
RUNOPT(1) = 'MD........'
RUNOPT(17) = 'AMORPHOUS '
END IF
IF (RUNOP1.EQ.'XD..' .or.
* 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.'stop ' .OR.
* RUNOPT(2).EQ.'END ' ) GO TO 888
IF (RUNOPT(2).EQ.'CONT. ' .or.
* runopt(2).eq.'cont. ' .or.
* runopt(2).eq.'CONTIMUE ' .or.
* runopt(2).eq.'continue ')
* RUNOPT(2) = 'CONTINUE '
GO TO 50
C
888 INOEND = -1
RETURN
C
C -------------------------------- Read file07.dat, file08.dat, etc.
c Enter subroutine f07f08
50 CALL F07F08 (INOEND)
C ---------------------------------------------- Title on file06.dat
CALL TITLET (1,1)
C ------------------------------------------- Economy, normal detail
READ (15,1000) RUNOPT(3), AREC1, AREC2, AREC3, AREC4, AREC5
if (runopt(3).eq.'economy ') runopt(3)='ECONOMY '
if (runopt(3).eq.'normal ') runopt(3)='NORMAL '
if (runopt(3).eq.'detail ') runopt(3)='DETAIL '
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 ------------------------------------------------------ Temperature
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 NO-SCALE') RUNOPT(5) = 'T NO-CNTL '
IF (RUNRUN.EQ.'T SCALING ' .or.
* RUNRUN.EQ.'T scaling ') THEN
RUNOPT(5) = 'T SCALING '
NTSTEP = STEMP0
IF (NTSTEP.LE.0) NTSTEP = 10
END IF
IF (RUNRUN.EQ.'T SCALE-A ') THEN
RUNOPT(5) = 'T SCALE-A '
NTSTEP = STEMP0
IF (NTSTEP.LE.0) NTSTEP = 10
END IF
IF (RUNRUN.EQ.'T NOSE ' .or.
* runrun.eq.'T Nose ') RUNOPT(5) = 'T NOSE '
IF (RUNRUN.EQ.'T GRAD ' .or.
* runrun.eq.'T grad ') RUNOPT(5) = 'T GRAD '
IF (NTSTEP.LE.0) NTSTEP = 1
DELTMP = DELT
TMPGET = TARGT
IF (TDUMP.LE.0.0001) TDUMP = 0.5
IF (RUNOPT(5) .NE.'T NOSE ' .OR.
* RUNOPT(2) .NE.'CONTINUE ' .OR.
* RUNOPT(51).NE.'THERMOSTAT' ) THEN
STEMP = STEMP0
VSTEMP = 0.0
END IF
C --------------------------------------------------------- Pressure
READ (15,1000) RUNRUN, SPRES, VIRM(1), VIRM(2), VIRM(3)
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 ' .or.
* runrun.eq.'P scaling ') then
RUNOPT(6) = 'P SCALING '
pdump = virm(1)
if (pdump.lt.0.001) pdump = 1.0
end if
IF (RUNRUN.EQ.'P ANDERSEN' .OR.
* runrun.eq.'P Andersen' .OR.
* RUNRUN.EQ.'P ANDERS-C' ) THEN
if (RUNRUN.EQ.'P ANDERSEN') RUNOPT(6) = 'P ANDERSEN'
if (RUNRUN.EQ.'P Andersen') RUNOPT(6) = 'P ANDERSEN'
if (RUNRUN.EQ.'P ANDERS-C') RUNOPT(6) = 'P ANDERS-C'
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
C ----------------------------------------------------------- Volume
READ (15,1000) RUNRUN, BOXA
IF (RUNRUN.EQ.' ') RUNOPT(7) = 'V FREE '
IF (RUNRUN.EQ.'V CONST. ' .or.
* runrun.eq.'V const. ' .or.
* runrun.eq.'V CONSTANT' .or.
* runrun.eq.'V constant') 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 ' .or.
* 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
C ----------------------------------------- Change density
ELSE IF (RUNRUN.EQ.'V DENSITY ' .or.
* 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
C ---------------------------------------- Uniaxizl change
ELSE IF (RUNRUN.EQ.'V CHANGE ' .or.
* runrun.eq.'V change ') THEN
RUNOPT(7) = 'V CHANGE '
ICAXIS = BOXA(1)
BTAGET = BOXA(2)
BCNGR = BOXA(3)
if (ABS(BCNGR).le.1.0E-6)
* BCNGR = sign(1.0,BCNGR)*1.0E-6
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-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.'PAIR-P ' .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)
IION(I) = 0
IF (I.NE.1) ZSUM = ZSUM + ZJ * ANJ
IF (ATY.EQ.'-') IION(I) = -1 ! P-fixed
IF (ATY.eq.'x') IION(I) = -11 ! fix x-coordinate
IF (ATY.eq.'y') IION(I) = -12 ! fix y-coordinate
IF (ATY.eq.'z') IION(I) = -13 ! fix z-coordinate
IF (ATY.EQ.'*') IION(I) = -999 ! dummy atom
IF (ATY.EQ.'=') IION(I) = 1 ! Morse only
IF (ATY.EQ.'/') IION(i) = 2 ! no T-control
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
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 (RUNOP1.EQ.'MD..'.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
IF (RUNOPT(6).NE.'P NO-CNTL ') THEN
IF (RUNOPT(6).EQ.'P SCALING ') WRITE (16,2020) RUNOPT(6),SPRES
IF (RUNOPT(6).EQ.'P ANDERSEN') WRITE (16,2027) RUNOPT(6),
* SPRES,(VIRM(LL),LL=1,3)
IF (RUNOPT(6).EQ.'P ANDERS-C') WRITE (16,2027) RUNOPT(6),
* SPRES,(VIRM(LL),LL=1,3)
END IF
if (RUNOPT(7).NE.'V '.and.RUNOPT(7).NE.'V CONST ') then
if (RUNOPT(7).eq.'V CHANGE ') write (16,2031) runopt(7),
* ICAXIS, BTAGET, BCNGR
end if
C
CALL TABLER (1)
C
C ------------------------------------------------- Read RUNOPT(9),...,(33)
write (16,2040)
lentab = lst
IPRDF(1) = 2
IPRDF(2) = 9999
520 READ (15,1000) RUNRUN,PARAM1,PARAM2,PARAM3,PARAM4,PARAM5,PARAM6
IF (RUNRUN.NE.' ') THEN
write (6,*) runrun
IF (RUNRUN.EQ.'STRUCTURE ' .or. ! Structure =================
* RUNRUN.EQ.'structure ') then
RUNOPT(9) = 'STRUCTURE '
lentab = param1
if (lentab.lt.1) lentab = LST
if (lentab.gt.LST) lentab = LST
end if
IF (RUNRUN.EQ.'NETWORK ' .or. ! Network ===================
* runrun.eq.'network ') THEN
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 (*,*) 'Network forming cation(s) is(are)',
* (i,atmnet(i),i=1,natx)
END IF
C
IF (RUNRUN.EQ.'VELOCITY ' .or. ! Velocity ===================
* runrun.eq.'velocity ') THEN
RUNOPT(11) = 'VELOCITY '
IRECRD(9) = PARAM1
PVMULT = 50000.0
IF (PARAM2.GT.0) PVMULT = PARAM2
IF (IRECRD(9).LE.0) IRECRD(9) = 1
END IF
IF (RUNRUN.EQ.'POSITION ' .or. ! Position =================
* runrun.eq.'position ') THEN
RUNOPT(11) = 'POSITION '
IRECRD(9) = PARAM1
PVMULT = 90000.0
IF (PARAM2.GT.0) PVMULT = PARAM2
IF (IRECRD(9).LE.1) IRECRD(9) = 1
END IF
IF (RUNRUN.EQ.'ENERGY ' .or. ! Energy ===================
* 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.'FORCE ' .or. ! Force =====================
* runrun.eq.'force ') THEN
RUNOPT(11) = 'FORCE '
IRECRD(9) = PARAM1
PVMULT = 1.0E0
IF (PARAM2.GT.0) PVMULT = PARAM2
IF (IRECRD(9).LE.1) IRECRD(9) = 1
END IF
IF (RUNRUN.EQ.'POSVELENE ' .or. ! POS.Vel.Ene. ==============
* 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 ' .or. ! Quantum correction =======
* runrun.eq.'quantum ') THEN
RUNOPT(12) = 'QUANTUM '
CALL QCTABL
END IF
IF (RUNRUN.EQ.'PCF '.OR. ! PCF table =================
* RUNRUN.EQ.'RDF '.or.
* runrun.eq.'pcf '.or.
* 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 ' .or. ! Dipole ===================
* runrun.eq.'dipole ') THEN
RUNOPT(14) = 'DIPOLE '
END IF
IF (RUNRUN.EQ.'CENTER '.OR. ! Center ===================
* RUNRUN.EQ.'CENTRE '.or.
* runrun.eq.'center '.or.
* runrun.eq.'centre ') THEN
RUNOPT(15) = 'CENTER '
END IF
IF (RUNRUN.EQ.'CENTERING ' .or. ! Centering ================
* runrun.eq.'centering ') THEN
RUNOPT(15) = 'CENTERING '
iaxcen = PARAM1
END IF
IF (RUNRUN.EQ.'NO(MV=0) ') THEN ! No sigma(mv)=0 =========
RUNOPT(16) = 'NO(MV=0) '
END IF
IF (RUNRUN.EQ.'AM(MV=0) ') THEN ! Atom sp sigma(MV)=0 ====
RUNOPT(16) = 'AM(MV=0) '
Iamv = param1
Namv = param2
if (Namv.gt.nion(Iamv)) Namv= nion(Iamv)
if (Namv.le.0) Namv = nion(Iamv)
END IF
IF (RUNRUN.EQ.'CRYSTAL ') THEN ! Crystal ================
RUNOPT(17) = 'CRYSTAL '
END IF
IF (RUNRUN.EQ.'BINARY ') THEN ! Binary output ==========
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 tensol ========
RUNOPT(19) = 'PRESSURE '
OPEN (27, FILE=FLNAME(13), STATUS='UNKNOWN',
* ACCESS='SEQUENTIAL', FORM='FORMATTED' )
REWIND 27
END IF
IF (RUNRUN.EQ.'ELEC.FIELD') THEN ! Electric field ==========
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 ! Gravity field ===========
runopt(21) = 'GRAV.FIELD'
gfd(1) = param1
gfd(2) = param2
gfd(3) = param3
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 ! Toriatomic 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
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.'CUBE ') then ! Cube basic cell ===========
runopt(24) = 'CUBE '
end if
if (runrun.eq.'CUBE-F ') then ! Forced Cube basic cell ====
runopt(24) = 'CUBE-F '
end if
if (runrun.eq.'CONVEC ') then ! Convection ================
runopt(25) = 'CONVECTION'
fconvc = param1
write (6,*) '[CONVECTION] option is set'
end if
if (runrun.eq.'MOLECULE ') then ! Molecule ==================
runopt(26) = 'MOLECULE '
dMOLintra = param1
MOLstart = param2
MOLend = param3
call MOLECULE
end if
if (runrun.eq.'EXCLUSION ') then ! Exclusion ===============
runopt(27) = 'EXCLUSION '
READ (15,1000) RUNRUN,PARAM1,PARAM2,PARAM3,PARAM4,
* PARAM5,PARAM6
if (RUNRUN.eq.'COLUMN '.or.
* RUNRUN.eq.'SLAB ' ) then
iextype = 1 !---------------- COLUMN
if (RUNRUN.eq.'SLUB ') iextype = 2 !-------- SLAB
iaex = param1
REXCL = param2
Fexcl = param3
c write (6,*) iextype, iaex,rexcl,fexcl
end if
if (RUNRUN.eq.'CUBE ') then
iextype = 3 ! ----------------- CUBE
rexcl = param1
Fexcl = param2
end if
if (RUNRUN.eq.'SPHERE ') then
iextype = 4 ! --------------- SPHERE
REXCL = param1
Fexcl = param2
end if
if (RUNRUN.eq.'HONEYCOMB ') then
iextype = 5 ! ------------ HONEYCOMB
iaex = param1
rexcl = param2
fexcl = param3
end if
if (Fexcl.lt.1.0E-9) Fexcl = 1.0E-5
end if
if (runrun.eq.'WALL ') then ! Wall in cell =============
runopt(28) = 'WALL '
WALLa = param1
WALLb = param2
end if
if (runrun.eq.'POLYATOMS ') then ! PolyAtomic molecule ======
runopt(29) = 'POLYATOMS '
dMOLintra = param1
MOLstart = param2
MOLend = param3
call MOLECULE
end if
if (runrun.eq.'REMOVE ') then ! Remove atom(s) ===========
runopt(30) = 'REMOVE '
RMZL = param1
RMZH = param2
RMVZ = param3
end if
if (runrun.eq.'T GRAD ') then ! Temperature gradient =====
runopt(31) = 'T GRAD '
IAXTDR = param1
T000 = param2
T050 = param3
end if
if (runrun.eq.'CELL CHAN ') then ! Cell size change with time ==
runopt(32) = 'CELL CHAN '
ICFIX(1) = param1+0.00001
ICFIX(2) = param2+0.00001
ICFIX(3) = param3+0.00001
end if
if (runrun.eq.'MOLTABLE ') then ! Create.revice molecule table =
runopt(32) = 'MOLTABLE '
end if
GOTO 520
END IF
WRITE (16,2030) (I,RUNOPT(I),I=1,40)
c
c ------------------------------ End of single job data read from file05.dat
c
C ------------------------------------------------------------ Check P and V
CALL CHECKP (DTMO)
C -------------------------------------------------------------- file09p.dat
IF (RUNOPT(2).EQ.'START ') THEN
IF (RUNOP1.EQ.'MD..') THEN
IF (TITLE(1).NE.'BENC' .OR.
* TITLE(2).NE. 'HMAR' ) THEN
NRECRD(4) = 1
IF (RUNOPT(18).EQ.'BINARY ') THEN
WRITE (19) NRECRD(4), 0, BOX(1), 0.0, 0.0,
* 0.0, BOX(2), 0.0, 0.0, 0.0, BOX(3)
WRITE (19) ((SNGL(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,9001) NRECRD(4), 0, BOX(1),
* 0.0, 0.0, 0.0, BOX(2),
* 0.0, 0.0, 0.0, BOX(3)
WRITE (19,9002) ((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 +50000.0
560 CONTINUE
WRITE (28,9001) NRECRD(1),IRECRD(9)
WRITE (28,9002) ((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),
* BOX(1), 0.0, 0.0, 0.0, BOX(2),
* 0.0, 0.0, 0.0, BOX(3)
WRITE (28) ((SNGL(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,9001) NRECRD(1),IRECRD(9),
* BOX(1), 0.0, 0.0, 0.0, BOX(2),
* 0.0, 0.0, 0.0, BOX(3)
WRITE (28,9002) ((IPV(J,I),J=1,3),I=1,NTION)
END IF
END IF
9001 FORMAT (I7,i3,9F7.3)
9002 FORMAT (18I5)
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' )
2020 FORMAT ('I [ ',A10,' ] Pressure is controlled at ',3F9.4,
* 'GPa using forced scaling of cell dimensions.',14X,
* 'I')
2027 FORMAT ('I [ ',A10,' ] Pressure is controlled at ',3F9.4,
* ' GPa by Andersen''s mass ',3(1X,G9.2E3),
* ' g I')
2031 format ('I [ ',A10,' ] Cell size of axis ',i1, ' is canged to ',
* F10.5, ' Angstroms by rate of ',F8.5,
* ' Angstroms/step',24x,'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' )
2040 format ('I',130('-'),'I' )
END
C
C
C ==========
C============================================================== MOLECULE
SUBROUTINE MOLECULE
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C ======================================recognize diatomic molecules
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /CARTES/ H(3,3),HINV(3,3),FTOQ(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
CT * ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,FTOQ,RBOX,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 rx, ry, rz, dx, dy, dz
integer mi(lni), ndistr(48)
c
cut2 = dMOLintra**2
do 10 I = 1, ntion
mi(i) = 0
10 continue
do 20 n = 1, 48
ndistr(n) = 0
20 continue
!
if (dMOLintra.lt.0) then !--------------- decomposition of ixmole
nmole=0
do 100 i=1, ntion
mm=mod(ixmole(i),900000)
mi(i)=mm
mmole(mm)=mmole(mm)+1
IMOLE(mmole(mm),mm) = i
if (mm.gt.nmole) nmole=mm
100 continue
do 110 n=1, nmole
ndistr(mmole(n))=ndistr(mmole(n))+1
110 continue
go to 800
end if
!
nnn = 1 ! No. of molecules
imole(1,nnn) = ions(1,MOLstart) ! 1st atom of 1st molecule
mi(ions(1,MOLstart)) = 1
mmole(nnn) = 1 ! No. of atoms in the molecule
!------------------------------------------- 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
!
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
! --------- delete these if-statements for triclinic
! IF (ABS(RX).GT.0.5) RX = RX - SIGN(1.0D0,RX)
! IF (ABS(RY).GT.0.5) RY = RY - SIGN(1.0D0,RY)
! 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
DX = RX * BOX(1)
DY = RY * BOX(2)
DZ = RZ * BOX(3)
RIJ2 = DX*DX + DY*DY + DZ*DZ
IF (RIJ2.le.CUT2) then
mmole(n) = mmole(n) + 1
IMOLE(mmole(n),n) = i
mi(i) = 1
go to 510
end if
400 CONTINUE
500 continue
nnn=nnn+1
imole(1,nnn) = i
mi(i)=1
mmole(nnn) = 1
510 CONTINUE
590 continue
!
write (6,*) 'nnn=',nnn
write (6,9999) (mmole(n),n=1,nnn) !!!!!!!!
9999 format (20I4)
!
do 660 n2=2, nnn ! nnn : number of molecules
mm2=mmole(n2) ! mm2 : number of atoms in n2-th molecule
if (mm2.le.0) go to 660
do 650 n1 = 1, n2-1 ! molecules n1-n2
mm1=mmole(n1) ! number of atoms in n1-th molecule
mm2=mmole(n2) ! number of atoms in n2-th molecule
if (mm1.le.0) go to 650
do 630 m1=1, mm1
do 640 m2=1, mm2
i=imole(m1,n1) ! atom i=m1 in n1
j=imole(m2,n2) ! atom j=m2 in 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
! --------- delete these if-statements for triclinic
! IF (ABS(RX).GT.0.5) RX = RX - SIGN(1.0D0,RX)
! IF (ABS(RY).GT.0.5) RY = RY - SIGN(1.0D0,RY)
! 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
DX = RX * BOX(1)
DY = RY * BOX(2)
DZ = RZ * BOX(3)
RIJ2 = DX*DX + DY*DY + DZ*DZ ! distance**2 i-j
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
!
!
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
do 770 n = 1, nmole
mm = mmole(n)
do 760 i=1, mm
j=imole(i,n)
ixmole(j)=900000 + n
760 continue
770 continue
c
800 write (6,1001) nmole
1001 format (' Total number of molecules is',I6)
c write (6,1002) (n,n=1,30), (ndistr(n),n=1,40)
write (6,1003) (ndistr(n),n,n=1,40)
1002 format ('N.A',15I5 / 3X,15I5 / 'N.M',15I5 / 3x,15I5)
1003 format (8(I5,'[',I2,'] '))
c
RETURN
END
C
C
C ========
C================================================================ DIATOM
SUBROUTINE DIATOM
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C ======================================recognize diatomic molecules
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /CARTES/ H(3,3),HINV(3,3),FTOQ(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
CT * ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,FTOQ,RBOX,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)
c DX = H(1,1)*RX + H(1,2)*RY + H(1,3)*RZ
c DY = H(2,1)*RX + H(2,2)*RY + H(2,3)*RZ
c DZ = H(3,1)*RX + H(3,2)*RY + H(3,3)*RZ
DX = RX * BOX(1)
DY = RY * BOX(2)
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
c
do 770 n = 1, nmole2
do 760 i=1, 2
j=idmole2(i,n)
ixmole(j)=200000 + n
760 continue
770 continue
c
RETURN
END
C ==========
C=============================================================== TRIATOM
SUBROUTINE TRIATOM
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C =====================================recognize triatomic molecules
c H2O, CO2, ...
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /CARTES/ H(3,3),HINV(3,3),FTOQ(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
CT * ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,FTOQ,RBOX,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
c
C --------------------------------------------calc distance of atoms
nnn = 0
do 900 iii = 1, 2
cut2 = dintra3(iii)**2
if (cut2.le.0.1) goto 900
io = iatom3(iii,1) ! Central atom of 3 atom molecule
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)
c DX = H(1,1)*RX + H(1,2)*RY + H(1,3)*RZ
c DY = H(2,1)*RX + H(2,2)*RY + H(2,3)*RZ
c DZ = H(3,1)*RX + H(3,2)*RY + H(3,3)*RZ
DX = RX * BOX(1)
DY = RY * BOX(2)
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
c
do 770 n = 1, nmole3
do 760 i=1, 3
j=idmole3(i,n)
ixmole(j)=300000 + n
760 continue
770 continue
RETURN
END
C
C
C ========
C================================================================ PREPAR
SUBROUTINE PREPAR (FORMUL)
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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/ TVAL(LVA),TVALL(LVA),VALMAX(LVA),VAL0(LVA),
* SVAL(LVA),SVALL(LVA),VALMIN(LVA),
* VAL(LVA),AVA(LVA,L50), NAV,NAVT
REAL *8 TVAL,SVAL,TVALL,SVALL,VALMAX,VALMIN,VAL0,VAL,AVA
COMMON /RADIAL/ NRDF(LTB,LEE),IPRDF(2)
INTEGER *4 NRDF
COMMON /GEOMET/ DTO(2),AVTHT(12),MBR(8,8,2),NRG(13,2),
* RTO(2),SVTHT(12),NBR(8,8,2),MEB(13,2),
* NTO(2),NVTHT(12),ANGL(3,12),TTAB(LST),NTT(121,12),
* ANCN(7,2),NTBL, ITBR(121,12)
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, NCOMPO
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) NELM, ntion
4455 FORMAT (' ***** THE NUMBER OF PARTICLES IN FILE05 IS MORE THAN ',
* 'THAT IN FILE07 *****'/20x,'file05:',I6 / 20x,
* 'file07:',i6 )
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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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 RL ,TT,FV,DL,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.0) 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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VECTOR/ FNV(LNV),UNV(LNV),PNV(6,LNV),ZIA(LEM),UCSELF,
* ALPHA,UCSELFI(LEM), MODE, NVN, NVEC(3,LNV)
REAL *8 FNV,UNV,PNV,ZIA,UCSELF,ALPHA,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 /
*' ******* ************************** ',
*' **** *********** ******** ',
*' ***** ********* ******** ',
*' ****** ********** ********* ',
*' ******* *********** *********',
*' **** *** ************ *********',
*' *** *** *** ********* *********',
*' *** *** *** ********* *********',
*' *** *** *** ********* *********',
*' *** *** *** ********* ******** ',
*' *** ******* ********* ******* ',
*' **** ***** ********* ******* ',
*' ***** *** ********* ******* ',
*' ***** * ********* ******* ',
*' ******* ********* ****** ',
*' ******** *********** ****** ',
*'*********** ************************ R',
*' '/
DATA LOGO2 /
*'************ ************************* ',
*' ********* ************ ******* ',
*' ******** *********** ******* ',
*' ******* *** ******** ******** ',
*' ****** *** ******** ******** ',
*' ****** *** ******** ********',
*' ****** *** ******** ********',
*' ******** ******** ********',
*' ****** ******** ********',
*' ******** ******** ******* ',
*' *** ****** ******** ******* ',
*' *** ****** ******** ******* ',
*' *** ****** ******** ******* ',
*' *** ****** ******** ****** ',
*' **** ****** ******** ****** ',
*' ****** ******* ********** ****** ',
*'********** *************************** 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
CALL TMATRX
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, LOGO1(4), LOGO1(5)
WRITE (16,5002) RUNOPT(8),MODE,NVN, LOGO1(6),
* 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',29X,'I ',A63,
* ' I' /
* 'I--',60('-'),'I ',A63, ' I' )
5002 FORMAT('I ',A8,' I Mode =',I3, 13X, 'No.of Nv=',I5,
* 9X,'I ',A63,' I' /
* 'I ',8X,' I Alpha=',F6.3,' A-1 Rcut(L) =',
* F7.3,' A', 5X,'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 ',A50,A13,' 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.'MORSE-AT ') CALL MORSEP
if (runopt(8).eq.'MORSEQ ') CALL MORSEQ
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
IF (RUNOPT(8).EQ.'PAIR-P ') CALL PAIRP
C
IF (RUNOPT(3).EQ.'DETAIL ') THEN
DO 200 I = 10, 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 = 10, 300, 10
RIJ = I * 0.01
WRITE (16,6666) RIJ,F0(I),(F1(I,J),J=1,NPAIR)
210 CONTINUE
WRITE (16,6666)
DO 220 I = 10, 300, 10
RIJ = I * 0.01
WRITE (16,6666) RIJ,F0(I),
* (F1(I,J)+zij(j)*F0(i),J=1,NPAIR)
220 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-AT ' .OR. RUNOPT(8).EQ.'BMH-EXP ' .OR.
* runopt(8).eq.'MORSEQ ' .or. RUNOPT(8).EQ.'BMH-EXP* ' .OR.
* runopt(8).eq.'BMH-EXPQ ' .or.
* RUNOPT(8).EQ.'BELONO ' .OR. RUNOPT(8).EQ.'PAIR-P ' .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
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C
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,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /CARTES/ H(3,3),HINV(3,3),FTOQ(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
CT * ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,FTOQ,RBOX,TRANSX,TRANSY,TRANSZ
C
REAL *8 SINA(3), COSA(3), DET, GG
C ---------------------------- cos and sin of alpha, beta, and gamma
DO 20 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)
20 CONTINUE
C
C ------------------ Transformation matrix from crystal to Cartesian
C
H(1,3) = 0.D0
H(2,3) = 0.D0
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)
H(2,1) = BOX(1)*COSA(3)*SINA(1)
H(1,1) = BOX(1)*SQRT(1-COSA(2)**2-(COSA(3)*SINA(1))**2)
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
DENSTY = TWEGHT / (ANA * VOL * 1.0D-24)
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
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
C ---------------------------------------------------- Metric tensor
DO 80 I = 1, 3
DO 80 J = 1, 3
GG = 0.0D0
DO 70 K = 1, 3
GG = GG + H(K,J) * H(K,I)
70 CONTINUE
G(J,I) = GG
80 CONTINUE
CALL INVERS (G, DET, GINV)
C -------------------------- Trans. of reciprocal force to cartesian
C
FTOQ(1,1) = H(1,1) / BOX(1)
FTOQ(2,1) = H(2,1) / BOX(1)
FTOQ(3,1) = H(3,1) / BOX(1)
FTOQ(1,2) = H(1,2) / BOX(2)
FTOQ(2,2) = H(2,2) / BOX(2)
FTOQ(2,3) = H(3,2) / BOX(2)
FTOQ(1,2) = H(1,3) / BOX(3)
FTOQ(2,2) = H(2,3) / BOX(3)
FTOQ(2,3) = H(3,3) / BOX(3)
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)
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 110 I = 0, 1
DO 110 J = 0, 1
DO 110 K = 0, 1
N = N + 1
TRANSX(N) = I
TRANSY(N) = J
TRANSZ(N) = K
110 CONTINUE
RETURN
END
C
C
C ========
C================================================================ INVERS
SUBROUTINE INVERS (X, DET, XINV)
C -------------------------------------------- Given 3 by 3 matrix X
C Store determinant at DET 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
c SUBROUTINE PTOXYZ (I)
c PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
c * LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
c * LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
c * LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C
c COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
c * UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
c * NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
c * NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
c REAL *8 P,V,VP,P0,UI,AU,AV3BP
c COMMON /CARTES/ H(3,3),HINV(3,3),FTOQ(3,3),RBOX(6),
c * G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
CT * ,Q(3,LNI),Q0(3,LNI)
c REAL *8 H,HINV,G,GINV,FTOQ,RBOX,TRANSX,TRANSY,TRANSZ
C
CT REAL *8 PX,PY,PZ
C
C -------------------------------- TRANSFORMATION OF ION COORDINATES
C FROM CRYSTAL TO CARTESIAN (X,Y,Z)
C
CT PX = P(1,I)
CT PY = P(2,I)
CT PZ = P(3,I)
CT Q(1,I) = H(1,1)*PX + H(1,2)*PY + H(1,3)*PZ
CT Q(2,I) = H(2,1)*PX + H(2,2)*PY + H(2,3)*PZ
CT Q(3,I) = H(3,1)*PX + H(3,2)*PY + H(3,3)*PZ
C
CT PX = P0(1,I)
CT PY = P0(2,I)
CT PZ = P0(3,I)
CT Q0(1,I) = H(1,1)*PX + H(1,2)*PY + H(1,3)*PZ
CT Q0(2,I) = H(2,1)*PX + H(2,2)*PY + H(2,3)*PZ
CT Q0(3,I) = H(3,1)*PX + H(3,2)*PY + H(3,3)*PZ
C RETURN
c END
C
C
C ========
C================================================================ XYZTOP
SUBROUTINE XYZTOP
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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),FTOQ(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
CT * ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,FTOQ,RBOX,TRANSX,TRANSY,TRANSZ
C
CT 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
C QX = Q(1,I)
C QY = Q(2,I)
C QZ = Q(3,I)
C P(1,I) = HINV(1,1)*QX + HINV(1,2)*QY + HINV(1,3)*QZ
C P(2,I) = HINV(2,1)*QX + HINV(2,2)*QY + HINV(2,3)*QZ
C P(3,I) = HINV(3,1)*QX + HINV(3,2)*QY + HINV(3,3)*QZ
C
C QX = Q0(1,I)
C QY = Q0(2,I)
C QZ = Q0(3,I)
C P0(1,I) = HINV(1,1)*QX + HINV(1,2)*QY + HINV(1,3)*QZ
C P0(2,I) = HINV(2,1)*QX + HINV(2,2)*QY + HINV(2,3)*QZ
C 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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* NTIOND,NIOND(LEM), IOND(LNI), NPAIR, IION(LEM)
REAL *8 P,V,VP,P0,UI,AU,AV3BP
COMMON /VECTOR/ FNV(LNV),UNV(LNV),PNV(6,LNV),ZIA(LEM),UCSELF,
* ALPHA,UCSELFI(LEM), MODE, NVN, NVEC(3,LNV)
REAL *8 FNV,UNV,PNV,ZIA,UCSELF,ALPHA,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),FTOQ(3,3),RBOX(6),
* G(3,3),GINV(3,3),TRANSX(8),TRANSY(8),TRANSZ(8)
CT * ,Q(3,LNI),Q0(3,LNI)
REAL *8 H,HINV,G,GINV,FTOQ,RBOX,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, PAA2,ELC2,ASP,ERFC,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 = 1, 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
I = IL + 1 - II
XN = I * DBLE(RBOX(1))
DO 260 JJ = 1, JL2
J = JL + 1 - JJ
YN = J * DBLE(RBOX(2))
DO 250 KK = 1, KL2
K = KK - 1
ZN = K * DBLE(RBOX(3))
IF (K.GT.0) GO TO 230
IF (J.LT.0) GO TO 250
IF (J.EQ.0 .AND. I.LE.0) GO TO 250
230 VN2 = XN**2 + YN**2 + ZN**2 +
* 2*(XN*YN*RBOX(6) + YN*ZN*RBOX(4) +
* XN*ZN*RBOX(5))
IF (VN2.GT.ABC2) GO TO 250
NVN = NVN + 1
IF (NVN.GT.LNV) THEN
WRITE (*,9901) ABS(MODE),lnv
9901 FORMAT (' ******* SET [MODE] LESS THAN ',I2,
* ' (LNV=',i5,') *******')
STOP
END IF
NVEC(1,NVN) = I
NVEC(2,NVN) = J
NVEC(3,NVN) = K
EXPVN = EXP(- VN2 * PIAL2) / VN2
FNV(NVN) = FCT * EXPVN
UNV(NVN) = UCT * EXPVN
PAA2 = 2.0D0 * (PIAL2 + 1.0D0/VN2)
PNV(1,NVN) = PCT * (1.0D0 - PAA2 * XN**2) * EXPVN
PNV(2,NVN) = PCT * (1.0D0 - PAA2 * YN**2) * EXPVN
PNV(3,NVN) = PCT * (1.0D0 - PAA2 * ZN**2) * EXPVN
PNV(4,NVN) = PCT * (0.0D0 - PAA2 * XN*YN) * EXPVN
PNV(5,NVN) = PCT * (0.0D0 - PAA2 * XN*ZN) * EXPVN
PNV(6,NVN) = PCT * (0.0D0 - PAA2 * YN*ZN) * EXPVN
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) = ZIO(IO)*ZIO(IO)*ASP*2.0
310 CONTINUE
RETURN
END
C
C
C ========
C================================================================ VWCORR
SUBROUTINE VWCORR
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
C
C --------- Correction of energy and pressur for Van der Waals terms
C
COMMON /TEMPRS/ DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
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 /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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 /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.0D0
VCORR = 0.0D0
N = 0
DO 230 I = 1, NCOMPO
DO 220 J = 1, I
N = N + 1
SATOMS = NION(I) * NION(J) / VOL * PI4
C SATOMS = NION(I) * NION(J) / VOL * PI4 * BETA
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 (*,*) RCUT(2), RCUT(1)
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================================================================ MORSEP
SUBROUTINE MORSEP
PARAMETER (LNI=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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 EALPHA, BETA, RIJ,ARIJ, E1M,F1M, AM1,AM2,
* EX, ARB, ZFORML(LEM), epsij(lef), 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
N = N + 1
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.'BELONO ' ) 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
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
EALPHA = DIJ(J)*ARIJ**4*EXP(-RIJ/4.43)*1.6022E-12
E1(I,J) = BETA * BIJ(J)*EX*EPSIJ(J)
C * - CIJ(J)*ARIJ**6 )
C * + EALPHA
F1(I,J) = BETA * EX*EPSIJ(J)
C * - 6.0*CIJ(J)*ARIJ**7)
C * + 4.0*EALPHA*ARIJ + EALPHA/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 -
* 2.0*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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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 BETA, RIJ,ARIJ, E1M,F1M, AM1,AM2,
* EX, ZFORML(LEM), epsij(lef), 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=62387,LTB=10004, LEL=8, LEM=10, LCT=5000000,
* LSR=1254, LEE=LEL*(LEL+1)/2, L50=LCT/50+1,
* LAA= 512, LNV=29876, LEF=LEM*(LEM+1)/2, LST=32,
* LAT=LAA*4,LVA=24+LEM*2, L3P=17, LRG=LNI*5 )
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,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES(3),PPXYZ(7),FJMOL,
* T000,T050, IAXTGR, NTSTEP
REAL *8 DTIME,TEMP,DELTMP,TMPGET,TINT,TPRE,STEMP,VSTEMP,
* TDUMP,PDUMP,SPRES,PPXYZ,FJMOL
COMMON /ABOXOF/ BOX(6),VBOX(6),VOL,DENSTY,VIRM(6),
* RCUT(2),NRCUT(2),MXCUT, NFORML,IAXCEN,ICFIX(3)
REAL *8 BOX, VBOX, VOL, DENSTY, VIRM, RCUT
COMMON /ATOMSI/ P(3,LNI), V(3,LNI), VP(3,LNI), P0(3,LNI),
* UI(LNI), AU(LNI), AV3BP(2,L3P), ixmole(LNI),
* NTION, NION(LEM), IONS(2,LEM),NCOMPO, Iamv,Namv,
* 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 EALPHA, BETA, RIJ,ARIJ, E1M,F1M, AM1,AM2,
* EX, ARB, epsij(lef), sepij(lef)
real *8 am3, dm3ij(lef), be3ij(lef), r03ij(lef)
integer ipara(2,10), npara
real *4 apara(8,10)
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
N = N + 1
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
ZIJ(N) = ZIO(I)*ZIO(J)
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)
c write (6,*) IP,JP, KP, ijkl,
c * D1, BE1, D2, BE2, RSIJP, GGG
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 ! 3-body term
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 !------------------ j-i-j
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 !------------------- J-i-k
N3BP = N3BP +1
c write (6,*) ip,jp,kp
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 ',24x,'DM1ij BE1ij DM2ij ',
* ' BE2ij DM3ij BE3ij R03ij ',
* ' Rswch',26x, '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 ',A2,'(',i2,') -- ',A2,'(',i2,') ',
* 3(F11.2, F10.3),F10.3,F10.3, 26X,'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 ',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)
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
EALPHA = DIJ(J)*ARIJ**4*EXP(-RIJ/4.43)*1.6022E-12
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)
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 ! R