SUBROUTINE GETGEO(IREAD,LABELS,GEO,LOPT,NA,NB,NC,AMS,NATOMS,INT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE "SIZES" DIMENSION GEO(3,*),NA(*),NB(*),NC(*),AMS(*), LOPT(3,*) 1,LABELS(*) LOGICAL INT ************************************************************************ * * GETGEO READS IN THE GEOMETRY. THE ELEMENT IS SPECIFIED BY IT'S * CHEMICAL SYMBOL, OR, OPTIONALLY, BY IT'S ATOMIC NUMBER. * * ON INPUT IREAD = CHANNEL NUMBER FOR READ, NORMALLY 5 * AMS = DEFAULT ATOMIC MASSES. * * ON OUTPUT LABELS = ATOMIC NUMBERS OF ALL ATOMS, INCLUDING DUMMIES. * GEO = INTERNAL COORDINATES, IN ANGSTROMS, AND DEGREES. * LOPT = INTEGER ARRAY, A '1' MEANS OPTIMIZE THIS PARAMETER, * '0' MEANS DO NOT OPTIMIZE, AND A '-1' LABELS THE * REACTION COORDINATE. * NA = INTEGER ARRAY OF ATOMS (SEE DATA INPUT) * NB = INTEGER ARRAY OF ATOMS (SEE DATA INPUT) * NC = INTEGER ARRAY OF ATOMS (SEE DATA INPUT) * ATMASS = ATOMIC MASSES OF ATOMS. ************************************************************************ COMMON /PATH / DUMMY,REACT(3,66), DUMM1,DUMM2 COMMON /ATMASS/ ATMASS(NUMATM) COMMON /KEYWRD/ KEYWRD DIMENSION ISTART(40), XYZ(3,NUMATM), VALUE(4) LOGICAL LEADSP, IRCDRC CHARACTER*80 KEYWRD CHARACTER ELEMNT(107)*2, LINE*80, SPACE*1, NINE*1,ZERO*1, 1TAB*1, COMMA*1, STRING*80, ELE*2, TURN*1 DATA (ELEMNT(I),I=1,107)/'H','HE', 1 'LI','BE','B','C','N','O','F','NE', 2 'NA','MG','AL','SI','P','S','CL','AR', 3 'K','CA','SC','TI','V','CR','MN','FE','CO','NI','CU', 4 'ZN','GA','GE','AS','SE','BR','KR', 5 'RB','SR','Y','ZR','NB','MO','TC','RU','RH','PD','AG', 6 'CD','IN','SN','SB','TE','I','XE', 7 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD','TB','DY', 8 'HO','ER','TM','YB','LU','HF','TA','W','RE','OS','IR','PT', 9 'AU','HG','TL','PB','BI','PO','AT','RN', 1 'FR','RA','AC','TH','PA','U','NP','PU','AM','CM','BK','CF','XX', 2 'FM','MD','CB','++','+','--','-','TV'/ DATA COMMA,SPACE,NINE,ZERO/',',' ','9','0'/ TAB=CHAR(9) IRCDRC=(INDEX(KEYWRD,'IRC')+INDEX(KEYWRD,'DRC') .NE.0) ILOWA = ICHAR('a') ILOWZ = ICHAR('z') ICAPA = ICHAR('A') NATOMS=0 NUMAT=0 10 READ(IREAD,'(A)',END=90,ERR=180)LINE IF(LINE.EQ.' ') GO TO 90 NATOMS=NATOMS+1 IF(NATOMS.GT.NUMATM)THEN WRITE(6,'(//10X,''**** MAX. NUMBER OF ATOMS ALLOWED:'',I4)') 1NUMATM STOP ENDIF * CLEAN THE INPUT DATA ************************************************************************ DO 20 I=1,80 ILINE=ICHAR(LINE(I:I)) IF(ILINE.GE.ILOWA.AND.ILINE.LE.ILOWZ) THEN LINE(I:I)=CHAR(ILINE+ICAPA-ILOWA) ENDIF 20 CONTINUE ************************************************************************ ICOMMA=ICHAR(COMMA) ITAB=ICHAR(TAB) DO 30 I=1,80 KHAR=ICHAR(LINE(I:I)) IF(KHAR.EQ.ICOMMA.OR.KHAR.EQ.ITAB)LINE(I:I)=SPACE 30 CONTINUE * * INITIALIZE ISTART TO INTERPRET BLANKS AS ZERO'S DO 40 I=1,10 40 ISTART(I)=80 * FIND INITIAL DIGIT OF ALL NUMBERS, CHECK FOR LEADING SPACES FOLLOWED * BY A CHARACTER AND STORE IN ISTART LEADSP=.TRUE. NVALUE=0 DO 50 I=1,80 IF (LEADSP.AND.LINE(I:I).NE.SPACE) THEN NVALUE=NVALUE+1 ISTART(NVALUE)=I END IF LEADSP=(LINE(I:I).EQ.SPACE) 50 CONTINUE * * ESTABLISH THE ELEMENT'S NAME AND ISOTOPE, CHECK FOR ERRORS OR E.O.DATA * WEIGHT=0.D0 STRING=LINE(ISTART(1):ISTART(2)-1) IF( STRING(1:1) .GE. ZERO .AND. STRING(1:1) .LE. NINE) THEN * ATOMIC NUMBER USED: NO ISOTOPE ALLOWED LABEL=READA(STRING,1) IF (LABEL.EQ.0) GO TO 80 IF (LABEL.LT.0.OR.LABEL.GT.107) THEN WRITE(6,'('' ILLEGAL ATOMIC NUMBER'')') GO TO 190 END IF GO TO 70 END IF * ATOMIC SYMBOL USED REAL=ABS(READA(STRING,1)) IF (REAL.LT.1.D-15) THEN * NO ISOTOPE ELE=STRING(1:2) ELSE WEIGHT=REAL IF( STRING(2:2) .GE. ZERO .AND. STRING(2:2) .LE. NINE) THEN ELE=STRING(1:1) ELSE ELE=STRING(1:2) END IF END IF * CHECK FOR ERROR IN ATOMIC SYMBOL IF(ELE(1:1).EQ.'-'.AND.ELE(2:2).NE.'-')ELE(2:2)=' ' DO 60 I=1,107 IF(ELE.EQ.ELEMNT(I)) THEN LABEL=I GO TO 70 END IF 60 CONTINUE WRITE(6,'('' UNRECOGNIZED ELEMENT NAME: ('',A,'')'')')ELE GOTO 190 * * ALL O.K. * 70 IF (LABEL.NE.99) NUMAT=NUMAT+1 IF(WEIGHT.NE.0.D0)THEN WRITE(6,'('' FOR ATOM'',I4,'' ISOTOPIC MASS:'' 1 ,F15.5)')NATOMS, WEIGHT ATMASS(NUMAT)=WEIGHT ELSE IF(LABEL .NE. 99) ATMASS(NUMAT)=AMS(LABEL) ENDIF LABELS(NATOMS) =LABEL GEO(1,NATOMS) =READA(LINE,ISTART(2)) GEO(2,NATOMS) =READA(LINE,ISTART(4)) GEO(3,NATOMS) =READA(LINE,ISTART(6)) IF(IRCDRC)THEN TURN=LINE(ISTART(3):ISTART(3)) IF(TURN.EQ.'T')THEN LOPT(1,NATOMS)=1 IF(NATOMS.EQ.1)WRITE(6,'(A)')' IN DRC MONITOR POTENTIAL ENER 1GY'//' TURNING POINTS' ELSE LOPT(1,NATOMS)=0 ENDIF TURN=LINE(ISTART(5):ISTART(5)) IF(TURN.EQ.'T')THEN LOPT(2,NATOMS)=1 ELSE LOPT(2,NATOMS)=0 ENDIF TURN=LINE(ISTART(7):ISTART(7)) IF(TURN.EQ.'T')THEN LOPT(3,NATOMS)=1 ELSE LOPT(3,NATOMS)=0 ENDIF ELSE LOPT(1,NATOMS) =READA(LINE,ISTART(3)) LOPT(2,NATOMS) =READA(LINE,ISTART(5)) LOPT(3,NATOMS) =READA(LINE,ISTART(7)) ENDIF NA(NATOMS) =READA(LINE,ISTART(8)) NB(NATOMS) =READA(LINE,ISTART(9)) NC(NATOMS) =READA(LINE,ISTART(10)) GOTO 10 * * ALL DATA READ IN, CLEAN UP AND RETURN * 80 NATOMS=NATOMS-1 90 NA(2)=1 IF(NATOMS.GT.3)THEN INT=(NA(4).NE.0) ELSE IF(GEO(2,3).LT.10.AND.NATOMS.EQ.3) 1WRITE(6,'(//10X,'' WARNING: INTERNAL COORDINATES ARE ASSUMED -'',/ 210X,'' FOR THREE-ATOM SYSTEMS '',//)') INT=.TRUE. ENDIF C C READ IN VELOCITY VECTOR, IF PRESENT C IF(INDEX(KEYWRD,'VELO').GT.0)THEN IF(INT)THEN WRITE(6,'(A)')' COORDINATES MUST BE CARTESIAN WHEN VELOCITY 1VECTOR IS USED.' STOP ENDIF C# WRITE(6,'(/10X,A)')'INITIAL VELOCITY VECTOR FOR DRC' DO 110 I=1,NATOMS READ(5,'(A)') LINE CALL NUCHAR(LINE,VALUE,NDMY) IF(NDMY.NE.3)THEN WRITE(6,'(/10X,A)') 1' THERE MUST BE EXACTLY THREE VELOCITY DATA PER LINE' STOP ENDIF DO 100 J=1,3 100 REACT(J,I+2)=VALUE(J) C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(VALUE(J),J=1,3) 110 CONTINUE DO 120 I=1,3 DO 120 J=1,2 120 REACT(I,J)=GEO(I,J+1)-GEO(I,1) C C NOW TO ROTATE VELOCITY VECTOR TO SUIT INTERNAL COORDINATE DEFINITION C C C ROTATE AROUND THE 1-2 X-AXIS TO AS TO ELIMINATE REACT(3,2) C (PUT ATOM 2 IN X-Y PLANE) SA=REACT(3,1)/SQRT(REACT(2,1)**2+REACT(3,1)**2+1.D-20) CA=SIGN(SQRT(1.D0-SA**2),REACT(2,1)) C# LABELS(NATOMS+1)=1 C# LABELS(NATOMS+2)=1 C# WRITE(6,*)' FIRST ROTATION, ABOUT 1-2 X-AXIS' DO 130 I=1,NATOMS+2 TEMP1= REACT(2,I)*CA+REACT(3,I)*SA TEMP2=-REACT(2,I)*SA+REACT(3,I)*CA REACT(2,I)=TEMP1 REACT(3,I)=TEMP2 C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(REACT(J,I),J=1,3) 130 CONTINUE C ROTATE AROUND THE 1-2 Z-AXIS TO AS TO ELIMINATE REACT(2,2) C (PUT ATOM 2 ON X AXIS) CA=REACT(1,1)/SQRT(REACT(2,1)**2+REACT(1,1)**2+1.D-20) SA=SIGN(SQRT(1.D0-CA**2),REACT(2,1)) C# WRITE(6,*)' SECOND ROTATION, ABOUT 1-2 Z-AXIS' DO 140 I=1,NATOMS+2 TEMP1= REACT(1,I)*CA+REACT(2,I)*SA TEMP2=-REACT(1,I)*SA+REACT(2,I)*CA REACT(1,I)=TEMP1 REACT(2,I)=TEMP2 C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(REACT(J,I),J=1,3) 140 CONTINUE C ROTATE AROUND THE 2-3 X-AXIS TO AS TO ELIMINATE REACT(3,3) C (PUT ATOM 3 ON X-Y PLANE) SA=REACT(3,2)/SQRT(REACT(2,2)**2+REACT(3,2)**2+1.D-20) CA=SIGN(SQRT(1.D0-SA**2),REACT(2,2)) C# WRITE(6,*)' THIRD ROTATION, ABOUT 2-3 X-AXIS' DO 150 I=1,NATOMS+2 TEMP1= REACT(2,I)*CA+REACT(3,I)*SA TEMP2=-REACT(2,I)*SA+REACT(3,I)*CA REACT(2,I)=TEMP1 REACT(3,I)=TEMP2 C# WRITE(6,'(2X,A2,2X,3F13.5)')ELEMNT(LABELS(I)),(REACT(J,I),J=1,3) 150 CONTINUE C C STRIP OFF FIRST TWO COORDINATES; THESE WERE THE COORDINATE AXIS C DEFINITIONS C DO 160 I=1,NATOMS DO 160 J=1,3 160 REACT(J,I)=REACT(J,I+2) ENDIF IF( .NOT. INT ) THEN DO 170 I=1,NATOMS DO 170 J=1,3 170 XYZ(J,I)=GEO(J,I) DEGREE=90.D0/ASIN(1.D0) CALL XYZINT(XYZ,NATOMS,NA,NB,NC,DEGREE,GEO) IF(ABS(GEO(2,3)-180.D0).LT.1.D-4.OR.ABS(GEO(2,3)).LT.1.D-4) 1THEN WRITE(6,'(A)')' DUE TO PROGRAM BUG, THE FIRST THREE ATOMS MU 1ST NOT LIE IN A STRAIGHT LINE.' STOP ENDIF ELSEIF (.NOT.IRCDRC) THEN IF(LOPT(1,1)+LOPT(2,1)+LOPT(3,1)+LOPT(2,2)+LOPT(3,2)+ 1 LOPT(3,3) .GT. 0)THEN LOPT(1,1)=0 LOPT(2,1)=0 LOPT(3,1)=0 LOPT(2,2)=0 LOPT(3,2)=0 LOPT(3,3)=0 WRITE(6,'(//10X,'' AN UNOPTIMIZABLE GEOMETRIC PARAMETER HAS 1'',/10X,'' BEEN MARKED FOR OPTIMIZATION. THIS IS A NON-FATAL '' 2,''ERROR'')') ENDIF ENDIF IF(NA(3).EQ.0) THEN NB(3)=1 NA(3)=2 ENDIF RETURN * ERROR CONDITIONS 180 IF(IREAD.EQ.5) THEN WRITE(6,'( '' ERROR DURING READ AT ATOM NUMBER '', I3 )')NATOMS ELSE NATOMS=0 RETURN ENDIF 190 J=NATOMS-1 WRITE(6,'('' DATA CURRENTLY READ IN ARE '')') DO 200 K=1,J 200 WRITE(6,210)LABELS(K),(GEO(JJ,K),LOPT(JJ,K),JJ=1,3), 1NA(K),NB(K),NC(K) 210 FORMAT(I4,2X,3(F10.5,2X,I2,2X),3(I2,1X)) STOP END