SUBROUTINE POWSQ(XPARAM, NVAR, FUNCT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE 'SIZES' DIMENSION XPARAM(*) COMMON /MESAGE/ IFLEPO,ISCF ********************************************************************** * * POWSQ OPTIMIZES THE GEOMETRY BY MINIMISING THE GRADIENT NORM. * THUS BOTH GROUND AND TRANSITION STATE GEOMETRIES CAN BE * CALCULATED. IT IS ROUGHLY EQUIVALENT TO FLEPO, FLEPO MINIMIZES * THE ENERGY, POWSQ MINIMIZES THE GRADIENT NORM. * * ON ENTRY XPARAM = VALUES OF PARAMETERS TO BE OPTIMIZED. * NVAR = NUMBER OF PARAMETERS TO BE OPTIMIZED. * * ON EXIT XPARAM = OPTIMIZED PARAMETERS. * FUNCT = HEAT OF FORMATION IN KCALS. * ********************************************************************** C ***** ROUTINE PERFORMS A LEAST SQUARES MINIMIZATION ***** C ***** OF A FUNCTION WHICH IS A SUM OF SQUARES. ***** C ***** INITIALLY WRITTEN BY J.W. MCIVER JR. AT SUNY/ ***** C ***** BUFFALO, SUMMER 1971. REWRITTEN AND MODIFIED ***** C ***** BY A.K. AT SUNY BUFFALO AND THE UNIVERSITY OF ***** C ***** TEXAS. DECEMBER 1973 ***** C COMMON /GEOVAR/ NDUM,LOC(2,MAXPAR), IDUMY, XARAM(MAXPAR) COMMON /GEOM / GEO(3,NUMATM) COMMON /LAST / LAST COMMON /KEYWRD/ KEYWRD COMMON /TIME / TIME0 COMMON /NUMSCF/ NSCF COMMON /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR), 1 LOCDEP(MAXPAR) COMMON /GRADNT/ GRAD(MAXPAR),GNFINA COMMON /TIMDMP/ TLEFT, TDUMP COMMON /GEOKST/ NATOMS,LABELS(NUMATM), 1 NA(NUMATM),NB(NUMATM),NC(NUMATM) COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 2 NCLOSE,NOPEN,NDUMY,FRACT COMMON /NUMCAL/ NUMCAL COMMON /SIGMA1/ GNEXT, AMIN, ANEXT COMMON /SIGMA2/ GNEXT1(MAXPAR), GMIN1(MAXPAR) COMMON /NLLCOM/ HESS(MAXPAR,MAXPAR),BMAT(MAXPAR,MAXPAR), 1PMAT(MAXPAR*MAXPAR) COMMON /SCRACH/ PVEC DIMENSION IPOW(9), SIG(MAXPAR), 1 E1(MAXPAR), E2(MAXPAR), 2 P(MAXPAR), WORK(MAXPAR), 3 PVEC(MAXPAR*MAXPAR), EIG(MAXPAR), Q(MAXPAR) LOGICAL DEBUG, RESTRT, TIMES, OKF, ROUGH, SCF1, RESFIL, LOG CHARACTER*241 KEYWRD SAVE ICALCN DATA ICALCN /0/ IF(ICALCN.NE.NUMCAL) THEN ICALCN=NUMCAL RESTRT=(INDEX(KEYWRD,'RESTART') .NE. 0) LOG=(INDEX(KEYWRD,'NOLOG') .EQ. 0) SCF1=(INDEX(KEYWRD,'1SCF') .NE. 0) ROUGH=.FALSE. TIME1=SECOND() TIME2=TIME1 ICYC=0 TIMES=(INDEX(KEYWRD,'TIME') .NE. 0) TLAST=TLEFT RESFIL=.FALSE. STEP=0.02D0 LAST=0 ILOOP=1 NAT3=NUMAT*3 XINC=0.00529167D0 RHO2=1.D-4 TOL2=4.D-1 IF(INDEX(KEYWRD,'PREC') .NE. 0) TOL2=1.D-2 IF(INDEX(KEYWRD,'GNORM') .NE. 0) THEN TOL2=READA(KEYWRD,INDEX(KEYWRD,'GNORM')) IF(TOL2.LT.0.01.AND.INDEX(KEYWRD,' LET').EQ.0)THEN WRITE(6,'(/,A)')' GNORM HAS BEEN SET TOO LOW, RESET TO 0 1.01' TOL2=0.01D0 ENDIF ENDIF DEBUG = (INDEX(KEYWRD,'POWSQ') .NE. 0) IF(RESTRT) THEN C C RESTORE STORED DATA C IPOW(9)=0 CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,ILOOP,BMAT,IPOW) IF(SCF1) GOTO 390 NSCF=IPOW(8) DO 10 I=1,NVAR GRAD(I)=GMIN1(I) 10 GNEXT1(I)=GMIN1(I) WRITE(6,'('' XPARAM'',6F10.6)')(XPARAM(I),I=1,NVAR) IF(ILOOP .GT. 0) THEN C# ILOOP=ILOOP+1 WRITE(6,'(//10X,'' RESTARTING AT POINT'',I3)')ILOOP ELSE WRITE(6,'(//10X,''RESTARTING IN OPTIMISATION'', 1 '' ROUTINES'')') ENDIF ENDIF * * DEFINITIONS: NVAR = NUMBER OF GEOMETRIC VARIABLES = 3*NUMAT-6 * ENDIF NVAR=ABS(NVAR) IF(DEBUG) THEN WRITE(6,'('' XPARAM'')') WRITE(6,'(5(2I3,F10.4))')(LOC(1,I),LOC(2,I),XPARAM(I),I=1,NVAR) ENDIF IF( .NOT. RESTRT) THEN DO 20 I=1,NVAR 20 GRAD(I)=0.D0 CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.) ENDIF IF(DEBUG) THEN WRITE(6,'('' STARTING GRADIENTS'')') WRITE(6,'(3X,8F9.4)')(GRAD(I),I=1,NVAR) ENDIF GMIN=SQRT(DOT(GRAD,GRAD,NVAR)) GLAST=GMIN DO 30 I=1,NVAR GNEXT1(I)=GRAD(I) GMIN1(I)=GNEXT1(I) 30 CONTINUE C C NOW TO CALCULATE THE HESSIAN MATRIX. C IF(ILOOP.LT.0) GOTO 140 C C CHECK THAT HESSIAN HAS NOT ALREADY BEEN CALCULATED. C ILPR=ILOOP DO 50 ILOOP=ILPR,NVAR TIME1=SECOND() XPARAM(ILOOP)=XPARAM(ILOOP) + XINC CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.) IF(SCF1) GOTO 390 IF(DEBUG)WRITE(6,'(I3,12(8F9.4,/3X))') 1 ILOOP,(GRAD(IF),IF=1,NVAR) GRAD(ILOOP)=GRAD(ILOOP)+1.D-5 XPARAM(ILOOP)=XPARAM(ILOOP) - XINC DO 40 J=1,NVAR 40 HESS(ILOOP,J)=-(GRAD(J)-GNEXT1(J))/XINC TIME2=SECOND() TSTEP=TIME2-TIME1 IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)') 1 TSTEP, TLEFT IF(TLAST-TLEFT.GT.TDUMP)THEN TLAST=TLEFT RESFIL=.TRUE. IPOW(9)=2 I=ILOOP IPOW(8)=NSCF CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) ENDIF IF( TLEFT .LT. TSTEP*2.D0) THEN C C STORE RESULTS TO DATE. C IPOW(9)=1 I=ILOOP IPOW(8)=NSCF CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) STOP ENDIF 50 CONTINUE C ***** SCALE -HESSIAN- MATRIX ***** IF( DEBUG) THEN WRITE(6,'(//10X,''UN-NORMALIZED HESSIAN MATRIX'')') DO 60 I=1,NVAR 60 WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR) ENDIF DO 80 I=1,NVAR SUM = 0.0D0 DO 70 J=1,NVAR 70 SUM = SUM+HESS(I,J)**2 80 WORK(I) = 1.0D0/SQRT(SUM) DO 90 I=1,NVAR DO 90 J=1,NVAR 90 HESS(I,J) = HESS(I,J)*WORK(I) IF( DEBUG) THEN WRITE(6,'(//10X,''HESSIAN MATRIX'')') DO 100 I=1,NVAR 100 WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR) ENDIF C ***** INITIALIZE B MATIRX ***** DO 120 I=1,NVAR DO 110 J=1,NVAR 110 BMAT(I,J) = 0.0D0 120 BMAT(I,I) = WORK(I)*2.D0 ************************************************************************ * * THIS IS THE START OF THE BIG LOOP TO OPTIMIZE THE GEOMETRY * ************************************************************************ ILOOP=-99 TSTEP=TSTEP*4 130 CONTINUE IF(TLAST-TLEFT.GT.TDUMP)THEN TLAST=TLEFT RESFIL=.TRUE. IPOW(9)=2 I=ILOOP IPOW(8)=NSCF CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) ENDIF IF( TLEFT .LT. TSTEP*2.D0) THEN C C STORE RESULTS TO DATE. C IPOW(9)=1 I=ILOOP IPOW(8)=NSCF CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) IFLEPO=-1 RETURN ENDIF 140 CONTINUE C ***** FORM-A- DAGGER-A- IN PA SLONG WITH -P- ***** IJ=0 DO 160 J=1,NVAR DO 160 I=1,J IJ=IJ+1 SUM = 0.0D0 DO 150 K=1,NVAR 150 SUM = SUM + HESS(I,K)*HESS(J,K) 160 PMAT(IJ) = SUM DO 180 I=1,NVAR SUM = 0.0D0 DO 170 K=1,NVAR 170 SUM = SUM-HESS(I,K)*GMIN1(K) 180 P(I) = -SUM L=0 IF(DEBUG) THEN WRITE(6,'(/10X,''P MATRIX IN POWSQ'')') CALL VECPRT(PMAT,NVAR) ENDIF CALL RSP(PMAT,NVAR,NVAR,EIG,PVEC) C ***** CHECK FOR ZERO EIGENVALUE ***** C# WRITE(6,'('' EIGS IN POWSQ:'')') C# WRITE(6,'(6F13.8)')(EIG(I),I=1,NVAR) IF(EIG(1).LT.RHO2) GO TO 240 INDC = 2 C ***** IF MATRIX IS NOT SINGULAR FORM INVERSE ***** C ***** BY BACK TRANSFORMING THE EIGENVECTORS ***** IJ=0 DO 200 I=1,NVAR DO 200 J=1,I IJ=IJ+1 SUM = 0.0D0 DO 190 K=1,NVAR 190 SUM = SUM+PVEC((K-1)*NVAR+J)*PVEC((K-1)*NVAR+I)/EIG(K) 200 PMAT(IJ) = SUM C ***** FIND -Q- VECTOR ***** L=0 IL=L+1 L=IL+I-1 DO 230 I=1,NVAR SUM = 0.0D0 DO 210 K=1,I IK=(I*(I-1))/2+K 210 SUM = SUM+PMAT(IK)*P(K) IP1=I+1 DO 220 K=IP1,NVAR IK=(K*(K-1))/2+I 220 SUM=SUM+PMAT(IK)*P(K) 230 Q(I) = SUM GO TO 260 240 CONTINUE C ***** TAKE -Q- VECTOR AS EIGENVECTOR OF ZERO ***** C ***** EIGENVALUE ***** DO 250 I=1,NVAR 250 Q(I) = PVEC(I) 260 CONTINUE C ***** FIND SEARCH DIRECTION ***** DO 270 I=1,NVAR SIG(I) = 0.0D0 DO 270 J=1,NVAR 270 SIG(I) = SIG(I) + Q(J)*BMAT(I,J) C ***** DO A ONE DIMENSIONAL SEARCH ***** IF (DEBUG) THEN WRITE(6,'('' SEARCH VECTOR'')') WRITE(6,'(8F10.5)')(SIG(I),I=1,NVAR) ENDIF CALL SEARCH(XPARAM, ALPHA, SIG, NVAR, GMIN, OKF, FUNCT) IF( NVAR .EQ. 1) GOTO 390 C C FIRST WE ATTEMPT TO OPTIMIZE GEOMETRY USING SEARCH. C IF THIS DOES NOT WORK, THEN SWITCH TO LINMIN, WHICH ALWAYS WORKS, C BUT IS TWICE AS SLOW AS SEARCH. C ROUGH= ( .NOT. OKF) RMX = 0.0D0 DO 280 K=1,NVAR RT = ABS(GMIN1(K)) IF(RT.GT.RMX)RMX = RT 280 CONTINUE IF(RMX.LT.TOL2) GO TO 390 C ***** TWO STEP ESTIMATION OF DERIVATIVES ***** DO 290 K=1,NVAR 290 E1(K) = (GMIN1(K)-GNEXT1(K))/(AMIN-ANEXT) RMU = DOT(E1,GMIN1,NVAR)/DOT(GMIN1,GMIN1,NVAR) DO 300 K=1,NVAR 300 E2(K) = E1(K) - RMU*GMIN1(K) C ***** SCALE -E2- AND -SIG- ***** SK = 1.0D0/SQRT(DOT(E2,E2,NVAR)) DO 310 K=1,NVAR 310 SIG(K) = SK*SIG(K) DO 320 K=1,NVAR 320 E2(K) = SK*E2(K) C ***** FIND INDEX OF REPLACEMENT DIRECTION ***** PMAX = -1.0D+20 DO 330 I=1,NVAR IF(ABS(P(I)*Q(I)).LE.PMAX) GO TO 330 PMAX = ABS(P(I)*Q(I)) ID = I 330 CONTINUE C ***** REPLACE APPROPRIATE DIRECTION AND DERIVATIVE *** DO 340 K=1,NVAR 340 HESS(ID,K) = -E2(K) C ***** REPLACE STARTING POINT ***** DO 350 K=1,NVAR 350 BMAT(K,ID) = SIG(K)/0.529167D0 DO 360 K=1,NVAR 360 GNEXT1(K) = GMIN1(K) GLAST = GMIN INDC = 1 TIME1=TIME2 TIME2=SECOND() TLEFT=TLEFT-TIME2+TIME0 TSTEP=TIME2-TIME1 ICYC=ICYC+1 IF(RESFIL)THEN WRITE(6,370)MIN(TLEFT,9999999.9D0), 1MIN(GMIN,999999.999D0),FUNCT IF(LOG)WRITE(11,370)MIN(TLEFT,9999999.9D0), 1MIN(GMIN,999999.999D0),FUNCT 370 FORMAT(' RESTART FILE WRITTEN, TIME LEFT:',F9.1, 1' GRAD.:',F10.3,' HEAT:',G14.7) RESFIL=.FALSE. ELSE WRITE(6,380)ICYC,MIN(TSTEP,9999.99D0), 1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT IF(LOG)WRITE(11,380)ICYC,MIN(TSTEP,9999.99D0), 1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT 380 FORMAT(' CYCLE:',I5,' TIME:',F6.1,' TIME LEFT:',F9.1, 1' GRAD.:',F10.3,' HEAT:',G14.7) ENDIF IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)') 1TSTEP, TLEFT GO TO 130 390 CONTINUE DO 400 I=1,NVAR 400 GRAD(I)=0.D0 LAST=1 CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.) DO 410 I=1,NVAR 410 GRAD(I)=GMIN1(I) GNFINA=SQRT(DOT(GRAD,GRAD,NVAR)) IFLEPO=11 IF(SCF1)IFLEPO=13 RETURN END