SUBROUTINE GMETRY(GEO,COORD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE "SIZES" COMMON /GEOKST/ NATOMS,LABELS(NUMATM), 1NA(NUMATM),NB(NUMATM),NC(NUMATM) 2 /EULER / TVEC(3,3), ID COMMON /REACTN/ STEP, GEOA(3,NUMATM), GEOVEC(3,NUMATM),COLCST COMMON /GEOOK/ IGEOOK DIMENSION GEO(3,*),COORD(3,*) CHARACTER *15 NDIMEN(4) LOGICAL FIRST, GEOOK DATA FIRST/.TRUE./, NDIMEN/' MOLECULE ',' POLYMER ', 1'LAYER STRUCTURE',' SOLID '/ C*********************************************************************** C C GMETRY COMPUTES COORDINATES FROM BOND-ANGLES AND LENGTHS. C *** IT IS ADAPTED FROM THE PROGRAM WRITTEN BY M.J.S. DEWAR. C C THREE SEPARATE OPTIONS EXIST WITHIN GMETRY. THESE ARE: C (A) IF NA(1) IS EQUAL TO 99 (IMPOSSIBLE UNDER NORMAL CIRCUMSTANCES) C THEN GEO IS ASSUMED TO BE IN CARTESIAN RATHER THAN INTERNAL C COORDINATES, AND COORD IS THEN SET EQUAL TO GEO. C (B) IF STEP IS NON-ZERO (THIS IS THE CASE WHEN "SADDLE" IS USED) C THEN GEO IS FIRST MODIFIED BY SHIFTING THE INTERNAL COORDINATES C ALONG A RADIUS FROM GEOA TO PLACE GEO AT ADISTANCESTEPFROMGEOA. C (C) NORMAL CONVERSION FROM INTERNAL TO CARTESIAN COORDINATESISDONE. C C ON INPUT: C GEO = ARRAY OF INTERNAL COORDINATES. C NATOMS = NUMBER OF ATOMS, INCLUDING DUMMIES. C NA = ARRAY OF ATOM LABELS FOR BOND LENGTHS. C C ON OUTPUT: C COORD = ARRAY OF CARTESIAN COORDINATES C C*********************************************************************** C OPTION (B) GEOOK=(IGEOOK.EQ.99) IF(ABS(STEP) .GT. 1.D-4) THEN SUM=0.D0 DO 10 I=1,NATOMS DO 10 J=1,3 GEOVEC(J,I)=GEO(J,I)-GEOA(J,I) 10 SUM=SUM+GEOVEC(J,I)**2 SUM=SQRT(SUM) ERROR=(SUM-STEP)/SUM ELSE ERROR=0.D0 ENDIF DO 20 I=1,NATOMS DO 20 J=1,3 20 GEO(J,I)=GEO(J,I)-ERROR*GEOVEC(J,I) C OPTION (A) IF(NA(1).EQ.99) THEN DO 30 I=1,3 DO 30 J=1,NATOMS 30 COORD(I,J)=GEO(I,J) GOTO 110 ENDIF C OPTION (C) COORD(1,1)=0.0D00 COORD(2,1)=0.0D00 COORD(3,1)=0.0D00 COORD(1,2)=GEO(1,2) COORD(2,2)=0.0D00 COORD(3,2)=0.0D00 IF(NATOMS.EQ.2) GOTO 110 CCOS=COS(GEO(2,3)) IF(NA(3).EQ.1)THEN COORD(1,3)=COORD(1,1)+GEO(1,3)*CCOS ELSE COORD(1,3)=COORD(1,2)-GEO(1,3)*CCOS ENDIF COORD(2,3)=GEO(1,3)*SIN(GEO(2,3)) COORD(3,3)=0.0D00 DO 100 I=4,NATOMS COSA=COS(GEO(2,I)) MB=NB(I) MC=NA(I) XB=COORD(1,MB)-COORD(1,MC) YB=COORD(2,MB)-COORD(2,MC) ZB=COORD(3,MB)-COORD(3,MC) RBC=1.0D00/SQRT(XB*XB+YB*YB+ZB*ZB) IF (ABS(COSA).LT.1.D0-1.D-12) GO TO 40 C C ATOMS MC, MB, AND (I) ARE COLLINEAR C RBC=GEO(1,I)*RBC*COSA COORD(1,I)=COORD(1,MC)+XB*RBC COORD(2,I)=COORD(2,MC)+YB*RBC COORD(3,I)=COORD(3,MC)+ZB*RBC GO TO 100 C C THE ATOMS ARE NOT COLLINEAR C 40 MA=NC(I) XA=COORD(1,MA)-COORD(1,MC) YA=COORD(2,MA)-COORD(2,MC) ZA=COORD(3,MA)-COORD(3,MC) C C ROTATE ABOUT THE Z-AXIS TO MAKE YB=0, AND XB POSITIVE. IF XYB IS C TOO SMALL, FIRST ROTATE THE Y-AXIS BY 90 DEGREES. C XYB=SQRT(XB*XB+YB*YB) K=-1 IF (XYB.GT.0.1D00) GO TO 50 XPA=ZA ZA=-XA XA=XPA XPB=ZB ZB=-XB XB=XPB XYB=SQRT(XB*XB+YB*YB) K=+1 C C ROTATE ABOUT THE Y-AXIS TO MAKE ZB VANISH C 50 COSTH=XB/XYB SINTH=YB/XYB XPA=XA*COSTH+YA*SINTH YPA=YA*COSTH-XA*SINTH SINPH=ZB*RBC COSPH=SQRT(ABS(1.D00-SINPH*SINPH)) XQA=XPA*COSPH+ZA*SINPH ZQA=ZA*COSPH-XPA*SINPH C C ROTATE ABOUT THE X-AXIS TO MAKE ZA=0, AND YA POSITIVE. C YZA=SQRT(YPA**2+ZQA**2) IF(YZA.LT.1.D-4)GOTO 70 IF(YZA.LT.2.D-2 .AND. .NOT.GEOOK)THEN WRITE(6,'(//20X,'' CALCULATION ABANDONED AT THIS POINT'')') WRITE(6,'(//10X,'' THREE ATOMS BEING USED TO DEFINE THE'',/ 110X,'' COORDINATES OF A FOURTH ATOM, WHOSE BOND-ANGLE IS'')') WRITE(6,'(10X,'' NOT ZERO OR 180 DEGREEES, ARE '', 1''IN AN ALMOST STRAIGHT'')') WRITE(6,'(10X,'' LINE. THERE IS A HIGH PROBABILITY THAT THE 1'',/10X,'' COORDINATES OF THE ATOM WILL BE INCORRECT.'')') WRITE(6,'(//20X,''THE FAULTY ATOM IS ATOM NUMBER'',I4)')I CALL GEOUT WRITE(6,'(//20X,''CARTESIAN COORDINATES UP TO FAULTY ATOM'') 1') WRITE(6,'(//5X,''I'',12X,''X'',12X,''Y'',12X,''Z'')') DO 60 J=1,I 60 WRITE(6,'(I6,F16.5,2F13.5)')J,(COORD(K,J),K=1,3) WRITE(6,'(//6X,'' ATOMS'',I3,'','',I3,'', AND'',I3, 1'' ARE WITHIN'',F7.4,'' ANGSTROMS OF A STRAIGHT LINE'')') 2MC,MB,MA,YZA STOP ENDIF COSKH=YPA/YZA SINKH=ZQA/YZA GOTO 80 70 CONTINUE C C ANGLE TOO SMALL TO BE IMPORTANT C COSKH=1.D0 SINKH=0.D0 80 CONTINUE C C COORDINATES :- A=(XQA,YZA,0), B=(RBC,0,0), C=(0,0,0) C NONE ARE NEGATIVE. C THE COORDINATES OF I ARE EVALUATED IN THE NEW FRAME. C SINA=SIN(GEO(2,I)) SIND=-SIN(GEO(3,I)) COSD=COS(GEO(3,I)) XD=GEO(1,I)*COSA YD=GEO(1,I)*SINA*COSD ZD=GEO(1,I)*SINA*SIND C C TRANSFORM THE COORDINATES BACK TO THE ORIGINAL SYSTEM. C YPD=YD*COSKH-ZD*SINKH ZPD=ZD*COSKH+YD*SINKH XPD=XD*COSPH-ZPD*SINPH ZQD=ZPD*COSPH+XD*SINPH XQD=XPD*COSTH-YPD*SINTH YQD=YPD*COSTH+XPD*SINTH IF (K.LT.1) GO TO 90 XRD=-ZQD ZQD=XQD XQD=XRD 90 COORD(1,I)=XQD+COORD(1,MC) COORD(2,I)=YQD+COORD(2,MC) COORD(3,I)=ZQD+COORD(3,MC) 100 CONTINUE C C *** NOW REMOVE THE TRANSLATION VECTORS, IF ANY, FROM THE ARRAY COOR C 110 CONTINUE K=NATOMS 120 IF(LABELS(K).NE.107) GOTO 130 K=K-1 GOTO 120 130 K=K+1 IF(K.GT.NATOMS) GOTO 180 C C SYSTEM IS A SOLID, OF DIMENSION NATOMS+1-K C L=0 DO 140 I=K,NATOMS L=L+1 MC=NA(I) DO 140 J=1,3 TVEC(J,L)=COORD(J,I)-COORD(J,MC) 140 CONTINUE ID=L IF (FIRST) THEN FIRST=.FALSE. WRITE(6,150)NDIMEN(ID+1) 150 FORMAT(/10X,' THE SYSTEM IS A ',A15,/) IF(ID.EQ.0) GOTO 180 WRITE(6,160) WRITE(6,170)(I,(TVEC(J,I),J=1,3),I=1,ID) 160 FORMAT(/,' UNIT CELL TRANSLATION VECTORS',/ 1/,' X Y Z') 170 FORMAT(' T',I1,' = ',F11.7,' ',F11.7,' ',F11.7) ENDIF 180 CONTINUE J=0 DO 200 I=1,NATOMS IF (LABELS(I).EQ.99.OR.LABELS(I).EQ.107) GO TO 200 J=J+1 DO 190 K=1,3 190 COORD(K,J)=COORD(K,I) 200 CONTINUE CUTOFF=200.D0 RETURN END