1 SUBROUTINE POWSQ(XPARAM, NVAR, FUNCT) 2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 3 INCLUDE 'SIZES' 4 DIMENSION XPARAM(*) 5 COMMON /MESAGE/ IFLEPO,ISCF 6********************************************************************** 7* 8* POWSQ OPTIMIZES THE GEOMETRY BY MINIMISING THE GRADIENT NORM. 9* THUS BOTH GROUND AND TRANSITION STATE GEOMETRIES CAN BE 10* CALCULATED. IT IS ROUGHLY EQUIVALENT TO FLEPO, FLEPO MINIMIZES 11* THE ENERGY, POWSQ MINIMIZES THE GRADIENT NORM. 12* 13* ON ENTRY XPARAM = VALUES OF PARAMETERS TO BE OPTIMIZED. 14* NVAR = NUMBER OF PARAMETERS TO BE OPTIMIZED. 15* 16* ON EXIT XPARAM = OPTIMIZED PARAMETERS. 17* FUNCT = HEAT OF FORMATION IN KCALS. 18* 19********************************************************************** 20C ***** ROUTINE PERFORMS A LEAST SQUARES MINIMIZATION ***** 21C ***** OF A FUNCTION WHICH IS A SUM OF SQUARES. ***** 22C ***** INITIALLY WRITTEN BY J.W. MCIVER JR. AT SUNY/ ***** 23C ***** BUFFALO, SUMMER 1971. REWRITTEN AND MODIFIED ***** 24C ***** BY A.K. AT SUNY BUFFALO AND THE UNIVERSITY OF ***** 25C ***** TEXAS. DECEMBER 1973 ***** 26C 27 COMMON /GEOVAR/ NDUM,LOC(2,MAXPAR), IDUMY, XARAM(MAXPAR) 28 COMMON /GEOM / GEO(3,NUMATM), XCOORD(3,NUMATM) 29 COMMON /LAST / LAST 30 COMMON /KEYWRD/ KEYWRD 31C ***** Modified by Jiro Toyoda at 1994-05-25 ***** 32C COMMON /TIME / TIME0 33 COMMON /TIMEC / TIME0 34C ***************************** at 1994-05-25 ***** 35 COMMON /NUMSCF/ NSCF 36 COMMON /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR), 37 1 LOCDEP(MAXPAR) 38 COMMON /GRADNT/ GRAD(MAXPAR),GNFINA 39 COMMON /TIMDMP/ TLEFT, TDUMP 40 COMMON /GEOKST/ NATOMS,LABELS(NUMATM), 41 1 NA(NUMATM),NB(NUMATM),NC(NUMATM) 42 COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 43 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 44 2 NCLOSE,NOPEN,NDUMY,FRACT 45 COMMON /NUMCAL/ NUMCAL 46 COMMON /SIGMA1/ GNEXT, AMIN, ANEXT 47 COMMON /SIGMA2/ GNEXT1(MAXPAR), GMIN1(MAXPAR) 48 COMMON /NLLCOM/ HESS(MAXPAR,MAXPAR),BMAT(MAXPAR,MAXPAR), 49 1PMAT(MAXPAR*MAXPAR) 50 COMMON /SCRACH/ PVEC 51 DIMENSION IPOW(9), SIG(MAXPAR), 52 1 E1(MAXPAR), E2(MAXPAR), 53 2 P(MAXPAR), WORK(MAXPAR), 54 3 PVEC(MAXPAR*MAXPAR), EIG(MAXPAR), Q(MAXPAR) 55 LOGICAL DEBUG, RESTRT, TIMES, OKF, SCF1, RESFIL, LOG 56 CHARACTER*241 KEYWRD 57 SAVE ICALCN 58 DATA ICALCN /0/ 59 IF(ICALCN.NE.NUMCAL) THEN 60 ICALCN=NUMCAL 61 RESTRT=(INDEX(KEYWRD,'RESTART') .NE. 0) 62 LOG=(INDEX(KEYWRD,'NOLOG') .EQ. 0) 63 SCF1=(INDEX(KEYWRD,'1SCF') .NE. 0) 64 TIME1=SECOND() 65 TIME2=TIME1 66 ICYC=0 67 TIMES=(INDEX(KEYWRD,'TIME') .NE. 0) 68 TLAST=TLEFT 69 RESFIL=.FALSE. 70 LAST=0 71 ILOOP=1 72 XINC=0.00529167D0 73 RHO2=1.D-4 74 TOL2=4.D-1 75 IF(INDEX(KEYWRD,'PREC') .NE. 0) TOL2=1.D-2 76 IF(INDEX(KEYWRD,'GNORM') .NE. 0) THEN 77 TOL2=READA(KEYWRD,INDEX(KEYWRD,'GNORM')) 78 IF(TOL2.LT.0.01.AND.INDEX(KEYWRD,' LET').EQ.0)THEN 79 WRITE(6,'(/,A)')' GNORM HAS BEEN SET TOO LOW, RESET TO 0 80 1.01' 81 TOL2=0.01D0 82 ENDIF 83 ENDIF 84 DEBUG = (INDEX(KEYWRD,'POWSQ') .NE. 0) 85 IF(RESTRT) THEN 86C 87C RESTORE STORED DATA 88C 89 IPOW(9)=0 90 CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,ILOOP,BMAT,IPOW) 91 IF(SCF1) GOTO 390 92 NSCF=IPOW(8) 93 DO 10 I=1,NVAR 94 GRAD(I)=GMIN1(I) 95 10 GNEXT1(I)=GMIN1(I) 96 WRITE(6,'('' XPARAM'',6F10.6)')(XPARAM(I),I=1,NVAR) 97 IF(ILOOP .GT. 0) THEN 98C# ILOOP=ILOOP+1 99 WRITE(6,'(//10X,'' RESTARTING AT POINT'',I3)')ILOOP 100 ELSE 101 WRITE(6,'(//10X,''RESTARTING IN OPTIMISATION'', 102 1 '' ROUTINES'')') 103 ENDIF 104 ENDIF 105* 106* DEFINITIONS: NVAR = NUMBER OF GEOMETRIC VARIABLES = 3*NUMAT-6 107* 108 ENDIF 109 NVAR=ABS(NVAR) 110 IF(DEBUG) THEN 111 WRITE(6,'('' XPARAM'')') 112 WRITE(6,'(5(2I3,F10.4))')(LOC(1,I),LOC(2,I),XPARAM(I),I=1,NVAR) 113 ENDIF 114 IF( .NOT. RESTRT) THEN 115 DO 20 I=1,NVAR 116 20 GRAD(I)=0.D0 117 CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.) 118 ENDIF 119 IF(DEBUG) THEN 120 WRITE(6,'('' STARTING GRADIENTS'')') 121 WRITE(6,'(3X,8F9.4)')(GRAD(I),I=1,NVAR) 122 ENDIF 123 GMIN=SQRT(DOT(GRAD,GRAD,NVAR)) 124 DO 30 I=1,NVAR 125 GNEXT1(I)=GRAD(I) 126 GMIN1(I)=GNEXT1(I) 127 30 CONTINUE 128C 129C NOW TO CALCULATE THE HESSIAN MATRIX. 130C 131 IF(ILOOP.LT.0) GOTO 140 132C 133C CHECK THAT HESSIAN HAS NOT ALREADY BEEN CALCULATED. 134C 135 ILPR=ILOOP 136 DO 50 ILOOP=ILPR,NVAR 137 TIME1=SECOND() 138 XPARAM(ILOOP)=XPARAM(ILOOP) + XINC 139 CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.) 140 IF(SCF1) GOTO 390 141 IF(DEBUG)WRITE(6,'(I3,12(8F9.4,/3X))') 142 1 ILOOP,(GRAD(IF),IF=1,NVAR) 143 GRAD(ILOOP)=GRAD(ILOOP)+1.D-5 144 XPARAM(ILOOP)=XPARAM(ILOOP) - XINC 145 DO 40 J=1,NVAR 146 40 HESS(ILOOP,J)=-(GRAD(J)-GNEXT1(J))/XINC 147 TIME2=SECOND() 148 TSTEP=TIME2-TIME1 149 IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)') 150 1 TSTEP, TLEFT 151 IF(TLAST-TLEFT.GT.TDUMP)THEN 152 TLAST=TLEFT 153 RESFIL=.TRUE. 154 IPOW(9)=2 155 I=ILOOP 156 IPOW(8)=NSCF 157 CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) 158 ENDIF 159 IF( TLEFT .LT. TSTEP*2.D0) THEN 160C 161C STORE RESULTS TO DATE. 162C 163 IPOW(9)=1 164 I=ILOOP 165 IPOW(8)=NSCF 166 CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) 167 STOP 168 ENDIF 169 50 CONTINUE 170C ***** SCALE -HESSIAN- MATRIX ***** 171 IF( DEBUG) THEN 172 WRITE(6,'(//10X,''UN-NORMALIZED HESSIAN MATRIX'')') 173 DO 60 I=1,NVAR 174 60 WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR) 175 ENDIF 176 DO 80 I=1,NVAR 177 SUM = 0.0D0 178 DO 70 J=1,NVAR 179 70 SUM = SUM+HESS(I,J)**2 180 80 WORK(I) = 1.0D0/SQRT(SUM) 181 DO 90 I=1,NVAR 182 DO 90 J=1,NVAR 183 90 HESS(I,J) = HESS(I,J)*WORK(I) 184 IF( DEBUG) THEN 185 WRITE(6,'(//10X,''HESSIAN MATRIX'')') 186 DO 100 I=1,NVAR 187 100 WRITE(6,'(8F10.4)')(HESS(J,I),J=1,NVAR) 188 ENDIF 189C ***** INITIALIZE B MATIRX ***** 190 DO 120 I=1,NVAR 191 DO 110 J=1,NVAR 192 110 BMAT(I,J) = 0.0D0 193 120 BMAT(I,I) = WORK(I)*2.D0 194************************************************************************ 195* 196* THIS IS THE START OF THE BIG LOOP TO OPTIMIZE THE GEOMETRY 197* 198************************************************************************ 199 ILOOP=-99 200 TSTEP=TSTEP*4 201 130 CONTINUE 202 IF(TLAST-TLEFT.GT.TDUMP)THEN 203 TLAST=TLEFT 204 RESFIL=.TRUE. 205 IPOW(9)=2 206 I=ILOOP 207 IPOW(8)=NSCF 208 CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) 209 ENDIF 210 IF( TLEFT .LT. TSTEP*2.D0) THEN 211C 212C STORE RESULTS TO DATE. 213C 214 IPOW(9)=1 215 I=ILOOP 216 IPOW(8)=NSCF 217 CALL POWSAV(HESS,GMIN1,XPARAM,PMAT,I,BMAT,IPOW) 218 IFLEPO=-1 219 RETURN 220 ENDIF 221 140 CONTINUE 222C ***** FORM-A- DAGGER-A- IN PA SLONG WITH -P- ***** 223 IJ=0 224 DO 160 J=1,NVAR 225 DO 160 I=1,J 226 IJ=IJ+1 227 SUM = 0.0D0 228 DO 150 K=1,NVAR 229 150 SUM = SUM + HESS(I,K)*HESS(J,K) 230 160 PMAT(IJ) = SUM 231 DO 180 I=1,NVAR 232 SUM = 0.0D0 233 DO 170 K=1,NVAR 234 170 SUM = SUM-HESS(I,K)*GMIN1(K) 235 180 P(I) = -SUM 236 L=0 237 IF(DEBUG) THEN 238 WRITE(6,'(/10X,''P MATRIX IN POWSQ'')') 239 CALL VECPRT(PMAT,NVAR) 240 ENDIF 241 CALL RSP(PMAT,NVAR,NVAR,EIG,PVEC) 242C ***** CHECK FOR ZERO EIGENVALUE ***** 243C# WRITE(6,'('' EIGS IN POWSQ:'')') 244C# WRITE(6,'(6F13.8)')(EIG(I),I=1,NVAR) 245 IF(EIG(1).LT.RHO2) GO TO 240 246C ***** IF MATRIX IS NOT SINGULAR FORM INVERSE ***** 247C ***** BY BACK TRANSFORMING THE EIGENVECTORS ***** 248 IJ=0 249 DO 200 I=1,NVAR 250 DO 200 J=1,I 251 IJ=IJ+1 252 SUM = 0.0D0 253 DO 190 K=1,NVAR 254 190 SUM = SUM+PVEC((K-1)*NVAR+J)*PVEC((K-1)*NVAR+I)/EIG(K) 255 200 PMAT(IJ) = SUM 256C ***** FIND -Q- VECTOR ***** 257 L=0 258 IL=L+1 259 L=IL+I-1 260 DO 230 I=1,NVAR 261 SUM = 0.0D0 262 DO 210 K=1,I 263 IK=(I*(I-1))/2+K 264 210 SUM = SUM+PMAT(IK)*P(K) 265 IP1=I+1 266 DO 220 K=IP1,NVAR 267 IK=(K*(K-1))/2+I 268 220 SUM=SUM+PMAT(IK)*P(K) 269 230 Q(I) = SUM 270 GO TO 260 271 240 CONTINUE 272C ***** TAKE -Q- VECTOR AS EIGENVECTOR OF ZERO ***** 273C ***** EIGENVALUE ***** 274 DO 250 I=1,NVAR 275 250 Q(I) = PVEC(I) 276 260 CONTINUE 277C ***** FIND SEARCH DIRECTION ***** 278 DO 270 I=1,NVAR 279 SIG(I) = 0.0D0 280 DO 270 J=1,NVAR 281 270 SIG(I) = SIG(I) + Q(J)*BMAT(I,J) 282C ***** DO A ONE DIMENSIONAL SEARCH ***** 283 IF (DEBUG) THEN 284 WRITE(6,'('' SEARCH VECTOR'')') 285 WRITE(6,'(8F10.5)')(SIG(I),I=1,NVAR) 286 ENDIF 287 CALL SEARCH(XPARAM, ALPHA, SIG, NVAR, GMIN, OKF, FUNCT) 288 IF( NVAR .EQ. 1) GOTO 390 289C 290C FIRST WE ATTEMPT TO OPTIMIZE GEOMETRY USING SEARCH. 291C IF THIS DOES NOT WORK, THEN SWITCH TO LINMIN, WHICH ALWAYS WORKS, 292C BUT IS TWICE AS SLOW AS SEARCH. 293C 294 RMX = 0.0D0 295 DO 280 K=1,NVAR 296 RT = ABS(GMIN1(K)) 297 IF(RT.GT.RMX)RMX = RT 298 280 CONTINUE 299 IF(RMX.LT.TOL2) GO TO 390 300C ***** TWO STEP ESTIMATION OF DERIVATIVES ***** 301 DO 290 K=1,NVAR 302 290 E1(K) = (GMIN1(K)-GNEXT1(K))/(AMIN-ANEXT) 303 RMU = DOT(E1,GMIN1,NVAR)/DOT(GMIN1,GMIN1,NVAR) 304 DO 300 K=1,NVAR 305 300 E2(K) = E1(K) - RMU*GMIN1(K) 306C ***** SCALE -E2- AND -SIG- ***** 307 SK = 1.0D0/SQRT(DOT(E2,E2,NVAR)) 308 DO 310 K=1,NVAR 309 310 SIG(K) = SK*SIG(K) 310 DO 320 K=1,NVAR 311 320 E2(K) = SK*E2(K) 312C ***** FIND INDEX OF REPLACEMENT DIRECTION ***** 313 PMAX = -1.0D+20 314 DO 330 I=1,NVAR 315 IF(ABS(P(I)*Q(I)).LE.PMAX) GO TO 330 316 PMAX = ABS(P(I)*Q(I)) 317 ID = I 318 330 CONTINUE 319C ***** REPLACE APPROPRIATE DIRECTION AND DERIVATIVE *** 320 DO 340 K=1,NVAR 321 340 HESS(ID,K) = -E2(K) 322C ***** REPLACE STARTING POINT ***** 323 DO 350 K=1,NVAR 324 350 BMAT(K,ID) = SIG(K)/0.529167D0 325 DO 360 K=1,NVAR 326 360 GNEXT1(K) = GMIN1(K) 327 TIME1=TIME2 328 TIME2=SECOND() 329 TLEFT=TLEFT-TIME2+TIME0 330 TSTEP=TIME2-TIME1 331 ICYC=ICYC+1 332 IF(RESFIL)THEN 333 WRITE(6,370)MIN(TLEFT,9999999.9D0), 334 1MIN(GMIN,999999.999D0),FUNCT 335 IF(LOG)WRITE(11,370)MIN(TLEFT,9999999.9D0), 336 1MIN(GMIN,999999.999D0),FUNCT 337 370 FORMAT(' RESTART FILE WRITTEN, TIME LEFT:',F9.1, 338 1' GRAD.:',F10.3,' HEAT:',G14.7) 339 RESFIL=.FALSE. 340 ELSE 341 WRITE(6,380)ICYC,MIN(TSTEP,9999.99D0), 342 1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT 343 IF(LOG)WRITE(11,380)ICYC,MIN(TSTEP,9999.99D0), 344 1MIN(TLEFT,9999999.9D0),MIN(GMIN,999999.999D0),FUNCT 345 380 FORMAT(' CYCLE:',I5,' TIME:',F6.1,' TIME LEFT:',F9.1, 346 1' GRAD.:',F10.3,' HEAT:',G14.7) 347 ENDIF 348 IF(TIMES)WRITE(6,'('' TIME FOR STEP:'',F8.2,'' LEFT'',F8.2)') 349 1TSTEP, TLEFT 350 GO TO 130 351 390 CONTINUE 352 DO 400 I=1,NVAR 353 400 GRAD(I)=0.D0 354 LAST=1 355 CALL COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.) 356 DO 410 I=1,NVAR 357 410 GRAD(I)=GMIN1(I) 358 GNFINA=SQRT(DOT(GRAD,GRAD,NVAR)) 359 IFLEPO=11 360 IF(SCF1)IFLEPO=13 361 RETURN 362 END 363