SUBROUTINE REACT1(ESCF) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE "SIZES" COMMON /GEOM / GEO(3,NUMATM) DIMENSION GEOA(3,NUMATM), GEOVEC(3,NUMATM), 1 P1STOR(MPACK), P2STOR(MPACK), 2 P3STOR(MPACK), XOLD(MAXPAR), GROLD(MAXPAR) COMMON /GEOKST/ NATOMS,LABELS(NUMATM), 1 NA(NUMATM), NB(NUMATM), NC(NUMATM) COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK) COMMON /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),LOCDEP(MAXPAR) COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR) COMMON /GRADNT/ GRAD(MAXPAR),GNORM COMMON /ISTOPE/ AMS(107) COMMON /GRAVEC/ COSINE COMMON /KEYWRD/ KEYWRD 1 /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 2 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 3 NCLOSE,NOPEN,NDUMY,FRACT COMMON /REACTN/ STEP, GEOA, GEOVEC,CALCST LOGICAL GRADNT, FINISH, XYZ, INT, GOK(2) ************************************************************************ * * REACT1 DETERMINES THE TRANSITION STATE OF A CHEMICAL REACTION. * * REACT WORKS BY USING TWO SYSTEMS SIMULTANEOUSLY, THE HEATS OF * FORMATION OF BOTH ARE CALCULATED, THEN THE MORE STABLE ONE * IS MOVED IN THE DIRECTION OF THE OTHER. AFTER A STEP THE * ENERGIES ARE COMPARED, AND THE NOW LOWER-ENERGY FORM IS MOVED * IN THE DIRECTION OF THE HIGHER-ENERGY FORM. THIS IS REPEATED * UNTIL THE SADDLE POINT IS REACHED. * * IF ONE FORM IS MOVED 3 TIMES IN SUCCESSION, THEN THE HIGHER ENERGY * FORM IS RE-OPTIMIZED WITHOUT SHORTENING THE DISTANCE BETWEEN THE TWO * FORMS. THIS REDUCES THE CHANCE OF BEING CAUGHT ON THE SIDE OF A * TRANSITION STATE. * ************************************************************************ DIMENSION IDUM1(NUMATM), IDUM2(NUMATM), XSTORE(MAXPAR), +IDUM3(NUMATM), COORD(3,NUMATM), IROT(2,3) DIMENSION IDUMMY(3*NUMATM) CHARACTER*80 KEYWRD EQUIVALENCE (IDUMMY,COORD) DATA IROT/1,2,1,3,2,3/ GOLD=0.D0 LINEAR=0 GOK(1)=.FALSE. GOK(2)=.FALSE. XYZ=(INDEX(KEYWRD,' XYZ') .NE. 0) GRADNT=(INDEX(KEYWRD,'GRAD') .NE. 0) I=(INDEX(KEYWRD,' BAR')) STEPMX=0.15D0 IF(I.NE.0) STEPMX=READA(KEYWRD,I) MAXSTP=1000 C C READ IN THE SECOND GEOMETRY. C IF(XYZ) THEN CALL GETGEO(5,LABELS,GEOA,LOC,NA,NB,NC,AMS,NATOMS,INT) ELSE CALL GETGEO(5,IDUM1,GEOA,IDUMMY, 1 IDUM1,IDUM2,IDUM3,AMS,NATOMS,INT) C C IF INTERNAL COORDINATES ARE TO BE USED, CHECK THE CONNECTIVITY C L=0 DO 8 I=1,NATOMS IF(IDUM1(I).NE.NA(I))THEN L=L+1 IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY'')') WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE BOND LABELS ARE D +IFFERENT: ',IDUM1(I),' AND',NA(I) ENDIF IF(IDUM2(I).NE.NB(I))THEN L=L+1 IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY'')') WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE BOND ANGLE LABELS + ARE DIFFERENT:',IDUM2(I),' AND',NB(I) ENDIF IF(IDUM3(I).NE.NC(I))THEN L=L+1 IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY'')') WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE DIHEDRAL LABELS A +RE DIFFERENT: ',IDUM3(I),' AND',NC(I) ENDIF 8 CONTINUE IF(L.NE.0)WRITE(6,'(10X,A)')' CORRECT BEFORE RESUBMISSION' IF(L.NE.0)STOP ENDIF CLOSE (5) TIME0= SECOND() C C SWAP FIRST AND SECOND GEOMETRIES AROUND C SO THAT GEOUT CAN OUTPUT DATA ON SECOND GEOMETRY. C NUMAT2=0 DO 10 I=1,NATOMS IF(LABELS(I).NE.99) NUMAT2=NUMAT2+1 CONST=1.D0 DO 10 J=1,3 X=GEOA(J,I)*CONST CONST=0.0174532925D0 GEOA(J,I)=GEO(J,I) GEO(J,I)=X 10 CONTINUE IF(NUMAT2.NE.NUMAT) THEN WRITE(6,'(//10X,'' NUMBER OF ATOMS IN SECOND SYSTEM IS '', 1''INCORRECT'',/)') WRITE(6,'('' NUMBER OF ATOMS IN FIRST SYSTEM ='',I4)')NUMAT WRITE(6,'('' NUMBER OF ATOMS IN SECOND SYSTEM ='',I4)')NUMAT2 GOTO 270 ENDIF WRITE(6,'(//10X,'' GEOMETRY OF SECOND SYSTEM'',/)') IF(NDEP.NE.0) CALL SYMTRY CALL GEOUT C C CONVERT TO CARTESIAN, IF NECESSARY C IF( XYZ )THEN CALL GMETRY(GEO,COORD) SUMX=0.D0 SUMY=0.D0 SUMZ=0.D0 DO 20 J=1,NUMAT SUMX=SUMX+COORD(1,J) SUMY=SUMY+COORD(2,J) 20 SUMZ=SUMZ+COORD(3,J) SUMX=SUMX/NUMAT SUMY=SUMY/NUMAT SUMZ=SUMZ/NUMAT DO 30 J=1,NUMAT GEO(1,J)=COORD(1,J)-SUMX GEO(2,J)=COORD(2,J)-SUMY 30 GEO(3,J)=COORD(3,J)-SUMZ WRITE(6,'(//,'' CARTESIAN GEOMETRY OF FIRST SYSTEM'',//)') WRITE(6,'(3F14.5)')((GEO(J,I),J=1,3),I=1,NUMAT) SUMX=0.D0 SUMY=0.D0 SUMZ=0.D0 DO 40 J=1,NUMAT SUMX=SUMX+GEOA(1,J) SUMY=SUMY+GEOA(2,J) 40 SUMZ=SUMZ+GEOA(3,J) SUM=0.D0 SUMX=SUMX/NUMAT SUMY=SUMY/NUMAT SUMZ=SUMZ/NUMAT DO 50 J=1,NUMAT GEOA(1,J)=GEOA(1,J)-SUMX GEOA(2,J)=GEOA(2,J)-SUMY GEOA(3,J)=GEOA(3,J)-SUMZ SUM=SUM+(GEO(1,J)-GEOA(1,J))**2 1 +(GEO(2,J)-GEOA(2,J))**2 2 +(GEO(3,J)-GEOA(3,J))**2 50 CONTINUE DO 100 L=3,1,-1 C C DOCKING IS DONE IN STEPS OF 16, 4, AND 1 DEGREES AT A TIME. C CA=COS(4.D0**(L-1)*0.01745329D0) SA=SQRT(ABS(1.D0-CA**2)) DO 90 J=1,3 IR=IROT(1,J) JR=IROT(2,J) DO 80 I=1,10 SUMM=0.D0 DO 60 K=1,NUMAT X = CA*GEOA(IR,K)+SA*GEOA(JR,K) GEOA(JR,K)=-SA*GEOA(IR,K)+CA*GEOA(JR,K) GEOA(IR,K)=X SUMM=SUMM+(GEO(1,K)-GEOA(1,K))**2 1 +(GEO(2,K)-GEOA(2,K))**2 2 +(GEO(3,K)-GEOA(3,K))**2 60 CONTINUE IF(SUMM.GT.SUM) THEN IF(I.GT.1)THEN SA=-SA DO 70 K=1,NUMAT X = CA*GEOA(IR,K)+SA*GEOA(JR,K) GEOA(JR,K)=-SA*GEOA(IR,K)+CA*GEOA(JR,K) GEOA(IR,K)=X 70 CONTINUE GOTO 90 ENDIF SA=-SA ENDIF 80 SUM=SUMM 90 CONTINUE 100 CONTINUE WRITE(6,'(//,'' CARTESIAN GEOMETRY OF SECOND SYSTEM'',//)') WRITE(6,'(3F14.5)')((GEOA(J,I),J=1,3),I=1,NUMAT) WRITE(6,'(//,'' "DISTANCE":'',F13.6)')SUM WRITE(6,'(//,'' REACTION COORDINATE VECTOR'',//)') WRITE(6,'(3F14.5)')((GEOA(J,I)-GEO(J,I),J=1,3),I=1,NUMAT) NA(1)=99 J=0 NVAR=0 DO 120 I=1,NATOMS IF(LABELS(I).NE.99)THEN J=J+1 DO 110 K=1,3 NVAR=NVAR+1 LOC(2,NVAR)=K 110 LOC(1,NVAR)=J LABELS(J)=LABELS(I) ENDIF 120 CONTINUE NATOMS=NUMAT ENDIF C C XPARAM HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEO C XOLD HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEOA C IF(NVAR.EQ.0)THEN WRITE(6,'(///10X,''THERE ARE NO VARIABLES IN THE SADDLE'', 1'' CALCULATION!'')') STOP ENDIF SUM=0.D0 DO 130 I=1,NVAR GROLD(I)=1.D0 XPARAM(I)=GEO(LOC(2,I),LOC(1,I)) XOLD(I)=GEOA(LOC(2,I),LOC(1,I)) 130 SUM=SUM+(XPARAM(I)-XOLD(I))**2 STEP0=SQRT(SUM) ONE=1.D0 DELL=0.1D0 EOLD=-2000.D0 TIME1=SECOND() SWAP=0 DO 230 ILOOP=1,MAXSTP TIME2=SECOND() WRITE(6,'('' TIME='',F9.2)')TIME2-TIME1 TIME1=TIME2 C C THIS METHOD OF CALCULATING 'STEP' IS QUITE ARBITARY, AND NEEDS C TO BE IMPROVED BY INTELLIGENT GUESSWORK! C IF (GNORM.LT.1.D-3)GNORM=1.D-3 WRITE(6,'('' CURRENT BAR, STEPMX, GNORM'',3F12.7)') 1STEP0,STEPMX,GNORM STEP=MIN(SWAP,0.5D0, 6.D0/GNORM, DELL,STEPMX*STEP0+0.005D0) SWAP=SWAP+1.D0 DELL=DELL+0.1 STEP0=STEP0-STEP IF(STEP0.LT.0.01D0) GOTO 240 STEP=STEP0 DO 150 I=1,NVAR 150 XSTORE(I)=XPARAM(I) CALL FLEPO(XPARAM, NVAR, ESCF) IF(LINEAR.EQ.0)THEN LINEAR=(NORBS*(NORBS+1))/2 DO 140 I=1,LINEAR P1STOR(I)=P(I) P2STOR(I)=PA(I) 140 P3STOR(I)=PB(I) ENDIF DO 160 I=1,NVAR 160 XPARAM(I)=GEO(LOC(2,I),LOC(1,I)) WRITE(6,'(//10X,''FOR POINT'',I3)')ILOOP WRITE(6,'('' DISTANCE A - B '',F12.6)')STEP C C NOW TO CALCULATE THE "CORRECT" GRADIENTS, SWITCH OFF 'STEP'. C STEP=0.D0 DO 170 I=1,NVAR 170 GRAD(I)=GROLD(I) CALL COMPFG (XPARAM, .TRUE., FUNCT1,.FALSE.,GRAD,.TRUE.) DO 180 I=1,NVAR 180 GROLD(I)=GRAD(I) IF (GRADNT) THEN WRITE(6,'('' ACTUAL GRADIENTS OF THIS POINT'')') WRITE(6,'(8F10.4)')(GRAD(I),I=1,NVAR) ENDIF WRITE(6,'('' HEAT '',F12.6)')FUNCT1 GNORM=SQRT(DOT(GRAD,GRAD,NVAR)) WRITE(6,'('' GRADIENT NORM '',F12.6)')GNORM COSINE=COSINE*ONE WRITE(6,'('' DIRECTION COSINE'',F12.6)')COSINE CALL GEOUT IF(SWAP.GT.2.9D0 .OR. ILOOP .GT. 3 .AND. COSINE .LT. 0.D0 1 .OR. ESCF .GT. EOLD) 2 THEN IF(SWAP.GT.2.9D0)THEN SWAP=0.D0 ELSE SWAP=0.5D0 ENDIF C C SWAP REACTANT AND PRODUCT AROUND C FINISH=(GOK(1).AND.GOK(2) .AND. COSINE .LT. 0.D0) IF(FINISH) THEN WRITE(6,'(//10X,'' BOTH SYSTEMS ARE ON THE SAME SIDE OF T 1HE '',''TRANSITION STATE -'',/10X,'' GEOMETRIES OF THE SYSTEMS'', 2'' ON EACH SIDE OF THE T.S. ARE AS FOLLOWS'')') DO 190 I=1,NVAR 190 XPARAM(I)=XSTORE(I) CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.) WRITE(6,'(//10X,'' GEOMETRY ON ONE SIDE OF THE TRANSITION 1'','' STATE'')') CALL WRITE(TIME0,FUNCT1) ENDIF WRITE(6,'('' REACTANTS AND PRODUCTS SWAPPED AROUND'')') ONE=-1.D0 EOLD=ESCF SUM=GOLD GOLD=GNORM I=1.7+ONE*0.5 IF(GNORM.GT.10.D0)GOK(I)=.TRUE. GNORM=SUM DO 200 I=1,NUMAT DO 200 J=1,3 X=GEO(J,I) GEO(J,I)=GEOA(J,I) 200 GEOA(J,I)=X DO 210 I=1,NVAR X=XOLD(I) XOLD(I)=XPARAM(I) 210 XPARAM(I)=X C C C SWAP AROUND THE DENSITY MATRICES. C DO 220 I=1,LINEAR X=P1STOR(I) P1STOR(I)=P(I) P(I)=X X=P2STOR(I) P2STOR(I)=PA(I) PA(I)=X X=P3STOR(I) P3STOR(I)=PB(I) PB(I)=X 220 CONTINUE IF(FINISH) GOTO 240 ELSE ONE=1.D0 ENDIF 230 CONTINUE 240 CONTINUE WRITE(6,'('' AT END OF REACTION'')') GOLD=SQRT(DOT(GRAD,GRAD,NVAR)) CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.) GNORM=SQRT(DOT(GRAD,GRAD,NVAR)) DO 250 I=1,NVAR 250 GROLD(I)=XPARAM(I) CALL WRITE(TIME0,FUNCT1) * * THE GEOMETRIES HAVE (A) BEEN OPTIMIZED CORRECTLY, OR * (B) BOTH ENDED UP ON THE SAME SIDE OF THE T.S. * * TRANSITION STATE LIES BETWEEN THE TWO GEOMETRIES * C1=GOLD/(GOLD+GNORM) C2=1.D0-C1 WRITE(6,'('' BEST ESTIMATE GEOMETRY OF THE TRANSITION STATE'')') WRITE(6,'(//10X,'' C1='',F8.3,''C2='',F8.3)')C1,C2 DO 260 I=1,NVAR 260 XPARAM(I)=C1*GROLD(I)+C2*XOLD(I) STEP=0.D0 CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.) CALL WRITE(TIME0,FUNCT1) 270 STOP END