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