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
|
SUBROUTINE FMAT(FMATRX, NREAL, TSCF, TDER, DELDIP, HEAT)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
COMMON /SYMOPS/ R(14,120), NSYM, IPO(NUMATM,120), NENT
DIMENSION FMATRX(*), DELDIP(3,*)
***********************************************************************
*
* VALUE CALCULATES THE SECOND-ORDER OF THE ENERGY WITH
* RESPECT TO THE CARTESIAN COORDINATES I AND J AND PLACES IT
* IN FMATRX
*
* ON INPUT NATOMS = NUMBER OF ATOMS IN THE SYSTEM.
* XPARAM = INTERNAL COORDINATES OF MOLECULE STORED LINEARLY
*
* VARIABLES USED
* COORDL = ARRAY OF CARTESIAN COORDINATES, STORED LINEARLY.
* I = INDEX OF CARTESIAN COORDINATE.
* J = INDEX OF CARTESIAN COORDINATE.
*
* ON OUTPUT FMATRX = SECOND DERIVATIVE OF THE ENERGY WITH RESPECT TO
* CARTESIAN COORDINATES I AND J.
***********************************************************************
COMMON /KEYWRD/ KEYWRD
COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
1 NA(NUMATM),NB(NUMATM),NC(NUMATM)
COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY, DUMY(MAXPAR)
COMMON /DENSTY/ P(MPACK),PDUMY(2,MPACK)
COMMON /TIMDMP/ TLEFT, TDUMP
COMMON /ATMASS/ ATMASS(NUMATM)
C ***** Modified by Jiro Toyoda at 1994-05-25 *****
C COMMON /TIME / TIME0
COMMON /TIMEC / TIME0
C ***************************** at 1994-05-25 *****
COMMON /CORE / CORE(107)
COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
2 NCLOSE,NOPEN,NDUMY,FRACT
COMMON /COORD / COORD(3,NUMATM)
COMMON /NLLCOM/ EVECS(MAXPAR*MAXPAR),BMAT(MAXPAR,MAXPAR*2)
DIMENSION GRAD(MAXPAR),
1GROLD(MAXPAR), COORDL(MAXPAR), Q(NUMATM), DEL2(3), G2OLD(MAXPAR)
2, EIGS(MAXPAR), G2RAD(MAXPAR),
3 FCONST(MAXPAR)
CHARACTER*241 KEYWRD
SAVE FACT
LOGICAL DEBUG, RESTRT, PRNT, RESFIL, PRECIS, BIG, LOG, GROUP
EQUIVALENCE (COORD(1,1),COORDL(1))
DATA FACT/6.95125D-3/
C
C FACT IS THE CONVERSION FACTOR FROM KCAL/MOLE TO ERGS
C
C SET UP CONSTANTS AND FLAGS
NA(1)=99
C
C SET UP THE VARIABLES IN XPARAM ANDLOC,THESE ARE IN CARTESIAN COORDINA
C
NUMAT=0
C$DOIT ASIS
DO 10 I=1,NATOMS
IF(LABELS(I).NE.99.AND.LABELS(I).NE.107) THEN
NUMAT=NUMAT+1
LABELS(NUMAT)=LABELS(I)
ENDIF
10 CONTINUE
NATOMS=NUMAT
C
C THIS IS A QUICK, IF CLUMSY, WAY TO CALCULATE NUMAT, AND TO REMOVE
C THE DUMMY ATOMS FROM THE ARRAY LABELS.
C
NVAR=NUMAT*3
DO 20 I=1,NUMAT
LOC(1,(I-1)*3+1)=I
LOC(2,(I-1)*3+1)=1
C
LOC(1,(I-1)*3+2)=I
LOC(2,(I-1)*3+2)=2
C
LOC(1,(I-1)*3+3)=I
LOC(2,(I-1)*3+3)=3
20 CONTINUE
LIN=(NVAR*(NVAR+1))/2
DO 30 I=1,LIN
30 FMATRX(I)=0.D0
PRNT =(INDEX(KEYWRD,'IRC=') .EQ. 0)
LOG =(INDEX(KEYWRD,'NOLOG') .EQ. 0)
PRECIS =(INDEX(KEYWRD,'PREC') .NE. 0)
RESTRT =(INDEX(KEYWRD,'RESTART') .NE. 0)
GROUP =(INDEX(KEYWRD,' GROUP').NE.0)
IF(INDEX(KEYWRD,'NLLSQ') .NE. 0) RESTRT=.FALSE.
DEBUG =(INDEX(KEYWRD,'FMAT') .NE. 0)
BIG =(INDEX(KEYWRD,'LARGE') .NE. 0 .AND. DEBUG)
IF(PRNT)WRITE(6,'(//4X,''FIRST DERIVATIVES WILL BE USED IN THE''
1,'' CALCULATION OF SECOND DERIVATIVES'')')
TLAST=TLEFT
RESFIL=.FALSE.
IF(RESTRT) THEN
ISTART = 0
I=0
CALL FORSAV(TOTIME,DELDIP,ISTART,FMATRX, COORD, NVAR,HEAT,
1 EVECS,JSTART,FCONST)
KOUNTF=(ISTART*(ISTART+1))/2
ISTART=ISTART+1
JSTART=JSTART+1
TIME2 = SECOND()
ELSE
KOUNTF=0
TOTIME=0.D0
IF (TSCF.GT.0.D0)TLEFT=TLEFT-TSCF-TDER
ISTART=1
ENDIF
C CALCULATE FMATRX
IF(ISTART.GT.1) THEN
ESTIME=(NVAR-ISTART+1)*TOTIME/(ISTART-1.D0)
ELSE
ESTIME=NVAR*(TSCF+TDER)*2.D0
IF (PRECIS) ESTIME=ESTIME*2.D0
ENDIF
IF(TSCF.GT.0)
1WRITE(6,'(/10X,''ESTIMATED TIME TO COMPLETE CALCULATION =''
2,F9.2,'' SECONDS'')')ESTIME
IF(RESTRT) THEN
IF(ISTART.LE.NVAR)
1 WRITE(6,'(/10X,''STARTING AGAIN AT LINE'',18X,I4)')ISTART
WRITE(6,'(/10X,''TIME USED UP TO RESTART ='',F22.2)')TOTIME
ENDIF
LU=KOUNTF
NUMAT=NVAR/3
DO 40 I=1,NVAR
40 EIGS(I)=0.D0
C
C READ IN THE SYMMETRY OPERATIONS, IF PRESENT
C
IF(GROUP) CALL SYMR
ISKIP=0
DO 110 I=ISTART,NVAR
IF(GROUP .AND. ((I-1)/3)*3.EQ.I-1)THEN
C
C START OF A NEW ATOM. DOES A SYMMETRY OPERATION RELATE AN ALREADY
C CALCULATED ATOM TO THIS ONE
C
J=(I+2)/3
CALL SYMPOP(FMATRX, J, ISKIP, DELDIP)
ENDIF
IF(ISKIP.GT.0) THEN
WRITE(6,'('' STEP:'',I4,'' '',9X, '' INTEGRAL =
1'',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
ISKIP=ISKIP-1
LU=LU+I
GOTO 110
ENDIF
TIME2 = SECOND()
DELTA=1.D0/120.D0
IF(PRECIS)THEN
C
C DETERMINE A GOOD STEP SIZE
C
G2OLD(1)=100.D0
COORDL(I)=COORDL(I)+DELTA
CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.)
COORDL(I)=COORDL(I)-DELTA
DELTA=DELTA*10.D0/SQRT(DOT(G2OLD,G2OLD,NVAR))
C
C CONSTRAIN DELTA TO A 'REASONABLE' VALUE
C
DELTA=MIN(0.05D0,MAX(0.005D0,DELTA))
IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT
1A
G2OLD(1)=100.D0
COORDL(I)=COORDL(I)+DELTA
CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.)
IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +1.0*DELTA',
1SQRT(DOT(G2OLD,G2OLD,NVAR))
COORDL(I)=COORDL(I)-DELTA*2.D0
G2RAD(1)=100.D0
CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., G2RAD, .TRUE.)
COORDL(I)=COORDL(I)+DELTA
IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -1.0*DELTA',
1SQRT(DOT(G2RAD,G2RAD,NVAR))
ELSE
IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT
1A
ENDIF
COORDL(I)=COORDL(I)+0.5D0*DELTA
GROLD(1)=100.D0
CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., GROLD, .TRUE.)
IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +0.5*DELTA',
1SQRT(DOT(GROLD,GROLD,NVAR))
CALL CHRGE(P,Q)
DO 50 II=1,NUMAT
50 Q(II)=CORE(LABELS(II))-Q(II)
SUM = DIPOLE(P,Q,COORDL,DELDIP(1,I),0)
COORDL(I)=COORDL(I)-DELTA
GRAD(1)=100.D0
CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., GRAD, .TRUE.)
IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -0.5*DELTA',
1SQRT(DOT(GRAD,GRAD,NVAR))
CALL CHRGE(P,Q)
DO 60 II=1,NUMAT
60 Q(II)=CORE(LABELS(II))-Q(II)
SUM = DIPOLE(P,Q,COORDL,DEL2,0)
COORDL(I)=COORDL(I)+DELTA*0.5D0
DELDIP(1,I)=(DELDIP(1,I)-DEL2(1))*0.5D0/DELTA
DELDIP(2,I)=(DELDIP(2,I)-DEL2(2))*0.5D0/DELTA
DELDIP(3,I)=(DELDIP(3,I)-DEL2(3))*0.5D0/DELTA
LL=LU+1
LU=LL+I-1
L=0
IF(PRECIS)THEN
DO 70 KOUNTF=LL,LU
L=L+1
C
C G2OLD = X + 1.0*DELTA
C GROLD = X + 0.5*DELTA
C GRAD = X - 0.5*DELTA
C G2RAD = X - 1.0*DELTA
C
DUMY(L)= (8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
1 /DELTA*FACT/24.D0
EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
1 /DELTA**3*FACT/56.D0
C
C CORRECT FOR 4'TH ORDER CONTAMINATION
C
C# CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0)
C# DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L))
FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L)
70 CONTINUE
L=L-1
DO 80 K=I,NVAR
L=L+1
KK=(K*(K-1))/2+I
DUMY(L)=(8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
1 /DELTA*FACT/24.D0
EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L)))
1 /DELTA**3*FACT/56.D0
C
C CORRECT FOR 4'TH ORDER CONTAMINATION
C
C# CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0)
C# DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L))
FMATRX(KK)=FMATRX(KK)+DUMY(L)
80 CONTINUE
ELSE
DO 90 KOUNTF=LL,LU
L=L+1
DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT
FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L)
90 CONTINUE
L=L-1
DO 100 K=I,NVAR
L=L+1
KK=(K*(K-1))/2+I
DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT
FMATRX(KK)=FMATRX(KK)+DUMY(L)
100 CONTINUE
ENDIF
IF(BIG)THEN
WRITE(6,'(A)')' CONTRIBUTIONS TO F-MATRIX'
WRITE(6,'(A)')' ELEMENT +1.0*DELTA +0.5*DELTA -0.5*DEL'
1//'TA -1.0*DELTA 2''ND ORDER 4TH ORDER'
WRITE(6,'(I7,6F12.6)')(L,G2OLD(L),GROLD(L),GRAD(L),G2RAD(L),
1DUMY(L),EIGS(L),L=1,NVAR)
ENDIF
TIME3 = SECOND()
TSTEP=TIME3-TIME2
TLEFT= MAX(0.1D0,TLEFT-TSTEP)
IF(TSTEP.GT.1.D6)TSTEP=TSTEP-1.D6
TOTIME= TOTIME+TSTEP
IF(RESFIL)THEN
WRITE(6,'('' STEP:'',I4,'' RESTART FILE WRITTEN, INTEGRAL ='
1',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
IF(LOG)WRITE(11,'('' STEP:'',I4,'' RESTART FILE WRITTEN, '',
1''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT
RESFIL=.FALSE.
ELSE
WRITE(6,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL =
1'',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT
IF(LOG) WRITE(11,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, '
1',''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEF
2T
ENDIF
ESTIM = TOTIME/I
IF(TLAST-TLEFT.GT.TDUMP)THEN
TLAST=TLEFT
RESFIL=.TRUE.
JSTART=1
II=I
CALL FORSAV(TOTIME,DELDIP,II,FMATRX, COORD,NVAR,HEAT,
1 EVECS,JSTART,FCONST)
ENDIF
IF(I.NE.NVAR.AND.TLEFT-10.D0 .LT. ESTIM) THEN
WRITE(6,'(//10X,''- - - - - - - TIME UP - - - - - - -'',//)'
1)
WRITE(6,'(/10X,'' POINT REACHED ='',I4)')I
WRITE(6,'(/10X,'' RESTART USING KEY-WORD "RESTART"'')')
WRITE(6,'(10X,''ESTIMATED TIME FOR THE NEXT STEP ='',F8.2,
1'' SECONDS'')')ESTIM
JSTART=1
II=I
CALL FORSAV(TOTIME,DELDIP,II,FMATRX, COORD,NVAR,HEAT,
1 EVECS,JSTART,FCONST)
WRITE(6,'(//10X,''FORCE MATRIX WRITTEN TO DISK'')')
NREAL=-1
RETURN
ENDIF
110 CONTINUE
DO 120 I=1,NATOMS
IF(ATMASS(I).LT.1.D-20.AND.LABELS(I).LT.99)THEN
CALL FORSAV(TOTIME,DELDIP,NVAR,FMATRX, COORD,NVAR,HEAT,
1 EVECS,ILOOP,FCONST)
WRITE(6,'(A)')' AT LEAST ONE ATOM HAS A ZERO MASS. A RESTART
1'
WRITE(6,'(A)')' FILE HAS BEEN WRITTEN AND THE JOB STOPPED'
STOP
ENDIF
120 CONTINUE
IF(ISTART.LE.NVAR .AND. INDEX(KEYWRD,'ISOT') .NE. 0)
1CALL FORSAV(TOTIME,DELDIP,NVAR,FMATRX, COORD,NVAR,HEAT,
2 EVECS,ILOOP,FCONST)
RETURN
END
|
