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
363
364
365
|
SUBROUTINE DERI1(C,NORBS,COORD,NUMBER,WORK,GRAD
1 ,F,MINEAR,FD,WMAT,HMAT,FMAT)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
DIMENSION WMAT(MPACK),HMAT(MPACK*2),FMAT(MPACK*2)
*********************************************************************
*
* DERI1 COMPUTE THE NON-RELAXED DERIVATIVE OF THE NON-VARIATIONALLY
* OPTIMIZED WAVEFUNCTION ENERGY WITH RESPECT TO ONE CARTESIAN
* COORDINATE AT A TIME
* AND
* COMPUTE THE NON-RELAXED FOCK MATRIX DERIVATIVE IN M.O BASIS AS
* REQUIRED IN THE RELAXATION SECTION (ROUTINE 'DERI2').
*
* INPUT
* C(NORBS,NORBS) : M.O. COEFFICIENTS.
* COORD : CARTESIAN COORDINATES ARRAY.
* NUMBER : LOCATION OF THE REQUIRED VARIABLE IN COORD.
* WORK : WORK ARRAY OF SIZE N*N.
* WMAT : WORK ARRAYS FOR d<PQ|RS> (2-CENTERS A.O)
* OUTPUT
* C,COORD,NUMBER : NOT MODIFIED.
* GRAD : DERIVATIVE OF THE HEAT OF FORMATION WITH RESPECT TO
* COORD(NUMBER), WITHOUT RELAXATION CORRECTION.
* F(MINEAR) : NON-RELAXED FOCK MATRIX DERIVATIVE WITH RESPECT TO
* COORD(NUMBER), EXPRESSED IN M.O BASIS, SCALED AND
* PACKED, OFF-DIAGONAL BLOCKS ONLY.
* FD : IDEM BUT UNSCALED, DIAGONAL BLOCKS, C.I-ACTIVE ONLY.
*
************************************************************************
COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM)
1 ,NLAST(NUMATM), NDUMY1, NELECS,NALPHA,NBETA
2 ,NCLOSE,NOPEN,NDUMY,FRACT
3 /VECTOR/ CDUM(MORB2),EIGS(MAXORB),WDUM(MORB2),EIGB(MAXORB)
4 /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK)
COMMON /CIBITS/ NMOS,LAB,NELEC,NBO(3)
1 /HMATRX/ H(MPACK)
2 /XYIJKL/ XY(NMECI,NMECI,NMECI,NMECI)
3 /CIVECT/ CONF(NMECI**4+NMECI**2+1)
COMMON /FOKMAT/ FDUMY(MPACK), SCALAR(MPACK)
COMMON /NVOMAT/ DIAG(MPACK/2)
1 /KEYWRD/ KEYWRD
COMMON /NUMCAL/ NUMCAL
DIMENSION COORD(*),C(NORBS,NORBS),WORK(NORBS,NORBS),F(*),FD(*)
CHARACTER KEYWRD*241
LOGICAL DEBUG
DATA ICALCN /0/
C
IF(ICALCN.NE.NUMCAL) THEN
DEBUG=INDEX(KEYWRD,'DERI1').NE.0
IPRT=6
LINEAR=NORBS*(NORBS+1)/2
ICALCN=NUMCAL
ENDIF
IF(DEBUG) CALL TIMER('BEFORE DERI1')
STEP=1.D-3
C
C 2 POINTS FINITE DIFFERENCE TO GET THE INTEGRAL DERIVATIVES
C ----------------------------------------------------------
C STORED IN HMAT AND WMAT, WITHOUT DIVIDING BY THE STEP SIZE.
C
NATI=(NUMBER-1)/3+1
NATX=NUMBER-3*(NATI-1)
CALL DHCORE (COORD,HMAT,WMAT,ENUCL2,NATI,NATX,STEP)
C
C HMAT HOLDS THE ONE-ELECTRON DERIVATIVES OF ATOM NATI FOR DIRECTION
C NATX W.R.T. ALL OTHER ATOMS
C WMAT HOLDS THE TWO-ELECTRON DERIVATIVES OF ATOM NATI FOR DIRECTION
C NATX W.R.T. ALL OTHER ATOMS
STEP=0.5D0/STEP
C
C NON-RELAXED FOCK MATRIX DERIVATIVE IN A.O BASIS.
C ------------------------------------------------
C STORED IN FMAT, DIVIDED BY STEP.
C
CALL SCOPY (LINEAR,HMAT,1,FMAT,1)
CALL DFOCK2 (FMAT,P,PA,WMAT,NUMAT,NFIRST,NMIDLE,NLAST,NATI)
C
C FMAT HOLDS THE ONE PLUS TWO - ELECTRON DERIVATIVES OF ATOM NATI FOR
C DIRECTION NATX W.R.T. ALL OTHER ATOMS
C
C DERIVATIVE OF THE SCF-ONLY ENERGY (I.E BEFORE C.I CORRECTION)
C
GRAD=(HELECT(NORBS,P,HMAT,FMAT)+ENUCL2)*STEP
C TAKE STEP INTO ACCOUNT IN FMAT
DO 10 I=1,LINEAR
10 FMAT(I)=FMAT(I)*STEP
C
C RIGHT-HAND SIDE SUPER-VECTOR F = C' FMAT C USED IN RELAXATION
C -----------------------------------------------------------
C STORED IN NON-STANDARD PACKED FORM IN F(MINEAR) AND FD.
C THE SUPERVECTOR IS THE NON-RELAXED FOCK MATRIX DERIVATIVE IN
C M.O BASIS: F(IJ)= ( (C' * FOCK * C)(I,J) ) WITH I.GT.J .
C F IS SCALED AND PACKED IN SUPERVECTOR FORM WITH
C THE CONSECUTIVE FOLLOWING OFF-DIAGONAL BLOCKS:
C 1) OPEN-CLOSED I.E. F(IJ)=F(I,J) WITH I OPEN & J CLOSED
C AND I RUNNING FASTER THAN J,
C 2) VIRTUAL-CLOSED SAME RULE OF ORDERING,
C 3) VIRTUAL-OPEN SAME RULE OF ORDERING.
C FD IS PACKED OVER THE C.I-ACTIVE M.O WITH
C THE CONSECUTIVE DIAGONAL BLOCKS:
C 1) CLOSED-CLOSED IN CANONICAL ORDER, WITHOUT THE
C DIAGONAL ELEMENTS,
C 2) OPEN-OPEN SAME RULE OF ORDERING,
C 3) VIRTUAL-VIRTUAL SAME RULE OF ORDERING.
C
C PART 1 : WORK(N,N) = FMAT(N,N) * C(N,N)
DO 20 I=1,NORBS
20 CALL SUPDOT (WORK(1,I),FMAT,C(1,I),NORBS,1)
C
C PART 2 : F(IJ) = (C' * WORK)(I,J) ... OFF-DIAGONAL BLOCKS.
L=1
IF(NBO(2).NE.0 .AND. NBO(1).NE.0) THEN
C OPEN-CLOSED
CALL MTXM (C(1,NBO(1)+1),NBO(2),WORK,NORBS,F(L),NBO(1))
L=L+NBO(2)*NBO(1)
ENDIF
IF(NBO(3).NE.0 .AND. NBO(1).NE.0) THEN
C VIRTUAL-CLOSED
CALL MTXM (C(1,NOPEN+1),NBO(3),WORK,NORBS,F(L),NBO(1))
L=L+NBO(3)*NBO(1)
ENDIF
IF(NBO(3).NE.0 .AND. NBO(2).NE.0) THEN
C VIRTUAL-OPEN
CALL MTXM (C(1,NOPEN+1),NBO(3),WORK(1,NBO(1)+1),NORBS,F(L),NBO(
12))
ENDIF
C SCALE F ACCORDING TO THE DIAGONAL METRIC TENSOR 'SCALAR '.
DO 30 I=1,MINEAR
30 F(I)=F(I)*SCALAR(I)
IF(DEBUG)THEN
WRITE(6,*)' F IN DERI1'
J=MIN(20,MINEAR)
WRITE(6,'(5F12.6)')(F(I),I=1,J)
ENDIF
C
C PART 3 : SUPER-VECTOR FD, C.I-ACTIVE DIAGONAL BLOCKS, UNSCALED.
L=1
NEND=0
DO 50 LOOP=1,3
NINIT=NEND+1
NEND =NEND+NBO(LOOP)
N1=MAX(NINIT,NELEC+1 )
N2=MIN(NEND ,NELEC+NMOS)
IF(N2.LT.N1) GO TO 50
DO 40 I=N1,N2
IF(I.GT.NINIT) THEN
CALL MXM (C(1,I),1,WORK(1,NINIT),NORBS,FD(L),I-NINIT)
L=L+I-NINIT
ENDIF
40 CONTINUE
50 CONTINUE
C
C NON-RELAXED C.I CORRECTION TO THE ENERGY DERIVATIVE.
C ----------------------------------------------------
C
C C.I-ACTIVE FOCK EIGENVALUES DERIVATIVES, STORED IN FD(CONTINUED).
LCUT=L
DO 60 I=NELEC+1,NELEC+NMOS
FD(L)=DOT(C(1,I),WORK(1,I),NORBS)
60 L=L+1
C
C C.I-ACTIVE 2-ELECTRONS INTEGRALS DERIVATIVES. STORED IN XY.
C FMAT IS USED HERE AS SCRATCH SPACE
C
CALL DIJKL1 (C(1,NELEC+1),NORBS,NATI,WMAT,FMAT,HMAT,FMAT)
DO 70 I=1,NMOS
DO 70 J=1,NMOS
DO 70 K=1,NMOS
DO 70 L=1,NMOS
70 XY(I,J,K,L)=XY(I,J,K,L)*STEP
C
C BUILD THE C.I MATRIX DERIVATIVE, STORED IN WMAT.
CALL MECID (FD(LCUT-NELEC),GSE,EIGB,WORK)
IF(DEBUG)THEN
WRITE(6,*)' GSE:',GSE
C# WRITE(6,*)' EIGB:',(EIGB(I),I=1,10)
C# WRITE(6,*)' WORK:',(WORK(I,1),I=1,10)
ENDIF
CALL MECIH (WORK,WMAT,NMOS,LAB)
C
C NON-RELAXED C.I CONTRIBUTION TO THE ENERGY DERIVATIVE.
CALL SUPDOT (WORK,WMAT,CONF,LAB,1)
GRAD=(GRAD+DOT(CONF,WORK,LAB))*23.061D0
IF(DEBUG) THEN
WRITE(IPRT,'('' * * * GRADIENT COMPONENT NUMBER'',I4)')NUMBER
WRITE(IPRT,'('' NON-RELAXED C.I-ACTIVE FOCK EIGENVALUES '',
1 ''DERIVATIVES (E.V.)'')')
WRITE(IPRT,'(8F10.4)')(FD(LCUT-1+I),I=1,NMOS)
WRITE(IPRT,'('' NON-RELAXED 2-ELECTRONS DERIVATIVES (E.V.)''/
1'' I J K L d<I(1)J(1)|K(2)L(2)>'')')
DO 80 I=1,NMOS
DO 80 J=1,I
DO 80 K=1,I
LL=K
IF(K.EQ.I) LL=J
DO 80 L=1,LL
80 WRITE(IPRT,'(4I5,F20.10)')
1 NELEC+I,NELEC+J,NELEC+K,NELEC+L,XY(I,J,K,L)
WRITE(IPRT,'('' NON-RELAXED GRADIENT COMPONENT'',F10.4,
1'' KCAL/MOLE'')')GRAD
CALL TIMER('AFTER DERI1')
ENDIF
RETURN
END
SUBROUTINE DHCORE (COORD,H,W,ENUCLR,NATI,NATX,STEP)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'SIZES'
DIMENSION COORD(3,*),H(*),W(*)
C
C DHCORE GENERATES THE 1-ELECTRON AND 2-ELECTRON INTEGRALS DERIVATIVES
C WITH RESPECT TO THE CARTESIAN COORDINATE COORD (NATX,NATI).
C
C INPUT
C COORD : CARTESIAN COORDINATES OF THE MOLECULE.
C NATI,NATX : INDICES OF THE MOVING COORDINATE.
C STEP : STEP SIZE OF THE 2-POINTS FINITE DIFFERENCE.
C OUTPUT
C H : 1-ELECTRON INTEGRALS DERIVATIVES (PACKED CANONICAL).
C W : 2-ELECTRON INTEGRALS DERIVATIVES (ORDERED AS REQUIRED
C IN DFOCK2 AND DIJKL1).
C ENUCLR : NUCLEAR ENERGY DERIVATIVE.
C
COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
2 NCLOSE,NOPEN,NDUMY,FRACT
3 /MOLORB/ USPD(MAXORB),DUMY(MAXORB)
4 /KEYWRD/ KEYWRD
5 /WMATRX/ WDUMMY(N2ELEC*2)
CHARACTER*241 KEYWRD
LOGICAL FIRST,MINDO
DIMENSION E1B(10),DE1B(10),E2A(10),DE2A(10)
1 ,DI(9,9),DDI(9,9),WJD(101),DWJD(101)
2,NB(0:8)
DATA NB/1,0,0,10,0,0,0,0,45/
DATA FIRST/.TRUE./
IF (FIRST) THEN
CUTOFF=1.D10
FIRST=.FALSE.
MINDO=INDEX(KEYWRD,'MINDO') .NE. 0
ENDIF
DO 10 I=1,(NORBS*(NORBS+1))/2
10 H(I)=0
ENUCLR=0.D0
KR=1
NROW=0
I=NATI
CSAVE=COORD(NATX,NATI)
IA=NFIRST(NATI)
IB=NLAST(NATI)
IC=NMIDLE(NATI)
NI=NAT(NATI)
NROW=-NB(IB-IA)
DO 20 J=1,NUMAT
20 NROW=NROW+NB(NLAST(J)-NFIRST(J))
C# NCOL=NB(NLAST(NATI)-NFIRST(NATI))
NBAND2=0
DO 120 J=1,NUMAT
IF (J.EQ.NATI) GO TO 120
JA=NFIRST(J)
JB=NLAST(J)
JC=NMIDLE(J)
NJ=NAT(J)
COORD(NATX,NATI)=CSAVE+STEP
CALL H1ELEC(NI,NJ,COORD(1,NATI),COORD(1,J),DI)
C
C THE FOLLOWING STYLE WAS NECESSARY TO GET ROUND A BUG IN THE
C GOULD COMPILER
C
COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
CALL H1ELEC(NI,NJ,COORD(1,NATI),COORD(1,J),DDI)
C
C FILL THE ATOM-OTHER ATOM ONE-ELECTRON MATRIX.
C
I2=0
IF (IA.GT.JA) THEN
DO 30 I1=IA,IB
IJ=I1*(I1-1)/2+JA-1
I2=I2+1
J2=0
DO 30 J1=JA,JB
IJ=IJ+1
J2=J2+1
30 H(IJ)=H(IJ)+(DI(I2,J2)-DDI(I2,J2))
ELSE
DO 40 I1=JA,JB
IJ=I1*(I1-1)/2+IA-1
I2=I2+1
J2=0
DO 40 J1=IA,IB
IJ=IJ+1
J2=J2+1
40 H(IJ)=H(IJ)+(DI(J2,I2)-DDI(J2,I2))
ENDIF
C
C CALCULATE THE TWO-ELECTRON INTEGRALS, W; THE ELECTRON NUCLEAR TERM
C E1B AND E2A; AND THE NUCLEAR-NUCLEAR TERM ENUC.
C
KRO=KR
NBAND2=NBAND2+NB(NLAST(J)-NFIRST(J))
IF (MINDO) THEN
COORD(NATX,NATI)=CSAVE+STEP
CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
1 ,WJD,KR,E1B,E2A,ENUC,CUTOFF)
KR=KRO
COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
1 ,DWJD,KR,DE1B,DE2A,DENUC,CUTOFF)
IF (KR.GT.KRO) THEN
DO 50 K=1,KR-KRO+1
50 W(KRO+K-1)=WJD(K)-DWJD(K)
ENDIF
ELSE
COORD(NATX,NATI)=CSAVE+STEP
CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
1 ,WJD,KR,E1B,E2A,ENUC,CUTOFF)
KR=KRO
COORD(NATX,NATI)=CSAVE+STEP*(-1.D0)
CALL ROTATE (NI,NJ,COORD(1,NATI),COORD(1,J)
1 ,DWJD,KR,DE1B,DE2A,DENUC,CUTOFF)
IF (KR.GT.KRO) THEN
DO 60 K=1,KR-KRO+1
60 WJD(K)=WJD(K)-DWJD(K)
J7=0
DO 70 I1=KRO,KR
J7=J7+1
70 W(I1)=WJD(J7)
ENDIF
ENDIF
COORD(NATX,NATI)=CSAVE
ENUCLR = ENUCLR + ENUC-DENUC
C
C ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM I.
C
I2=0
DO 80 I1=IA,IC
II=I1*(I1-1)/2+IA-1
DO 80 J1=IA,I1
II=II+1
I2=I2+1
80 H(II)=H(II)+E1B(I2)-DE1B(I2)
C CONTRIB D, CNDO.
DO 90 I1=IC+1,IB
II=(I1*(I1+1))/2
90 H(II)=H(II)+E1B(1)-DE1B(1)
C
C ADD ON THE ELECTRON-NUCLEAR ATTRACTION TERM FOR ATOM J.
C
I2=0
DO 100 I1=JA,JC
II=I1*(I1-1)/2+JA-1
DO 100 J1=JA,I1
II=II+1
I2=I2+1
100 H(II)=H(II)+E2A(I2)-DE2A(I2)
C CONTRIB D, CNDO.
DO 110 I1=JC+1,JB
II=(I1*(I1+1))/2
110 H(II)=H(II)+E2A(1)-DE2A(1)
120 CONTINUE
C
C 'SIZE' OF H IS NROW * NCOL
C
RETURN
END
|
