1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
|
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), XCOORD(3,NUMATM)
COMMON /LAST / LAST
COMMON /KEYWRD/ KEYWRD
C ***** Modified by Jiro Toyoda at 1994-05-25 *****
C COMMON /TIME / TIME0
COMMON /TIMEC / TIME0
C ***************************** at 1994-05-25 *****
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, 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)
TIME1=SECOND()
TIME2=TIME1
ICYC=0
TIMES=(INDEX(KEYWRD,'TIME') .NE. 0)
TLAST=TLEFT
RESFIL=.FALSE.
LAST=0
ILOOP=1
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))
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
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
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)
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
|
