PROGRAM MXDINP C*************************************************************** C*** PROGRAM : MXDINPUT *** C*** [ CREATE FILE07, FILE05.DAT, AND FILE10 ] *** C*** VERSION 840923-2 BY KATS *** C*** PC-VERSION 880115-0 BY KATS *** C*** REVICED FOR JCPE 900414 BY KATS *** C*** Reviced for low symmetry structure *** C*** Prepare for large unit cell 910314 by Kats *** C*** Integrated version (MD and XD) 910718 by Kats *** C*** 911023 by Kats *** C*** Multi-componets 920324 by Kats *** C*** New file05.dat format 930104 by Kats *** C*** Diatomic molecule 971014 by Kats *** C*** Water film 980319 by Kats *** C*** sheet 991023 by Kats *** C*** Formats and LEL=10 010620 by Kats *** C*** Formats 021130 by Kats *** C*** Rod and drop 040816 by Kats *** C*************************************************************** C PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) C C LNI : Maximum number of atoms in a MD basic cell C LEL : Maximum number of atom species C LNA : Maximum number of atoms in a crystal asymmetric unit C LAT : Maximum number of atoms in a cristal unit cell C LSY : Maximum number of symmetry operations C COMMON /ATOMS/ MA,NA,P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT),NION(LEL) REAL *8 P COMMON /MDDATA/NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q COMMON /SYMMT/ NS,NL, RS(3,3,LSY),TS(3,LSY), TL(3,4) COMMON /TRAJEC/NPT,NPTP,P0C(3,LAT),JON(LAT),XYZ(3,LNA),XYZH(3,LNA) COMMON /MUDANA/ IXD(LNI), JXD(LNI) COMMON /RANDM3/ RD(3), IR,JR,KR, LR,MR,NR INTEGER *4 IR,JR,KR, LR,MR,NR COMMON /ANAME/ ATOM(LNA),ATM(LAT),TITLE(15),ADX(LEL),AOX(LEL), * FLNAME(19) CHARACTER *4 ATOM, ATM, TITLE, ADX, AOX CHARACTER *24 FLNAME C COMMON /ATOMSC/ ION(3,LAM) COMMON /ATOMSD/ WGT(LAM), CHG(LAM), AOI(LAM), BOI(LAM), COI(LAM), * RAD(LAM) CHARACTER *4 ION C INTEGER *4 IONI(LEL) CHARACTER *4 HEX, ANS, RANS, DUM, MXD CHARACTER *12 XNAME(300),YNAME, STRUCT C FLNAME( 1) = 'MXDINPUT.F ' FLNAME( 2) = '2004-08-16-00 ' C flname( 3) = 'WinPC ' C FLNAME( 3) = 'UNIX ' C FLNAME( 5) = 'file05.dat ' FLNAME( 6) = 'file06.dat ' FLNAME( 7) = 'file07.dat ' FLNAME( 8) = 'file08.dat ' FLNAME( 9) = 'file09p.dat ' FLNAME(11) = 'file09v.dat ' FLNAME(10) = 'file10.dat ' FLNAME(12) = ' ' C FLNAME(18) = 'c:\mxd\xtaldata.dat ' FLNAME(19) = '/usr/MXD/xtaldata.dat ' C C CALL ATMDAT CALL RNDMIZ C IF (FLNAME(3).EQ.'UNIX ') THEN FLNAME(18) = FLNAME(19) END IF C C OPEN ( 4, FILE='CON', STATUS='NEW') OPEN (16, FILE=FLNAME(6), STATUS='UNKNOWN', * ACCESS='SEQUENTIAL',FORM='FORMATTED') C 11 OPEN ( 8,FILE=FLNAME(18), STATUS='OLD', * ACCESS='SEQUENTIAL',FORM='FORMATTED') C C ------------------------------------------- Search for a structure MXD = 'XD' 111 WRITE (6,2022) 2022 FORMAT('Structure type (A12)? (see XTALDATA.DAT, or ', * 'type "LIST")' / * '............',12X,'"CHAOS" for liquid or glass') c READ (5,1002) STRUCT 1002 FORMAT (A12) IF (STRUCT.EQ.' ') STOP 888 IF (STRUCT.EQ.'CANCEL'.OR.STRUCT.EQ.'cancel') STOP 999 C C -------------------------------------------------- LIQUID OR GLASS IF (STRUCT.EQ.'CHAOS '.OR.STRUCT.EQ.'chaos ' .OR. * STRUCT.EQ.'LIQUID'.OR.STRUCT.EQ.'liquid' .OR. * STRUCT.EQ.'GLASS '.OR.STRUCT.EQ.'glass ') THEN MXD = 'MD' CALL LIQUID write (6,*) 'Random structure ',adx GO TO 777 END IF C C ----------------------------------------------- List of structures IF (STRUCT.EQ.'LIST '.OR.STRUCT.EQ.'list ') THEN XNAME(1) = 'CHAOS ' NX = 1 400 READ (8,1004,END=490) YNAME 1004 FORMAT (A12) IF (YNAME.NE.'::::::::::::') THEN IF (YNAME.NE.'............') GO TO 400 NX = NX + 1 READ (8,1004) XNAME(NX) GO TO 400 END IF 490 WRITE (6,4888) (XNAME(N),N=1,NX) 4888 FORMAT (1X,6A13) REWIND 8 GO TO 111 END IF C C ------------------------------------------------- Search structure REWIND 8 222 READ (8,1002,END=999) YNAME IF (YNAME.EQ.'::::::::::::') GO TO 111 IF (YNAME.NE.STRUCT) GO TO 222 IF (YNAME.NE.STRUCT) GO TO 222 IF (YNAME.NE.STRUCT) GO TO 222 C C ------------------------------------------------ Xtal initial data C Input from database, and display READ (8,1001) TITLE 1001 FORMAT (18A4) WRITE (6,2003) TITLE 2003 FORMAT (1X,18A4) READ (8,1010) (BOX(I),I=1,6) 1010 FORMAT (6F10.5) WRITE (6,2010) BOX 2010 FORMAT (5X,' A=',F7.4,' B=',F7.4,' C=',F7.4,' (',3F9.4,')') IF (ABS(BOX(4)).GT.1.01) BOX(4) = COS(BOX(4)*3.1415926/180.) IF (ABS(BOX(5)).GT.1.01) BOX(5) = COS(BOX(5)*3.1415926/180.) IF (ABS(BOX(6)).GT.1.01) BOX(6) = COS(BOX(6)*3.1415926/180.) C CALL SPACEG (MA) CALL INATOM (HEX) C READ (8,1020) AOX 1020 FORMAT (10(A4,1X)) write (6,*) AOX CLOSE (8) C WRITE (6,2031) AOX 2031 FORMAT (10(A4,1X),': Atoms in the structure'/ * 'O Si Al Mg Ca Na Cs X Y Z : ', * '(Exampl) Type in please') C ----------------------------- Key in atom names in proper sequence READ (5,1020) ADX C 41 WRITE (6,4402) BOX 4402 FORMAT ('A=',F7.3,' B=',F7.3,' C=',F7.3,'(',3F7.3,')'/ * 'How many cells?' / ' A B C') C ------------------------------ Key in numbers of stacking cells READ (5,4103) NX,NY,NZ 4103 FORMAT (3I5, A4) WRITE (6,*) 'A=',BOX(1)*NX, 'B=',BOX(2)*NY, 'C=',BOX(3)*NZ, * ' Ok?' READ (5,9110) ANS IF (ANS.EQ.'n'.OR.ANS.EQ.'N') GO TO 41 C CALL GENPOS CALL STRCHK C C -------------------------------------------------------- MD and XD 777 WRITE (6,3752) READ (5,4750) TEMP 3752 FORMAT ( 'I',8('-'),'I Temperature(K) ?') 4750 FORMAT (F10.4) IF (TEMP.LT.0.001) TEMP = 300.0 DO 70 I = 1, NTION DO 60 J = 1, 3 CALL RANDOM Q(J,I)= (RD(J)-0.5)*0.05 + 5.0 60 CONTINUE 70 CONTINUE C C ---------------------------------------------- Write on file07.dat DELTAT = -1.0 DTIME = 2.0E-15 NCOMPO = 0 DO 703 IO = 1, LEL ioni(io) = 0 IF (NION(IO).GT.0) THEN NCOMPO = IO DO 701 J = 1, LAM IF (ADX(IO).EQ.ION(1,J)) IONI(IO) = J IF (ADX(IO).EQ.ION(2,J)) IONI(IO) = J IF (ADX(IO).EQ.ION(3,J)) IONI(IO) = J 701 CONTINUE if (ioni(io).gt.0 .and. ioni(io).le.lam) * ADX(IO)= ion(1,ioni(io)) IF (IONI(IO).LE.0) IONI(IO) = LAM END IF 703 CONTINUE C DENSTY = 0.0 DO 708 I = 1, LEL if (nion(i).gt.0) then DENSTY = DENSTY + NION(I) * WGT(IONI(I)) end if 708 CONTINUE DENSTY = DENSTY / (BOX(1)*BOX(2)*BOX(3)*(1.0E-24*6.02214E23)) PRESS = 0.0001 C OPEN ( 7, FILE=FLNAME(7), STATUS='UNKNOWN', * ACCESS='SEQUENTIAL',FORM='FORMATTED') REWIND 7 WRITE (7,7007) TITLE, 0, 0, * NTION, NCOMPO, 0,0,0,0,0,0,0,0,0 WRITE (7,7017) (ADX(I), I=1,LEL) WRITE (7,7018) (NION(I), I=1,LEL) WRITE (7,7018) (IONS(1,I),I=1,LEL) WRITE (7,7018) (IONS(2,I),I=1,LEL) WRITE (7,7070) TEMP,DELTAT,TEMP, PRESS,PRESS,PRESS, * DTIME, (BOX(I),I=1,6), * DENSTY, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 do 717 io=1, ncompo do 717 i=ions(1,io), ions(2,io) WRITE (7,7770) (P(J,I),J=1,3),(Q(J,I),J=1,3), * (P(J,I),J=1,3), io 717 continue 7007 FORMAT (15A4,2I5 / I7,I3, 10I10) 7017 FORMAT ( 10(2X,A4) ) 7018 FORMAT ( 10I6 ) 7070 FORMAT (F10.2, F10.4, F10.2, 3F10.5 / * E10.2, 10X, 6F10.6 / * F10.6, 10X, 6F10.6 ) 7770 FORMAT (3F10.8,1X,3F9.7,1X,3F10.7,1x,i2) ENDFILE (7) REWIND 7 CLOSE (7) C IF (MXD.EQ.'XD ') THEN C --------------------------------------- Write on file10.dat OPEN (10, FILE=FLNAME(10), STATUS='UNKNOWN', * ACCESS='SEQUENTIAL',FORM='FORMATTED') REWIND 10 WRITE (10,8010) (BOX(I),I=1,6), * NX,NY,NZ,NPT,NPTP,NS,HEX,MA WRITE (10,8012) (ATM(J),J=1,MA) WRITE (10,8014) (NID(J),J=1,MA) WRITE (10,8020) (JON(N),(P0C(J,N),J=1,3),IXD(JON(N)), * N=1,NPTP) WRITE (10,8030) (((RS(J,I,N),J=1,3),I=1,3),N=1,NS) WRITE (10,8040) (NSYM(N),N=1,NTION) 8010 FORMAT (3F10.6,3F10.8 / 6I5,5X,A4,I6 ) 8012 FORMAT ( 18(A4) ) 8014 FORMAT ( 18I4 ) 8020 FORMAT (I5,3F10.7,I5) 8030 FORMAT (9F6.1) 8040 FORMAT (12I6) ENDFILE (10) REWIND 10 CLOSE (10) END IF C C ---------------------------------------------------- Confirmations WRITE (6,1201) (ADX(I), NION(I), I=1,LEL) 1201 FORMAT (2X,'*** NO. of ions : ',5(A2,':',I5,3X) / * 2X,' ',5(A2,':',I5,3X) ) IF (STRUCT.NE.'CHAOS '.AND.STRUCT.NE.'chaos ') THEN WRITE (6,5555) NX,NY,NZ, NPT,NPTP 5555 FORMAT (2X,'*** CELL IS ',3I4,6X,'NPT(NPTP) :',I4, * '(',I3,')') END IF WRITE (6,3333) (BOX(I),I=1,6) 3333 FORMAT (2X,'*** A=', F8.4,' B=', F8.4,' C=', F8.4, * ' CA=',F7.4,' CB=',F7.4,' CC=',F7.4) WRITE (6,4801) DENSTY 4801 FORMAT (2X,'*** Density is ',F8.4,' g/cm3') WRITE (6,5801) C C ------------------------------------------------------- File05.dat 5801 FORMAT (' Do you want a new file05.dat (y/n/r) ?') READ (5,9110) ANS 9110 FORMAT (A1) IF (ANS.EQ.'R'.OR. ANS.EQ.'r') GO TO 11 IF (ANS.NE.'Y'.AND.ANS.NE.'y') GO TO 909 C ANREC5 = 50. IF (MXD.EQ.'XD') ANREC5 = 5. DTIME = 2.0 DO 808 IO = 1, LEL IF (ADX(IO).EQ.'H ') DTIME = 0.4 808 CONTINUE C ---------------------------------------------- Write on file05.dat OPEN (15, FILE=FLNAME(5),STATUS='UNKNOWN', * ACCESS='SEQUENTIAL',FORM='FORMATTED') WRITE (15,5000) MXD 5000 FORMAT (A2,'.......I....:....I....:....I....:....I....:....I', * '....:....I....:....I:') WRITE (15,5001) (TITLE(I),I=1,15), ANREC5, DTIME, TEMP 5001 FORMAT ( * 'START ',15A4,':' / * 'ECONOMY 01000. 1000. 100. ',F8.0,2X, * ' 5. ',10X,':' / * 'NOACCUM ',F7.2,3X,' 1.0 0.0 ',30X,':'/ * 'T SCALING ',F9.2,1X,' -0.1 1. ',30X,':' / * 'P NO-CNTL 0.0001 0.0001 0.0001 ' ,30X,':' / * 'V CONST. ',60X,':' / * 'BUSING 4. 0.0 ',20X,':') DO 803 IO = 1, LEL IF (NION(IO).GT.0) THEN ACHG = 0.0 AWGT = 0.0 DO 801 J = 1, LAM IF (ADX(IO).EQ.ION(1,J).OR. * ADX(IO).EQ.ION(2,J).OR. * ADX(IO).EQ.ION(3,J) ) THEN ACHG = CHG(J) AWGT = WGT(J) AI = AOI(J) BI = BOI(J) CI = COI(J) END IF 801 CONTINUE WRITE (15,5007) IO,ADX(IO),REAL(NION(IO)),ACHG,AWGT, * AI, BI, CI 5007 FORMAT(I1,1X,A2,F6.0, 2X,F6.3,2X, 2X,F6.2,2X, F8.3,2X, * F8.3,2X, F8.3,2X, ' :') END IF 803 CONTINUE WRITE (15,5008) 5008 FORMAT (70X,':' / 70X,':') WRITE (15,5000) MXD WRITE (15,5009) 5009 FORMAT ('STOP ',60X,':' /) ENDFILE (15) REWIND 15 CLOSE (15) C 909 ENDFILE (16) CLOSE (16) 999 STOP END C C C ============ C============================================================= Atom data SUBROUTINE ATMDAT PARAMETER (LAM=42) COMMON /ATOMSC/ ION(3,LAM) COMMON /ATOMSD/ WGT(LAM), CHG(LAM), AOI(LAM), BOI(LAM), COI(LAM), * RAD(LAM) CHARACTER *4 ION C CHARACTER *4 AION(LAM), BION(LAM), CION(LAM) REAL *4 AWGT(LAM), ACHG(LAM), AAOI(LAM), ABOI(LAM), * ACOI(LAM), ARAD(LAM) C DATA AION / 'H ', * 'Li ', 'Be ', 'B ', 'C ', 'N ', 'O ', 'F ', * 'Na ', 'Mg ', 'Al ', 'Si ', 'P ', 'S ', 'Cl ', * 'K ', 'Ca ', 'Sc ', 'Ti ', 'Ga ', 'Ge ', 'As ', * 'Se ', 'Br ', * 'Rb ', 'Sr ', 'In ', 'Sn ', 'Sb ', 'Te ', 'I ', * 'Cs ', 'Ba ', 'Pb ', * 'Y ', 'Zr ', * 'He ', 'Ne ', 'Ar ', 'Kr ', 'Xe ', ' ' / C DATA BION / 'h ', * 'li ', 'be ', 'b ', 'c ', 'n ', 'o ', 'f ', * 'na ', 'mg ', 'al ', 'si ', 'p ', 's ', 'cl ', * 'k ', 'ca ', 'sc ', 'ti ', 'ga ', 'ge ', 'as ', * 'se ', 'br ', * 'rb ', 'sr ', 'in ', 'sn ', 'sb ', 'te ', 'i ', * 'cs ', 'ba ', 'pb ', * 'y ', 'zr ', * 'he ', 'ne ', 'ar ', 'kr ', 'xe ', ' ' / C DATA CION / 'H ', * 'LI ', 'BE ', 'B ', 'C ', 'N ', 'O ', 'F ', * 'NA ', 'MG ', 'AL ', 'SI ', 'P ', 'S ', 'CL ', * 'K ', 'CA ', 'SC ', 'TI ', * 'GA ', 'GE ', 'AS ', 'SE ', 'BR ', * 'RB ', 'SR ', 'IN ', 'SN ', 'SB ', 'TE ', 'I ', * 'CS ', 'BA ', 'PB ', * 'Y ', 'ZR ', * 'HE ', 'NE ', 'AR ', 'KR ', 'XE ', ' ' / C DATA AWGT / 1.008, * 6.941, 9.012, 10.81, 12.01, 14.01, 16.00, 19.00, * 22.99, 24.31, 26.98, 28.09, 30.97, 32.07, 35.45, * 39.10, 40.08, 44.96, 47.88, * 69.72, 72.61, 74.92, 78.96, 79.90, * 85.47, 87.62, 114.8, 118.7, 121.8, 127.6, 126.9, * 132.9, 137.3, 207.2, * 88.91, 91.22, * 4.002, 20.18, 39.95, 83.80, 131.29, 0.0 / C DATA ACHG / 1.00, * 1.00, 2.00, 3.00, 4.00, -3.00, -2.00, -1.00, * 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, -1.00, * 1.00, 2.00, 3.00, 4.00, * 3.00, 4.00, 5.00, 6.00, -1.00, * 1.00, 2.00, 3.00, 4.00, 5.00, 6.00, -1.00, * 1.00, 2.00, 4.00, * 3.00, 4.00, * 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 / C DATA AAOI / 0.0, * 0.0, 0.0, 0.720, 0.0, 1.713, 1.629, 1.565, * 1.260, 1.161, 1.064, 1.012, 0.0, 0.0, 1.950, * 1.595, 1.440, 0.0, 1.235, * 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 1.632, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 1.820, 0.0, * 0.0, 0.0, * 1.156, 1.397, 1.860, 2.025, 2.246, 0.0 / C DATA ABOI / 0.0, * 0.0, 0.0, 0.080, 0.0, 0.080, 0.085, 0.085, * 0.080, 0.080, 0.080, 0.080, 0.0, 0.0, 0.090, * 0.080, 0.080, 0.0, 0.080, * 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.080, 0.0, 0.0, 0.0, 0.0, 0.0, * 0.0, 0.080, 0.0, * 0.0, 0.0, * 0.115, 0.120, 0.145, 0.165, 0.180, 0.0 / C DATA ACOI / 0.0, * 0.0, 0.0, 0.0, 0.0, 0.0, 20.00, 20.0, * 10.0, 2.0, 0.0, 0.0, 0.0, 0.0, 30.0, * 15.0, 10.0, 0.0, 0.0, * 0.0, 0.0, 0.0, 0.0, 40.0, * 20.0, 15.0, 0.0, 0.0, 0.0, 0.0, 50.0, * 25.0, 20.0, 0.0, * 0.0, 0.0, * 4.76, 11.04, 38.54, 55.35, 85.57, 0.0 / C DATA ARAD / 0.10, * 0.60, 0.20, 0.10, 0.10, 1.40, 1.35, 1.30, * 1.00, 0.65, 0.25, 0.15, 0.15, 0.10, 1.40, * 1.05, 0.80, 0.50, 0.50, * 0.40, 0.20, 0.20, 0.15, 1.50, * 1.20, 0.95, 0.60, 0.60, 0.70, 0.75, 1.70, * 1.40, 1.10, 1.30, * 0.65, 0.60, * 1.20, 1.30, 1.80, 1.90, 2.00, 0.0 / C DO 100 I = 1, LAM ION(1,I) = AION(I) ION(2,I) = BION(I) ION(3,I) = CION(I) WGT(I) = AWGT(I) CHG(I) = ACHG(I) AOI(I) = AAOI(I) BOI(I) = ABOI(I) COI(I) = ACOI(I) RAD(I) = ARAD(I) 100 CONTINUE RETURN END C C C ========= C================================================================ RANDOM SUBROUTINE RANDOM COMMON /RANDM3/ R1, R2, R3, IR,JR,KR, LR,MR,NR INTEGER *4 IR,JR,KR, LR,MR,NR C INTEGER *4 II,JJ,KK, LL,MM,NN C II = ((JR/3) * (KR/3) + LR) / 2 JJ = ((IR/3) * (KR/3) + MR) / 2 KK = ((IR/3) * (JR/3) + NR) / 2 IR = MOD(II,100000) JR = MOD(JJ,100000) KR = MOD(KK,100000) R1 = FLOAT(IR) * 0.00001 R2 = FLOAT(JR) * 0.00001 R3 = FLOAT(KR) * 0.00001 LL = ((MR/3) * (NR/3) + IR) / 2 MM = ((LR/3) * (NR/3) + JR) / 2 NN = ((LR/3) * (MR/3) + KR) / 2 LR = MOD(LL,100000) MR = MOD(MM,100000) NR = MOD(NN,100000) RETURN C ENTRY RNDMIZ IR = 32723 JR = 23557 KR = 47979 LR = 54893 MR = 16617 NR = 79423 RETURN END C C C ======== C================================================================= CHAOS SUBROUTINE LIQUID PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /ATOMS/ MA,NA,P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT), * NION(LEL) REAL *8 P COMMON /MDDATA/ NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q COMMON /ANAME/ ATOM(LNA),ATM(LAT),TITLE(15),ADX(LEL),AOX(LEL), * FLNAME(19) CHARACTER *4 ATOM, ATM, TITLE, ADX, AOX CHARACTER *24 FLNAME C COMMON /ATOMSC/ ION(3,LAM) COMMON /ATOMSD/ WGT(LAM), CHG(LAM), AOI(LAM), BOI(LAM), COI(LAM), * RAD(LAM) CHARACTER *4 ION C COMMON /RANDM3/ RD(3), IR,JR,KR, LR,MR,NR INTEGER *4 IR,JR,KR, LR,MR,NR C INTEGER *4 IONI(LEL) REAL *4 R(20), rpx(2,3) CHARACTER *1 ANSWER C WRITE (6,3402) READ (5, *) KNK 3402 FORMAT ('KEY IN AN INTEGER FOR RANDOMIZE') DO 10 I = 1, KNK CALL RANDOM 10 CONTINUE C WRITE (6,3432) READ (5, 4430) TITLE 3432 FORMAT ( 'I',58('-'),'I TITLE ?') 4430 FORMAT (15A4) C WRITE (6,3442) (mod(I,10),I=1,LEL) READ (5, 4440) ADX 3442 FORMAT ( 10('ION',I1,' '),' : ATOM SPECIES ?') 4440 FORMAT (10(A4,1X)) DO 30 I = 1, LEL NION(I) = 0 IONI(I) = 0 DO 20 J = 1, LAM IF (ADX(I).EQ.ION(1,J)) IONI(I) = J IF (ADX(I).EQ.ION(2,J)) IONI(I) = J IF (ADX(I).EQ.ION(3,J)) IONI(I) = J 20 CONTINUE if (ioni(i).gt.0 .and. ioni(i).le.lam) adx(i)=ion(1,ioni(i)) IF (IONI(I).LE.0) IONI(I) = LAM 30 CONTINUE C WRITE (6,3462) (mod(I,10),I=1,LEL) READ (5, 4460) NION 3462 FORMAT ( 10('....',I1),' NO. OF IONS ?') 4460 FORMAT (10I5) JO = 0 DO 50 IO = 1, LEL IONS(1,IO) = JO + 1 IONS(2,IO) = IONS(1,IO) + NION(IO) -1 JO = IONS(2,IO) write (6,*) io,adx(io), nion(io),ions(1,io),ions(2,io) 50 continue C 77 WRITE (6,3452) READ (5,4450)(BOX(I),I=1,3) 3452 FORMAT ( 3('I',9('-')),'I BOX EDGE LENGTHS (1-3) ?') 4450 FORMAT (3F10.4) IF (BOX(2).LE.0.1) BOX(2) = BOX(1) IF (BOX(3).LE.0.1) BOX(3) = BOX(2) BOX(4) = 0.0 BOX(5) = 0.0 BOX(6) = 0.0 DENSTY = 0.0 DO 708 I = 1, LEL DENSTY = DENSTY + NION(I) * WGT(IONI(I)) 708 CONTINUE DENSTY = DENSTY / (BOX(1)*BOX(2)*BOX(3)*(1.0E-24*6.02214E23)) WRITE (6,7001) DENSTY 7001 FORMAT (' -----> The density is ', F7.4, 'g/cm3. OK ? (Y/N)') READ (5,7002) ANSWER 7002 FORMAT (A1) IF (ANSWER.NE.'Y' .AND. ANSWER.NE.'y') GO TO 77 C DT = 1.0E-15 C write (6,7201) 7201 format ('Molecule: 1:monoatomic, 2:diatomic, or 3:H2O/CO2 ', * '(1/2/3) ?') read (5,7203) imol 7203 format (i1) C do 5 i=1, 3 rpx(1,i)=0.0 rpx(2,i)=1.0 5 continue write (6,7301) 7301 format ('Structure: 1:bulk, 2:sheet, 3:rod, 4:cube, ', * '5:sphere ?') read (5,7303) istr 7303 format (i1) if (istr.eq.2) then write (6,7311) 7311 format (1x,'Sheet perpendicular to 1:x, 2:y, 3:z ?') read (5,7315) ipx 7315 format (i1) rpx(1,ipx) = 0.20 rpx(2,ipx) = 0.80 write (*,*) 'Range of sheet ',rpx(1,ipx),rpx(2,ipx), * ' OK (y/n)?' READ (5,7002) ANSWER IF (ANSWER.EQ.'N' .OR. ANSWER.EQ.'n') then write (6,*) 'Input range (2F, 0-1)' read (5,*) rpx(1,ipx), rpx(2,ipx) end if end if if (istr.eq.3) then write (6,7312) 7312 format (1x,'Rod along with 1:x, 2:y, 3:z ?') read (5,7315) ipx if (ipx.eq.1) then ipy=2 ipz=3 else if (ipx.eq.2) then ipy=1 ipz=3 else ipy=1 ipz=2 end if rpx(1,ipy) = 0.20 rpx(2,ipy) = 0.80 rpx(1,ipz) = 0.20 rpx(2,ipz) = 0.80 write (*,*) 'Range of axis ',ipy,rpx(1,ipy),rpx(2,ipy) write (*,*) ' axis ',ipz,rpx(1,ipz),rpx(2,ipz), * ' OK (y/n)?' READ (5,7002) ANSWER IF (ANSWER.EQ.'N' .OR. ANSWER.EQ.'n') then write (6,*) 'Input range (2F, 0-1)' read (5,*) rpx(1,ipy), rpx(2,ipy) rpx(1,ipz) = rpx(1,ipy) rpx(2,ipz) = rpx(2,ipy) end if end if if (istr.eq.4) then rpx(1,1) = 0.20 rpx(2,1) = 0.80 rpx(1,2) = 0.20 rpx(2,2) = 0.80 rpx(1,3) = 0.20 rpx(2,3) = 0.80 write (*,*) 'Range of axis x',rpx(1,1),rpx(2,1) write (*,*) 'Range of axis y',rpx(1,2),rpx(2,2) write (*,*) ' axis z',rpx(1,3),rpx(2,3), * ' OK (y/n)?' READ (5,7002) ANSWER IF (ANSWER.EQ.'N' .OR. ANSWER.EQ.'n') then write (6,*) 'Input range (2F, 0-1)' read (5,*) rpx(1,1), rpx(2,1) rpx(1,2) = rpx(1,1) rpx(2,2) = rpx(2,1) rpx(1,3) = rpx(1,1) rpx(2,3) = rpx(2,1) end if end if if (istr.eq.5) then write (*,*) 'Type radius of shpere in Angstrome ' read (5,*) rpx(1,1) end if c write (6,*) adx if (imol.le.1) CALL CHAOS (IONI, IER, istr, ipx, rpx) if (imol.eq.2) call chaos2 (ioni, ier, istr, ipx, rpx) if (imol.eq.3) call chaos3 (ioni, ier, istr, ipx, rpx) IF (IER.NE.0) then write (6,*) 'The volume is too small. Change cell lengths' GO TO 77 end if C DO 100 I = 1, NTION CALL RANDOM DO 100 J = 1, 3 P(J,I) = Q(J,I) / BOX(J) 100 CONTINUE C RETURN END C C C ======== C================================================================= CHAOS SUBROUTINE CHAOS (IONI, IER, istr, ipx, rpx) PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /ATOMS/ MA,NA,P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT), * NION(LEL) REAL *8 P COMMON /MDDATA/ NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q C COMMON /ATOMSC/ ION(3,LAM) COMMON /ATOMSD/ WGT(LAM), CHG(LAM), AOI(LAM), BOI(LAM), COI(LAM), * RAD(LAM) CHARACTER *4 ION C COMMON /RANDM3/ RD(3), IR,JR,KR, LR,MR,NR INTEGER *4 IR,JR,KR, LR,MR,NR C INTEGER *4 IONI(LEL) real *4 rpx(2,3) REAL *8 PX, PY, PZ C IER = 1 NTION = 0 DO 10 I = 1, LEL NTION = NTION + NION(I) 10 CONTINUE ALMIN = AMIN1(BOX(1),BOX(2),BOX(3)) RCUT = ALMIN / 2.0 DO 90 IO = 1, LEL DO 80 I = IONS(1,IO), IONS(2,IO) NTRIAL = 0 20 NTRIAL = NTRIAL + 1 IF (NTRIAL.GT.999) RETURN CALL RANDOM PX = RD(1) * BOX(1) PY = RD(2) * BOX(2) PZ = RD(3) * BOX(3) if (istr.eq.5) then r2=(px-box(1)/2)**2+(py-box(2)/2)**2+(pz-box(3)/2)**2 if (r2.le.rpx(1,1)**2) go to 25 go to 20 end if if (RD(1).lt.rpx(1,1).or.rd(1).gt.rpx(2,1)) go to 20 if (RD(2).lt.rpx(1,2).or.rd(2).gt.rpx(2,2)) go to 20 if (RD(3).lt.rpx(1,3).or.rd(3).gt.rpx(2,3)) go to 20 25 IF (I.LE.1) GO TO 70 RAS = RAD(IONI(IO)) ZAS = CHG(IONI(IO)) DO 40 JO = 1, IO RLIM = (RAS + RAD(IONI(JO))) * 0.7 IF (ZAS*CHG(IONI(JO)).GT.0.0) RLIM = 1.9 RLIM2 = RLIM * RLIM DO 30 J = IONS(1,JO), IONS(2,JO) IF (J.GE.I) GO TO 70 DX = ABS(PX - Q(1,J)) IF (DX.GT.RCUT) DX = BOX(1) - DX DY = ABS(PY - Q(2,J)) IF (DY.GT.RCUT) DY = BOX(2) - DY DZ = ABS(PZ - Q(3,J)) IF (DZ.GT.RCUT) DZ = BOX(3) - DZ IF (DX*DX+DY*DY+DZ*DZ.LT.RLIM2) GO TO 20 30 CONTINUE 40 CONTINUE C 70 Q(1,I) = PX Q(2,I) = PY Q(3,I) = PZ C WRITE (4,*) I,PX,PY,PZ IF (MOD(I, 100).EQ.0) WRITE (6,1001) I if (mod(i,1000).eq.0) write (6,1003) 80 CONTINUE 90 CONTINUE C WRITE (6,1002) IER = 0 RETURN 1001 FORMAT (1X,I5,$) 1002 FORMAT (/) 1003 format (1x) END C C C ========= C================================================================ CHAOS2 SUBROUTINE CHAOS2 (IONI, IER, istr, ipx, rpx) C Diatomic molecules PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /ATOMS/ MA,NA,P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT), * NION(LEL) REAL *8 P COMMON /MDDATA/ NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q C COMMON /ATOMSC/ ION(3,LAM) COMMON /ATOMSD/ WGT(LAM), CHG(LAM), AOI(LAM), BOI(LAM), COI(LAM), * RAD(LAM) CHARACTER *4 ION C COMMON /RANDM3/ RD(3), IR,JR,KR, LR,MR,NR INTEGER *4 IR,JR,KR, LR,MR,NR C INTEGER *4 IONI(LEL) real *4 rpx(2,3) REAL *8 PX, PY, PZ C write (6,2001) 2001 format (1x,'Input intramolecular atomic distance ', * '(O:1.191, N:1.095)') read (5,2002) Dintra 2002 format (f10.5) c IER = 1 NTION = 0 DO 10 I = 1, LEL NTION = NTION + NION(I) 10 CONTINUE ALMIN = AMIN1(BOX(1),BOX(2),BOX(3)) RCUT = ALMIN / 2.0 DO 90 IO = 1, LEL if (mod(nion(io),2).eq.1) then write (*,*) 'Odd number!',io,nion(io) stop end if DO 80 I = IONS(1,IO), IONS(2,IO), 2 NTRIAL = 0 20 NTRIAL = NTRIAL + 1 IF (NTRIAL.GT.999) RETURN CALL RANDOM if (RD(1).lt.rpx(1,1).or.rd(1).gt.rpx(2,1)) go to 20 if (RD(2).lt.rpx(1,2).or.rd(2).gt.rpx(2,2)) go to 20 if (RD(3).lt.rpx(1,3).or.rd(3).gt.rpx(2,3)) go to 20 PX1 = RD(1) * BOX(1) ! First atom PY1 = RD(2) * BOX(2) PZ1 = RD(3) * BOX(3) CALL RANDOM ! second atom srd = sqrt((rd(1)*box(1))**2 + (rd(2)*box(2))**2 + * (rd(3)*box(3))**2) PX2 = px1 + RD(1) * box(1) * dintra / srd PY2 = py1 + RD(2) * box(2) * dintra / srd PZ2 = pz1 + RD(3) * box(3) * dintra / srd c IF (I.LE.1) GO TO 70 RAS = RAD(IONI(IO)) ZAS = CHG(IONI(IO)) DO 40 JO = 1, IO RLIM = (RAS + RAD(IONI(JO))) * 0.7 IF (ZAS*CHG(IONI(JO)).GT.0.0) RLIM = 1.9 RLIM2 = RLIM * RLIM DO 30 J = IONS(1,JO), IONS(2,JO) IF (J.GE.I) GO TO 25 DX = ABS(PX1 - Q(1,J)) IF (DX.GT.RCUT) DX = BOX(1) - DX DY = ABS(PY1 - Q(2,J)) IF (DY.GT.RCUT) DY = BOX(2) - DY DZ = ABS(PZ1 - Q(3,J)) IF (DZ.GT.RCUT) DZ = BOX(3) - DZ IF (DX*DX+DY*DY+DZ*DZ.LT.RLIM2) GO TO 20 30 CONTINUE 40 CONTINUE c 25 RAS = RAD(IONI(IO)) ZAS = CHG(IONI(IO)) DO 45 JO = 1, IO RLIM = (RAS + RAD(IONI(JO))) * 0.7 IF (ZAS*CHG(IONI(JO)).GT.0.0) RLIM = 1.9 RLIM2 = RLIM * RLIM DO 35 J = IONS(1,JO), IONS(2,JO) IF (J.GE.I) GO TO 70 DX = ABS(PX2 - Q(1,J)) IF (DX.GT.RCUT) DX = BOX(1) - DX DY = ABS(PY2 - Q(2,J)) IF (DY.GT.RCUT) DY = BOX(2) - DY DZ = ABS(PZ2 - Q(3,J)) IF (DZ.GT.RCUT) DZ = BOX(3) - DZ IF (DX*DX+DY*DY+DZ*DZ.LT.RLIM2) GO TO 20 35 CONTINUE 45 CONTINUE C 70 Q(1,I) = PX1 Q(2,I) = PY1 Q(3,I) = PZ1 Q(1,I+1) = PX2 Q(2,I+1) = PY2 Q(3,I+1) = PZ2 C WRITE (4,*) I,PX,PY,PZ IF (MOD(I, 100).EQ.0) WRITE (6,1001) I if (mod(i,1000).eq.0) write (6,1003) 80 CONTINUE 90 CONTINUE C WRITE (6,1002) IER = 0 RETURN 1001 FORMAT (1X,I5,$) 1002 FORMAT (/) 1003 format (1x) END C C C ========= C================================================================ CHAOS3 SUBROUTINE CHAOS3 (IONI, IER, istr, ipx, rpx) C H2O or CO2 PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /ATOMS/ MA,NA,P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT), * NION(LEL) REAL *8 P COMMON /MDDATA/ NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q C COMMON /ATOMSC/ ION(3,LAM) COMMON /ATOMSD/ WGT(LAM), CHG(LAM), AOI(LAM), BOI(LAM), COI(LAM), * RAD(LAM) CHARACTER *4 ION C COMMON /ANAME/ ATOM(LNA),ATM(LAT),TITLE(15),ADX(LEL),AOX(LEL), * FLNAME(19) CHARACTER *4 ATOM, ATM, TITLE, ADX, AOX CHARACTER *24 FLNAME C COMMON /RANDM3/ RD(3), IR,JR,KR, LR,MR,NR INTEGER *4 IR,JR,KR, LR,MR,NR C INTEGER *4 IONI(LEL) REAL *8 PX, PY, PZ real *4 rpx(2,3) C c write (6,2001) c2001 format (1x,'Input intramolecular atomic distance ', c * '(O-H:0.95, C-O:1.13 A)'/) c read (5,2002) Dintra c2002 format (f10.5) c write (6,*) adx if (ADX(1).EQ.'O ') then io1 = 1 io2 = 2 rlim2 = 2.9**2 dintra = 0.95 end if if (ADX(2).eq.'C ') then io1 = 2 io2 = 1 rlim2 = 3.7**2 dintra = 1.13 end if write (6,*) 'rlim2=',rlim2,' dintra=',dintra c IER = 1 NTION = 0 DO 10 I = 1, LEL NTION = NTION + NION(I) 10 CONTINUE ALMIN = AMIN1(BOX(1),BOX(2),BOX(3)) RCUT = ALMIN / 2.0 C IO = io1 ! O of H2O or C of CO2 nio2 = ions(1,io2)-1 DO 80 I = IONS(1,IO), IONS(2,IO) NTRIAL = 0 20 NTRIAL = NTRIAL + 1 IF (NTRIAL.GT.2999) RETURN CALL RANDOM PX1 = RD(1) * BOX(1) ! First atom PY1 = RD(2) * BOX(2) PZ1 = RD(3) * BOX(3) if (istr.eq.5) then r2=(px1-box(1)/2)**2+(py1-box(2)/2)**2+(pz1-box(3)/2)**2 if (r2.le.rpx(1,1)**2) go to 25 go to 20 end if if (RD(1).lt.rpx(1,1).or.rd(1).gt.rpx(2,1)) go to 20 if (RD(2).lt.rpx(1,2).or.rd(2).gt.rpx(2,2)) go to 20 if (RD(3).lt.rpx(1,3).or.rd(3).gt.rpx(2,3)) go to 20 c c if (rd(2).lt.0.07 .or. rd(2).gt.0.93) go to 20 C 25 IF (I.LE.1) GO TO 40 DO 30 J = IONS(1,IO), IONS(1,IO)+I-1 DX = ABS(PX1 - Q(1,J)) IF (DX.GT.RCUT) DX = BOX(1) - DX DY = ABS(PY1 - Q(2,J)) IF (DY.GT.RCUT) DY = BOX(2) - DY DZ = ABS(PZ1 - Q(3,J)) IF (DZ.GT.RCUT) DZ = BOX(3) - DZ IF (DX*DX+DY*DY+DZ*DZ.LT.RLIM2) GO TO 20 30 CONTINUE C 40 CALL RANDOM ! second atom srd = sqrt(((2*rd(1)-1.0)*box(1))**2 + * ((2*rd(2)-1.0)*box(2))**2 + * ((2*rd(3)-1.0)*box(3))**2) PX2 = px1 + (2*RD(1)-1.0) * box(1) * dintra / srd PY2 = py1 + (2*RD(2)-1.0) * box(2) * dintra / srd PZ2 = pz1 + (2*RD(3)-1.0) * box(3) * dintra / srd 50 call random srd = sqrt(((2*rd(1)-1.0)*box(1))**2 + * ((2*rd(2)-1.0)*box(2))**2 + * ((2*rd(3)-1.0)*box(3))**2) PX3 = px1 + (2*RD(1)-1.0) * box(1) * dintra / srd !3rd atom PY3 = py1 + (2*RD(2)-1.0) * box(2) * dintra / srd PZ3 = pz1 + (2*RD(3)-1.0) * box(3) * dintra / srd vv = ( (px2-px1)*(px3-px1) + * (py2-py1)*(py3-py1) + * (pz2-pz1)*(pz3-pz1) ) / dintra**2 c write (6,*) vv if (io1.eq.1) then if (vv.ge.0.0 .or. vv.le.-0.5) go to 50 end if if (io1.eq.2 .and. vv.gt.-0.5) go to 50 c Q(1,I) = PX1 Q(2,I) = PY1 Q(3,I) = PZ1 nio2 = nio2 +1 Q(1,nio2) = PX2 Q(2,nio2) = PY2 Q(3,nio2) = PZ2 nio2 = nio2 + 1 Q(1,nio2) = PX3 Q(2,nio2) = PY3 Q(3,nio2) = PZ3 c WRITE (6,*) I,PX1,PY1,PZ1 IF (MOD(I, 100).EQ.0) WRITE (6,1001) I if (mod(i,1000).eq.0) write (6,1003) 80 CONTINUE C WRITE (6,1002) IER = 0 RETURN 1001 FORMAT (1X,I5,$) 1002 FORMAT (/) 1003 format (1X) END C C C ========= C================================================================ SPACEG SUBROUTINE SPACEG (MA) C Read reduced symmetry operations C Prepare full set of symmetry operations PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /SYMMT/ NS,NL, RS(3,3,LSY),TS(3,LSY), TL(3,4) DIMENSION TAL(3,5),ISQ(4),Q(6,7),MS(7),SG(7),PG(7),TR(2,3) CHARACTER *1 LTP CHARACTER *4 HMS, PG, SG DATA TAL/ .0,.0,.0, .0,.5,.5, .5,.0,.5, .5,.5,.0, .5,.5,.5 / DATA ISQ/ 1, 99, 3, 4/ DATA SG / '1', '2', '21', 'M', 'C', 'N', ' ' / DATA PG / '-1', '2/M', '21/C', '21/M', '21/C', '21/N', '2/C'/ DATA MS / 1, 2, 2, 2, 2, 2, 2 / DATA Q / 1., 1., 1., 0.0,0. ,0. , -1., 1.,-1., 0.0,0. ,0.0, * -1., 1.,-1., 0.0, .5,0. , 1.,-1., 1., 0.0,0. ,0.0, * 1.,-1., 1., 0.0, .5,0.5, 1.,-1., 1., 0.5, .5, .5, * -1., 1.,-1., 0.0,0.0,0.5/ DO 20 N = 1, LSY DO 20 I = 1, 3 TS(I,N) = 0.0 DO 10 J = 1, 3 10 RS(I,J,N) = 0.0 20 RS(I,I,N) = 1.0 READ (8,1010) LTP, HMS, NS, IC, MA1,MA2 MA = MA1 + MA2 C WRITE (6,2010) LTP, HMS C NL = 2 IF (LTP.EQ.'P') ISQ(2) = 0 IF (LTP.EQ.'A') ISQ(2) = 2 IF (LTP.EQ.'B') ISQ(2) = 3 IF (LTP.EQ.'C') ISQ(2) = 4 IF (LTP.EQ.'I') ISQ(2) = 5 IF (LTP.EQ.'R') THEN NL = 3 ISQ(2) = 2 TAL(1,2) = 1.0 / 3.0 TAL(2,2) = 2.0 / 3.0 TAL(3,2) = 2.0 / 3.0 TAL(1,3) = 2.0 / 3.0 TAL(2,3) = 1.0 / 3.0 TAL(3,3) = 1.0 / 3.0 END IF IF (ISQ(2).EQ.0) NL = 1 IF (ISQ(2).GE.0.AND.ISQ(2).LE.5) GO TO 100 IF (LTP.NE.'F ') STOP 111 ISQ(2) = 2 NL = 4 100 DO 110 J = 1, NL DO 110 I = 1, 3 ISQJ = ISQ(J) TL(I,J) = TAL(I,ISQJ) 110 CONTINUE DO 200 IG = 1, 7 JC = 0 IF (NS.EQ.0.AND.HMS.EQ.SG(IG)) GO TO 220 JC = 1 IF (NS.EQ.0.AND.HMS.EQ.PG(IG)) GO TO 220 200 CONTINUE NS = NS + 1 DO 30 I = 2, NS READ (8,3030) TR, ((RS(K,J,I),K=1,3),J=1,3) DO 30 J = 1,3 TS(J,I) = 0.0 IF (TR(2,J).GT.0.0) TS(J,I) = TR(1,J) / TR(2,J) 30 CONTINUE IF (IC.EQ.0) JC = 0 GO TO 400 C 220 NS = MS(IG) DO 300 I = 1, 3 RS(I,I,2) = Q(I,IG) TS(I,2) = Q(I+3,IG) 300 CONTINUE C 400 IF (JC.GE.1) THEN DO 320 N = 1, NS NSN = NS + N DO 320 J = 1, 3 TS(J,NSN) = -TS(J,N) DO 320 K = 1, 3 RS(K,J,NSN) = -RS(K,J,N) 320 CONTINUE NS = NS * 2 END IF WRITE (6,2020) LTP,HMS, NS, NL C WRITE (6,2030) (((RS(K,J,I),K=1,3),J=1,3),(TS(J,I),J=1,3),I=1,NS) C WRITE (6,2033) ((TL(J,I),J=1,3),I=1,NL) RETURN 1010 FORMAT (A1, A4, 5X, 4I5) C2010 FORMAT (/6X, A1, A4 /) 2020 FORMAT (5X,'Space group :',A1,A4 / * 5x,'No. of symmetry operations is',I3 / * 5x,'No. of lattice points is',I3) 2030 FORMAT (11X, 3F4.1,1X,3F4.1,1X,3F4.1, 2X, 3F6.3) 2033 FORMAT (11X, 3F4.1,1X,3F4.1,1X,3F4.1, 1X, 3F4.1) 3030 FORMAT (6F1.0,9F2.0) END C C C ========= C================================================================ INATOM SUBROUTINE INATOM (HEX) C Read atom data in an asymmetric unit C Expand atoms over a crystal unit cell PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /MDDATA/ NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q COMMON /ANAME/ ATOM(LNA),ATM(LAT),TITLE(15),ADX(LEL),AOX(LEL), * FLNAME(19) CHARACTER *4 ATOM, ATM, TITLE, ADX, AOX CHARACTER *24 FLNAME COMMON /SYMMT/ NS,NL, RS(3,3,LSY),TS(3,LSY), TL(3,4) COMMON /TRAJEC/ NPT,NPTP,P0C(3,LAT),JON(LAT),XYZ(3,LNA), * XYZH(3,LNA) COMMON /ATOMS/ MA,NA, P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT), * NION(LEL) REAL *8 P CHARACTER *4 HEX CHARACTER *1 ANS REAL *8 X(3) C ---------------------------- Read atom data in an asymmetric unit DO 110 M = 1, MA READ ( 8,1020) ATOM(M),ID(M),VA,(P(I,M),I=1,3) WRITE(16,2020) ATOM(M),ID(M),VA,(P(I,M),I=1,3) IDD(M) = M ISYM(M) = 1 DO 100 I = 1, 3 IF (P(I,M).LT.0.0) P(I,M) = P(I,M) + 1.0 IF (P(I,M).GE.1.0) P(I,M) = P(I,M) - 1.0 XYZ(I,M) = P(I,M) XYZH(I,M) = P(I,M) 100 CONTINUE write (6,2020) atom(m),id(m),va,(xyz(i,m),i=1,3) 110 CONTINUE NA = MA WRITE (6,2010) NA C ------------------------------ Expand atoms over crystal unit cell DO 700 N = 1, MA DO 700 M = 1, NS DO 700 L = 1, NL DO 200 J = 1,3 X(J) = P(1,N)*RS(1,J,M) + P(2,N)*RS(2,J,M) + * P(3,N)*RS(3,J,M) + TS(J,M) + TL(J,L) IF (X(J).LT.0.0) X(J) = X(J) - AINT(X(J)-1.0) IF (X(J).GE.1.0) X(J) = X(J) - AINT(X(J)) 200 CONTINUE DO 300 J = 1, NA IF (ID(N).EQ.ID(J)) THEN D = ABS(X(1)-P(1,J)) + ABS(X(2)-P(2,J)) + * ABS(X(3)-P(3,J)) IF (D.LT.0.001) GO TO 700 END IF 300 CONTINUE NA = NA + 1 ID(NA) = ID(N) IDD(NA) = IDD(N) ISYM(NA) = M + NS * (L - 1) DO 400 J = 1, 3 P(J,NA) = X(J) 400 CONTINUE 700 CONTINUE WRITE (6,2030) NA C C ------------------------------------ Hxagonal or Rhombohedral case HEX = ' ' IF (ABS(BOX(4))+ABS(BOX(5))+ABS(BOX(6)+0.5).GE.1.E-4) RETURN WRITE (6,*) ' The crystal system is hexagonal or trigonal.' WRITE (6,*) ' Transform to orthogonal system ?' READ (5,*) ANS IF (ANS.EQ.'n' .OR. ANS.EQ.'N') RETURN C HEX = 'HEX ' IF (NL.EQ.3) HEX = 'HEXR' BOX(6) = 0.0 BOX(2) = BOX(2) * SIN(2.0943951) * 2.0 N = NA DO 5 I = 1, NA P(1,I) = P(1,I) - P(2,I) * 0.5 P(2,I) = P(2,I) * 0.5 XYZH(1,I) = P(1,I) XYZH(2,I) = P(2,I) IF (P(1,I).LT.0.0) P(1,I) = P(1,I) + 1.0 N = N + 1 P(1,N) = P(1,I) - 0.5 IF (P(1,N).LT.0.0) P(1,N) = P(1,N) + 1.0 P(2,N) = P(2,I) + 0.5 P(3,N) = P(3,I) ID(N) = ID(I) IDD(N) = IDD(I) ISYM(N) = ISYM(I) + NS * NL 5 CONTINUE NA = N WRITE (6,3003) (I,ID(I),(P(J,I),J=1,3),I=1,N) WRITE (6,3333) (BOX(I),I=1,6) RETURN 1020 FORMAT (A4,1X, I5, F10.5,3F10.5, 3I5) 2010 FORMAT (5X,'The number of atoms in an asymmetric unit is ',I3) 2020 FORMAT (11X,A4, 3X, I5, F10.2,3F10.4, 3X, 3I7) 2030 FORMAT (5X,'The number of atoms in a cell is ',I3) 3003 FORMAT (3X,I3,I2,1X,3F6.3, I5,I2,1X,3F6.3, I5,I2,1X,3F6.3) 3333 FORMAT (1X, 'A=',F9.5,' B=',F9.5,' C=',F9.5,' ALPHA=',F6.2, * ' BETA=',F6.2,' GAMMA=',F6.2) END C C C ========= C================================================================ GENPOS SUBROUTINE GENPOS PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /ATOMS/ MA,NA,P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT), * NION(LEL) REAL *8 P COMMON /ANAME/ ATOM(LNA),ATM(LAT),TITLE(15),ADX(LEL),AOX(LEL), * FLNAME(19) CHARACTER *4 ATOM, ATM, TITLE, ADX, AOX CHARACTER *24 FLNAME COMMON /RANDM3/ RD(3), IR,JR,KR, LR,MR,NR INTEGER *4 IR,JR,KR, LR,MR,NR COMMON /MDDATA/ NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q COMMON /TRAJEC/ NPT,NPTP,P0C(3,LAT),JON(LAT),XYZ(3,LNA), * XYZH(3,LNA) COMMON /MUDANA/ IXD(LNI), JXD(LNI) REAL *4 CL(3) INTEGER *4 IH,IM,IS,IMM, IDV(LNI), NCL(3) CHARACTER *4 ATOMI REAL *8 QQ C DO 10 I = 1, LNA NID(I) = 0 ATM(I) = ' ' 10 CONTINUE C NNX = 2 NNY = 2 NNZ = 2 IF (NX.LE.1) NNX = 1 IF (NY.LE.1) NNY = 1 IF (NZ.LE.1) NNZ = 1 C NNS = 0 NTION = 0 JO = 0 DO 50 IO = 1, LEL NION(IO) = 0 IONS(1,IO) = JO + 1 IF (ADX(IO).EQ.' ') GO TO 45 DO 40 I = 1, NA IDI = ID(I) IF (ADX(IO).EQ.AOX(IDI)) THEN NS = 0 DO 20 IX = 0, NX-1 DO 20 IY = 0, NY-1 DO 20 IZ = 0, NZ-1 NS = NS + 1 NTION = NTION + 1 NION(IO) = NION(IO) + 1 IDV(NTION) = IDD(I) NSYM(NTION) = ISYM(I) + (NS-1) * 200 Q(1,NTION) = (P(1,I) + IX) * BOX(1) Q(2,NTION) = (P(2,I) + IY) * BOX(2) Q(3,NTION) = (P(3,I) + IZ) * BOX(3) IF (IX+1.EQ.NNX .AND. IY+1.EQ.NNY .AND. * IZ+1.EQ.NNZ ) NNS = NS 20 CONTINUE END IF 40 CONTINUE 45 IONS(2,IO) = IONS(1,IO) + NION(IO) - 1 JO = IONS(2,IO) 50 CONTINUE WRITE (6,*) 'The total number of atoms in a basic cell is',NTION C NU = 0 DO 530 IO = 1, LEL IF (NION(IO).LE.0) GO TO 530 N1 = IONS(1,IO) N2 = IONS(2,IO) NMAX = IDV(N1) NMIN = IDV(N1) N11 = N2 - 1 DO 520 I = N1, N11 N22 = I + 1 DO 515 J = N22, N2 IF (NMAX.LT.IDV(I)) NMAX = IDV(I) IF (NMIN.GT.IDV(I)) NMIN = IDV(I) IF (NMAX.LT.IDV(J)) NMAX = IDV(J) IF (NMIN.GT.IDV(J)) NMIN = IDV(J) IF (IDV(I).GT.IDV(J)) THEN IDI = IDV(I) IDV(I) = IDV(J) IDV(J) = IDI NSM = NSYM(I) NSYM(I) = NSYM(J) NSYM(J) = NSM DO 510 K = 1,3 QQ = Q(K,I) Q(K,I) = Q(K,J) Q(K,J) = QQ 510 CONTINUE END IF 515 CONTINUE 520 CONTINUE DO 540 I = N1, N2 IDVI = IDV(I) IDVJ = IDVI - (NMIN - NU) + 1 IDV(I) = IDVJ ATM(IDVJ) = ATOM(IDVI) IXD(I) = IDVJ JXD(I) = IDVI 540 CONTINUE NU = NU + (NMAX - NMIN) + 1 530 CONTINUE C CL(1) = BOX(1) CL(2) = BOX(2) CL(3) = BOX(3) NCL(1) = NX NCL(2) = NY NCL(3) = NZ BOX(1) = BOX(1) * NX BOX(2) = BOX(2) * NY BOX(3) = BOX(3) * NZ DO 70 I = 1, NTION Q(1,I) = Q(1,I) / BOX(1) Q(2,I) = Q(2,I) / BOX(2) Q(3,I) = Q(3,I) / BOX(3) 70 CONTINUE C 80 NPT = 0 DO 65 I = 1, NTION II = IDV(I) NID(II) = NID(II) + 1 CALL RANDOM DO 60 J = 1, 3 P(J,I) = Q(J,I) 60 CONTINUE IF (NSYM(I).GE.200) GO TO 65 NPT = NPT + 1 JON(NPT) = I DO 63 J = 1, 3 P0C(J,NPT) = P(J,I) * NCL(J) 63 CONTINUE 65 CONTINUE C WRITE (6,*) NPT,NID(1),NID(2),NID(3),NID(4) NPTP = NPT C IF (NX+NY+NZ.LE.3) RETURN IF (NPT.GT.64) GO TO 800 NNX = 2 NNY = 2 NNZ = 2 IF (NCL(1).LE.1) NNX = 1 IF (NCL(2).LE.1) NNY = 1 IF (NCL(3).LE.1) NNZ = 1 DO 95 I = 1, NTION MMS = NSYM(I) / 200 IF (MMS.LE.0) GO TO 95 IF (MMS.NE.(NNS-1)) GO TO 95 NPT = NPT + 1 JON(NPT) = I P0C(1,NPT) = P(1,I) * NCL(1) P0C(2,NPT) = P(2,I) * NCL(2) P0C(3,NPT) = P(3,I) * NCL(3) 95 CONTINUE C 800 NPTP = NPT DO 880 I = 1, NTION DO 810 J = 1, 3 IF (P(J,I)*NCL(J).GT.1.01) GO TO 880 810 CONTINUE DO 820 J = 1, 3 IF (P(J,I)*NCL(J).GT.0.99999) GO TO 830 820 CONTINUE GO TO 880 830 NPTP = NPTP + 1 JON(NPTP) = I DO 870 J = 1, 3 P0C(J,NPTP) = P(J,I) * NCL(J) 870 CONTINUE 880 CONTINUE C RETURN END C C C ========= C================================================================ STRCHK SUBROUTINE STRCHK PARAMETER (LNI=62783, LEL=10, LNA=682, LAT=LNA*8, LSY=192, * LAM=42 ) COMMON /ATOMS/ MA,NA,P(3,LNI),ID(LAT),IDD(LAT),ISYM(LAT), * NION(LEL) REAL *8 P COMMON /MDDATA/ NTION,BOX(6),NX,NY,NZ,IONS(2,LEL),Q(3,LNI), * NID(LNA),NSYM(LNI) REAL *8 Q COMMON /TRAJEC/ NPT,NPTP,P0C(3,LAT),JON(LAT),XYZ(3,LNA), * XYZH(3,LNA) COMMON /MUDANA/ IXD(LNI), JXD(LNI) COMMON /ANAME/ ATOM(LNA),ATM(LAT),TITLE(15),ADX(LEL),AOX(LEL), * FLNAME(19) CHARACTER *4 ATOM, ATM, TITLE, ADX, AOX CHARACTER *24 FLNAME REAL *4 PXYZ(3),DB(100) CHARACTER *4 AB(100),ABAB C ABCXY = BOX(1)*BOX(2)*box(6) ABCYZ = BOX(2)*BOX(3)*box(4) ABCZX = BOX(3)*BOX(1)*box(5) C DO 700 M = 1, MA PXYZ(1) = XYZH(1,M) / NX PXYZ(2) = XYZH(2,M) / NY PXYZ(3) = XYZH(3,M) / NZ NB = 0 DO 200 I = 1, 100 AB(I) = ' ' DB(I) = 0.0 200 CONTINUE DO 300 N = 1, NTION DX = PXYZ(1) - P(1,N) DY = PXYZ(2) - P(2,N) DZ = PXYZ(3) - P(3,N) IF (DX.GT. 0.5) DX = DX - 1.0 IF (DX.LT.-0.5) DX = DX + 1.0 IF (DY.GT. 0.5) DY = DY - 1.0 IF (DY.LT.-0.5) DY = DY + 1.0 IF (DZ.GT. 0.5) DZ = DZ - 1.0 IF (DZ.LT.-0.5) DZ = DZ + 1.0 DD = (DX*BOX(1))**2 + (DY*BOX(2))**2 + (DZ*BOX(3))**2 + * 2.0*DX*DY*ABCXY + 2.0*DY*DZ*ABCYZ + 2.0*DZ*DX*ABCZX IF (DD.GT.0.1.AND.DD.LE.16.00) THEN IF (NB.GE.99) GO TO 300 NB = NB + 1 DB(NB) = SQRT(DD) AB(NB) = ATOM(JXD(N)) END IF 300 CONTINUE DO 450 I1 = 1, NB-1 DO 400 I2 = I1+1, NB IF (DB(I1).GT.DB(I2)) THEN DDDDDD = DB(I1) DB(I1) = DB(I2) DB(I2) = DDDDDD ABAB = AB(I1) AB(I1) = AB(I2) AB(I2) = ABAB END IF 400 CONTINUE 450 CONTINUE WRITE (16,2223) M, ATOM(M), (DB(II),AB(II),II=1,NB) WRITE ( 6,2222) M, ATOM(M), (DB(II),AB(II),II=1,6) 700 CONTINUE WRITE (6,*) 'Structure check completed !' RETURN 2222 FORMAT (1X,I4, 1x,A4,2X,11(1X,F4.2,'(',A4,')') / * 12X, 11(1X,F4.2,'(',A4,')') / * 12X, 11(1X,F4.2,'(',A4,')') / * 12X, 11(1X,F4.2,'(',A4,')') ) 2223 FORMAT (1X,I4, 1x,A4,2X,11(1X,F5.3,'(',A4,')') / * 12X, 11(1X,F5.3,'(',A4,')') / * 12X, 11(1X,F5.3,'(',A4,')') / * 12X, 11(1X,F5.3,'(',A4,')') ) END