1 /* Copyright 2014-2018 The PySCF Developers. All Rights Reserved.
2
3 Licensed under the Apache License, Version 2.0 (the "License");
4 you may not use this file except in compliance with the License.
5 You may obtain a copy of the License at
6
7 http://www.apache.org/licenses/LICENSE-2.0
8
9 Unless required by applicable law or agreed to in writing, software
10 distributed under the License is distributed on an "AS IS" BASIS,
11 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 See the License for the specific language governing permissions and
13 limitations under the License.
14
15 *
16 * JKoperator
17 */
18
19 #include <stdlib.h>
20 #include <stdio.h>
21 #include "nr_direct.h"
22 #include "np_helper/np_helper.h"
23
24 #define ASSERT(expr, msg) \
25 if (!(expr)) { fprintf(stderr, "Fail at %s\n", msg); exit(1); }
26
27 #define MAXCGTO 64
28
29 #define ISH0 0
30 #define ISH1 1
31 #define JSH0 2
32 #define JSH1 3
33 #define KSH0 4
34 #define KSH1 5
35 #define LSH0 6
36 #define LSH1 7
37
38 #define JKOP_DATA_SIZE(obra, oket) \
39 static size_t JKOperator_data_size_##obra##oket(int *shls_slice, int *ao_loc) \
40 { \
41 int nbra = ao_loc[shls_slice[obra##SH1]] - ao_loc[shls_slice[obra##SH0]]; \
42 int nket = ao_loc[shls_slice[oket##SH1]] - ao_loc[shls_slice[oket##SH0]]; \
43 return nbra * nket; \
44 }
JKOP_DATA_SIZE(K,L)45 JKOP_DATA_SIZE(K, L)
46 JKOP_DATA_SIZE(L, K)
47 JKOP_DATA_SIZE(I, J)
48 JKOP_DATA_SIZE(J, I)
49 JKOP_DATA_SIZE(K, J)
50 JKOP_DATA_SIZE(J, K)
51 JKOP_DATA_SIZE(I, L)
52 JKOP_DATA_SIZE(L, I)
53 JKOP_DATA_SIZE(K, I)
54 JKOP_DATA_SIZE(I, K)
55 JKOP_DATA_SIZE(J, L)
56 JKOP_DATA_SIZE(L, J)
57
58 #define ADD_JKOP(fname, ibra, iket, obra, oket, type) \
59 JKOperator CVHF##fname = {ibra##SH0, iket##SH0, obra##SH0, oket##SH0, \
60 fname, JKOperator_data_size_##obra##oket, \
61 JKOperator_sanity_check_##type}
62
63 static void JKOperator_sanity_check_s1(int *shls_slice)
64 {
65 }
JKOperator_sanity_check_s2ij(int * shls_slice)66 static void JKOperator_sanity_check_s2ij(int *shls_slice)
67 {
68 ASSERT(((shls_slice[0] == shls_slice[2]) &&
69 (shls_slice[1] == shls_slice[3])), "s2ij");
70 }
JKOperator_sanity_check_s2kl(int * shls_slice)71 static void JKOperator_sanity_check_s2kl(int *shls_slice)
72 {
73 ASSERT(((shls_slice[4] == shls_slice[6]) &&
74 (shls_slice[5] == shls_slice[7])), "s2kl");
75 }
JKOperator_sanity_check_s4(int * shls_slice)76 static void JKOperator_sanity_check_s4(int *shls_slice)
77 {
78 ASSERT(((shls_slice[0] == shls_slice[2]) &&
79 (shls_slice[1] == shls_slice[3])), "s4 ij");
80 ASSERT(((shls_slice[4] == shls_slice[6]) &&
81 (shls_slice[5] == shls_slice[7])), "s4 kl");
82 }
JKOperator_sanity_check_s8(int * shls_slice)83 static void JKOperator_sanity_check_s8(int *shls_slice)
84 {
85 ASSERT(((shls_slice[0] == shls_slice[2]) &&
86 (shls_slice[1] == shls_slice[3])), "s8 ij");
87 ASSERT(((shls_slice[4] == shls_slice[6]) &&
88 (shls_slice[5] == shls_slice[7])), "s8 kl");
89 ASSERT(((shls_slice[0] == shls_slice[4]) &&
90 (shls_slice[1] == shls_slice[5])), "s8 ik");
91 }
92
93 #define iSH 0
94 #define jSH 1
95 #define kSH 2
96 #define lSH 3
97 #define LOCATE(v, i, j) \
98 int d##i##j = d##i * d##j; \
99 _poutptr = out->outptr + shls[i##SH]*out->v_ket_nsh+shls[j##SH] - out->offset0_outptr; \
100 if (*_poutptr == NOVALUE) { \
101 *_poutptr = out->stack_size; \
102 out->stack_size += d##i##j * ncomp; \
103 NPdset0(out->data+*_poutptr, d##i##j*ncomp); \
104 } \
105 double *v = out->data + *_poutptr;
106 #define DECLARE(v, i, j) \
107 int ncomp = out->ncomp; \
108 int ncol = out->dm_dims[1]; \
109 int di = i1 - i0; \
110 int dj = j1 - j0; \
111 int dk = k1 - k0; \
112 int dl = l1 - l0; \
113 int *_poutptr; \
114 LOCATE(v, i, j)
115
116 #define DEF_NRS1_CONTRACT(D1, D2, V1, V2) \
117 static void nrs1_##D1##D2##_s1##V1##V2(double *eri, double *dm, JKArray *out, int *shls, \
118 int i0, int i1, int j0, int j1, \
119 int k0, int k1, int l0, int l1) \
120 { \
121 DECLARE(v, V1, V2); \
122 int i, j, k, l, ijkl, icomp; \
123 dm += D1##0 * ncol + D2##0 * d##D1; \
124 \
125 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) { \
126 for (l = 0; l < dl; l++) { \
127 for (k = 0; k < dk; k++) { \
128 for (j = 0; j < dj; j++) { \
129 for (i = 0; i < di; i++, ijkl++) { \
130 v[V1*d##V2+V2] += eri[ijkl] * dm[D1*d##D2+D2]; \
131 } } } } \
132 v += d##V1##V2; \
133 } \
134 }
135
136 #define DEF_DM(I, J) \
137 double *dm##I##J = dm + I##0 * ncol + J##0 * d##I
138
139 /* eri in Fortran order; dm, out in C order */
140
nrs1_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)141 static void nrs1_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
142 int i0, int i1, int j0, int j1,
143 int k0, int k1, int l0, int l1)
144 {
145 DECLARE(v, k, l);
146 int dij = di * dj;
147 DEF_DM(j, i);
148 int k, l, ij, icomp;
149 double s;
150
151 for (icomp = 0; icomp < ncomp; icomp++) {
152 for (l = 0; l < dl; l++) {
153 for (k = 0; k < dk; k++) {
154 s = v[k*dl+l];
155 for (ij = 0; ij < dij; ij++) {
156 s += eri[ij] * dmji[ij];
157 }
158 v[k*dl+l] = s;
159 eri += dij;
160 } }
161 v += dkl;
162 }
163 }
164 ADD_JKOP(nrs1_ji_s1kl, J, I, K, L, s1);
165
nrs1_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)166 static void nrs1_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
167 int i0, int i1, int j0, int j1,
168 int k0, int k1, int l0, int l1)
169 {
170 if (k0 >= l0) {
171 nrs1_ji_s1kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
172 }
173 }
174 ADD_JKOP(nrs1_ji_s2kl, J, I, K, L, s1);
175
176
nrs1_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)177 static void nrs1_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
178 int i0, int i1, int j0, int j1,
179 int k0, int k1, int l0, int l1)
180 {
181 DECLARE(v, i, j);
182 DEF_DM(l, k);
183 int i, j, k, l, ij, icomp;
184 double buf[MAXCGTO*MAXCGTO];
185 double s;
186
187 for (icomp = 0; icomp < ncomp; icomp++) {
188
189 for (i = 0; i < dij; i++) { buf[i] = 0; }
190 for (l = 0; l < dl; l++) {
191 for (k = 0; k < dk; k++) {
192 s = dmlk[l*dk+k];
193 for (ij = 0; ij < dij; ij++) {
194 buf[ij] += eri[ij] * s;
195 }
196 eri += dij;
197 } }
198
199 for (ij = 0, j = 0; j < dj; j++) {
200 for (i = 0; i < di; i++, ij++) {
201 v[i*dj+j] += buf[ij];
202 } }
203 v += dij;
204 }
205 }
206 ADD_JKOP(nrs1_lk_s1ij, L, K, I, J, s1);
207
nrs1_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)208 static void nrs1_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
209 int i0, int i1, int j0, int j1,
210 int k0, int k1, int l0, int l1)
211 {
212 if (i0 >= j0) {
213 nrs1_lk_s1ij (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
214 }
215 }
216 ADD_JKOP(nrs1_lk_s2ij, L, K, I, J, s1);
217
218
nrs1_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)219 static void nrs1_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
220 int i0, int i1, int j0, int j1,
221 int k0, int k1, int l0, int l1)
222 {
223 DECLARE(v, i, l);
224 DEF_DM(j, k);
225 int i, j, k, l, ijkl, icomp;
226 double s;
227
228 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
229 for (l = 0; l < dl; l++) {
230 for (k = 0; k < dk; k++) {
231 for (j = 0; j < dj; j++) {
232 s = dmjk[j*dk+k];
233 for (i = 0; i < di; i++, ijkl++) {
234 v[i*dl+l] += eri[ijkl] * s;
235 }
236 } } }
237 v += dil;
238 }
239 }
240 ADD_JKOP(nrs1_jk_s1il, J, K, I, L, s1);
241
242 //DEF_NRS1_CONTRACT(j, k, i, l); ADD_JKOP(nrs1_jk_s1il, J, K, I, L, s1);
243 DEF_NRS1_CONTRACT(j, k, l, i); ADD_JKOP(nrs1_jk_s1li, J, K, L, I, s1);
244 DEF_NRS1_CONTRACT(k, j, i, l); ADD_JKOP(nrs1_kj_s1il, K, J, I, L, s1);
245 DEF_NRS1_CONTRACT(k, j, l, i); ADD_JKOP(nrs1_kj_s1li, K, J, L, I, s1);
246 DEF_NRS1_CONTRACT(i, k, j, l); ADD_JKOP(nrs1_ik_s1jl, I, K, J, L, s1);
247 DEF_NRS1_CONTRACT(i, k, l, j); ADD_JKOP(nrs1_ik_s1lj, I, K, L, J, s1);
248 DEF_NRS1_CONTRACT(k, i, l, j); ADD_JKOP(nrs1_ki_s1lj, K, I, L, J, s1);
249 DEF_NRS1_CONTRACT(k, i, j, l); ADD_JKOP(nrs1_ki_s1jl, K, I, J, L, s1);
250 DEF_NRS1_CONTRACT(j, l, k, i); ADD_JKOP(nrs1_jl_s1ki, J, L, K, I, s1);
251 DEF_NRS1_CONTRACT(j, l, i, k); ADD_JKOP(nrs1_jl_s1ik, J, L, I, K, s1);
252 DEF_NRS1_CONTRACT(l, j, k, i); ADD_JKOP(nrs1_lj_s1ki, L, J, K, I, s1);
253 DEF_NRS1_CONTRACT(l, j, i, k); ADD_JKOP(nrs1_lj_s1ik, L, J, I, K, s1);
254 DEF_NRS1_CONTRACT(l, i, k, j); ADD_JKOP(nrs1_li_s1kj, L, I, K, J, s1);
255 DEF_NRS1_CONTRACT(l, i, j, k); ADD_JKOP(nrs1_li_s1jk, L, I, J, K, s1);
256 DEF_NRS1_CONTRACT(i, l, k, j); ADD_JKOP(nrs1_il_s1kj, I, L, K, J, s1);
257 DEF_NRS1_CONTRACT(i, l, j, k); ADD_JKOP(nrs1_il_s1jk, I, L, J, K, s1);
258
259 //DEF_NRS1_CONTRACT(j, i, k, l); ADD_JKOP(nrs1_ji_s1kl, J, I, K, L, s1);
260 //DEF_NRS1_CONTRACT(l, k, i, j); ADD_JKOP(nrs1_lk_s1ij, L, K, I, J, s1);
261 DEF_NRS1_CONTRACT(i, j, k, l); ADD_JKOP(nrs1_ij_s1kl, I, J, K, L, s1);
262 DEF_NRS1_CONTRACT(i, j, l, k); ADD_JKOP(nrs1_ij_s1lk, I, J, L, K, s1);
263 DEF_NRS1_CONTRACT(j, i, l, k); ADD_JKOP(nrs1_ji_s1lk, J, I, L, K, s1);
264 DEF_NRS1_CONTRACT(l, k, j, i); ADD_JKOP(nrs1_lk_s1ji, L, K, J, I, s1);
265 DEF_NRS1_CONTRACT(k, l, i, j); ADD_JKOP(nrs1_kl_s1ij, K, L, I, J, s1);
266 DEF_NRS1_CONTRACT(k, l, j, i); ADD_JKOP(nrs1_kl_s1ji, K, L, J, I, s1);
267
nrs1_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)268 static void nrs1_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
269 int i0, int i1, int j0, int j1,
270 int k0, int k1, int l0, int l1)
271 {
272 if (i0 >= l0) {
273 nrs1_jk_s1il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
274 }
275 }
276 ADD_JKOP(nrs1_jk_s2il, J, K, I, L, s1);
277
278
nrs1_kj_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)279 static void nrs1_kj_s2il(double *eri, double *dm, JKArray *out, int *shls,
280 int i0, int i1, int j0, int j1,
281 int k0, int k1, int l0, int l1)
282 {
283 if (i0 >= l0) {
284 nrs1_kj_s1il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
285 }
286 }
287 ADD_JKOP(nrs1_kj_s2il, J, K, I, L, s1);
288
289
nrs1_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)290 static void nrs1_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
291 int i0, int i1, int j0, int j1,
292 int k0, int k1, int l0, int l1)
293 {
294 if (k0 >= j0) {
295 nrs1_li_s1kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
296 }
297 }
298 ADD_JKOP(nrs1_li_s2kj, L, I, K, J, s1);
299
nrs2ij_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)300 static void nrs2ij_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
301 int i0, int i1, int j0, int j1,
302 int k0, int k1, int l0, int l1)
303 {
304 if (i0 > j0) {
305 DECLARE(v, k, l);
306 int dij = di * dj;
307 DEF_DM(i, j);
308 DEF_DM(j, i);
309 int i, j, k, l, ij, icomp;
310 double tdm[MAXCGTO*MAXCGTO];
311 double tmp;
312
313 for (ij = 0, j = 0; j < dj; j++) {
314 for (i = 0; i < di; i++, ij++) {
315 tdm[ij] = dmij[i*dj+j] + dmji[j*di+i];
316 } }
317
318 for (icomp = 0; icomp < ncomp; icomp++) {
319 for (l = 0; l < dl; l++) {
320 for (k = 0; k < dk; k++) {
321 tmp = 0;
322 for (ij = 0; ij < dij; ij++) {
323 tmp += eri[ij] * tdm[ij];
324 }
325 v[k*dl+l] += tmp;
326 eri += dij;
327 } }
328 v += dkl;
329 }
330 } else {
331 nrs1_ji_s1kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
332 }
333 }
334 ADD_JKOP(nrs2ij_ji_s1kl, J, I, K, L, s2ij);
335
nrs2ij_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)336 static void nrs2ij_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
337 int i0, int i1, int j0, int j1,
338 int k0, int k1, int l0, int l1)
339 {
340 if (k0 >= l0) {
341 nrs2ij_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
342 }
343 }
344 ADD_JKOP(nrs2ij_ji_s2kl, J, I, K, L, s2ij);
345
346
nrs2ij_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)347 static void nrs2ij_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
348 int i0, int i1, int j0, int j1,
349 int k0, int k1, int l0, int l1)
350 {
351 if (i0 > j0) {
352 DECLARE(vij, i, j);
353 LOCATE(vji, j, i);
354 DEF_DM(l, k);
355 int i, j, k, l, ij, icomp;
356 double buf[MAXCGTO*MAXCGTO];
357 double s;
358
359 for (icomp = 0; icomp < ncomp; icomp++) {
360
361 for (i = 0; i < dij; i++) { buf[i] = 0; }
362 for (l = 0; l < dl; l++) {
363 for (k = 0; k < dk; k++) {
364 s = dmlk[l*dk+k];
365 for (ij = 0; ij < dij; ij++) {
366 buf[ij] += eri[ij] * s;
367 }
368 eri += dij;
369 } }
370
371 for (ij = 0, j = 0; j < dj; j++) {
372 for (i = 0; i < di; i++, ij++) {
373 vij[i*dj+j] += buf[ij];
374 vji[ij] += buf[ij];
375 } }
376 vij += dij;
377 vji += dij;
378 }
379 } else {
380 nrs1_lk_s1ij (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
381 }
382 }
383 ADD_JKOP(nrs2ij_lk_s1ij, L, K, I, J, s2ij);
384
nrs2ij_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)385 static void nrs2ij_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
386 int i0, int i1, int j0, int j1,
387 int k0, int k1, int l0, int l1)
388 {
389 nrs2ij_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
390 }
391 ADD_JKOP(nrs2ij_lk_s2ij, L, K, I, J, s2ij);
392
393
nrs2ij_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)394 static void nrs2ij_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
395 int i0, int i1, int j0, int j1,
396 int k0, int k1, int l0, int l1)
397 {
398 if (i0 > j0) {
399 DECLARE(vil, i, l);
400 DEF_DM(i, k);
401 DEF_DM(j, k);
402 LOCATE(vjl, j, l);
403 int i, j, k, l, ijkl, icomp;
404 double s, tmp;
405
406 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
407 for (l = 0; l < dl; l++) {
408 for (k = 0; k < dk; k++) {
409 for (j = 0; j < dj; j++) {
410 s = dmjk[j*dk+k];
411 tmp = vjl[j*dl+l];
412 for (i = 0; i < di; i++, ijkl++) {
413 vil[i*dl+l] += eri[ijkl] * s;
414 tmp += eri[ijkl] * dmik[i*dk+k];
415 }
416 vjl[j*dl+l] = tmp;
417 } } }
418 vil += dil;
419 vjl += djl;
420 }
421 } else {
422 nrs1_jk_s1il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
423 }
424 }
425 ADD_JKOP(nrs2ij_jk_s1il, J, K, I, L, s2ij);
426
nrs2ij_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)427 static void nrs2ij_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
428 int i0, int i1, int j0, int j1,
429 int k0, int k1, int l0, int l1)
430 {
431 if (j0 >= l0) {
432 nrs2ij_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
433 } else {
434 nrs1_jk_s2il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
435 }
436 }
437 ADD_JKOP(nrs2ij_jk_s2il, J, K, I, L, s2ij);
438
439
nrs2ij_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)440 static void nrs2ij_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
441 int i0, int i1, int j0, int j1,
442 int k0, int k1, int l0, int l1)
443 {
444 if (i0 > j0) {
445 DECLARE(vkj, k, j);
446 DEF_DM(l, i);
447 DEF_DM(l, j);
448 LOCATE(vki, k, i);
449 int i, j, k, l, ijkl, icomp;
450 double s, tmp;
451
452 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
453 for (l = 0; l < dl; l++) {
454 for (k = 0; k < dk; k++) {
455 for (j = 0; j < dj; j++) {
456 s = dmlj[l*dj+j];
457 tmp = vkj[k*dj+j];
458 for (i = 0; i < di; i++, ijkl++) {
459 tmp += eri[ijkl] * dmli[l*di+i];
460 vki[k*di+i] += eri[ijkl] * s;
461 }
462 vkj[k*dj+j] = tmp;
463 } } }
464 vkj += dkj;
465 vki += dki;
466 }
467 } else {
468 nrs1_li_s1kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
469 }
470 }
471 ADD_JKOP(nrs2ij_li_s1kj, L, I, K, J, s2ij);
472
nrs2ij_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)473 static void nrs2ij_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
474 int i0, int i1, int j0, int j1,
475 int k0, int k1, int l0, int l1)
476 {
477 if (k0 >= i0) {
478 nrs2ij_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
479 } else {
480 nrs1_li_s2kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
481 }
482 }
483 ADD_JKOP(nrs2ij_li_s2kj, L, I, K, J, s2ij);
484
485
nrs2kl_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)486 static void nrs2kl_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
487 int i0, int i1, int j0, int j1,
488 int k0, int k1, int l0, int l1)
489 {
490 if (k0 > l0) {
491 DECLARE(vkl, k, l);
492 LOCATE(vlk, l, k);
493 int dij = di * dj;
494 DEF_DM(j, i);
495 int k, l, ij, icomp;
496 double tmp;
497
498 for (icomp = 0; icomp < ncomp; icomp++) {
499 for (l = 0; l < dl; l++) {
500 for (k = 0; k < dk; k++) {
501 tmp = 0;
502 for (ij = 0; ij < dij; ij++) {
503 tmp += eri[ij] * dmji[ij];
504 }
505 vkl[k*dl+l] += tmp;
506 vlk[l*dk+k] += tmp;
507 eri += dij;
508 } }
509 vkl += dkl;
510 vlk += dkl;
511 }
512 } else {
513 return nrs1_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
514 }
515 }
516 ADD_JKOP(nrs2kl_ji_s1kl, J, I, K, L, s2kl);
517
nrs2kl_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)518 static void nrs2kl_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
519 int i0, int i1, int j0, int j1,
520 int k0, int k1, int l0, int l1)
521 {
522 nrs1_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
523 }
524 ADD_JKOP(nrs2kl_ji_s2kl, J, I, K, L, s2kl);
525
526
nrs2kl_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)527 static void nrs2kl_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
528 int i0, int i1, int j0, int j1,
529 int k0, int k1, int l0, int l1)
530 {
531 if (k0 > l0) {
532 DECLARE(v, i, j);
533 DEF_DM(k, l);
534 DEF_DM(l, k);
535 int i, j, k, l, ij, icomp;
536 double buf[MAXCGTO*MAXCGTO];
537 double tdm;
538
539 for (icomp = 0; icomp < ncomp; icomp++) {
540
541 for (i = 0; i < dij; i++) { buf[i] = 0; }
542 for (l = 0; l < dl; l++) {
543 for (k = 0; k < dk; k++) {
544 tdm = dmkl[k*dl+l] + dmlk[l*dk+k];
545 for (ij = 0; ij < dij; ij++) {
546 buf[ij] += eri[ij] * tdm;
547 }
548 eri += dij;
549 } }
550
551 for (ij = 0, j = 0; j < dj; j++) {
552 for (i = 0; i < di; i++, ij++) {
553 v[i*dj+j] += buf[ij];
554 } }
555 v += dij;
556 }
557 } else {
558 nrs1_lk_s1ij (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
559 }
560 }
561 ADD_JKOP(nrs2kl_lk_s1ij, L, K, I, J, s2kl);
562
nrs2kl_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)563 static void nrs2kl_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
564 int i0, int i1, int j0, int j1,
565 int k0, int k1, int l0, int l1)
566 {
567 if (i0 >= j0) {
568 nrs2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
569 }
570 }
571 ADD_JKOP(nrs2kl_lk_s2ij, L, K, I, J, s2kl);
572
573
nrs2kl_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)574 static void nrs2kl_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
575 int i0, int i1, int j0, int j1,
576 int k0, int k1, int l0, int l1)
577 {
578 if (k0 > l0) {
579 DECLARE(vil, i, l);
580 DEF_DM(j, k);
581 DEF_DM(j, l);
582 LOCATE(vik, i, k);
583 int i, j, k, l, ijkl, icomp;
584 double s0, s1;
585
586 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
587 for (l = 0; l < dl; l++) {
588 for (k = 0; k < dk; k++) {
589 for (j = 0; j < dj; j++) {
590 s0 = dmjk[j*dk+k];
591 s1 = dmjl[j*dl+l];
592 for (i = 0; i < di; i++, ijkl++) {
593 vil[i*dl+l] += eri[ijkl] * s0;
594 vik[i*dk+k] += eri[ijkl] * s1;
595 }
596 } } }
597 vil += dil;
598 vik += dik;
599 }
600 } else {
601 nrs1_jk_s1il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
602 }
603 }
604 ADD_JKOP(nrs2kl_jk_s1il, J, K, I, L, s2kl);
605
nrs2kl_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)606 static void nrs2kl_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
607 int i0, int i1, int j0, int j1,
608 int k0, int k1, int l0, int l1)
609 {
610 if (i0 >= k0) {
611 nrs2kl_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
612 } else {
613 nrs1_jk_s2il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
614 }
615 }
616 ADD_JKOP(nrs2kl_jk_s2il, J, K, I, L, s2kl);
617
618
nrs2kl_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)619 static void nrs2kl_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
620 int i0, int i1, int j0, int j1,
621 int k0, int k1, int l0, int l1)
622 {
623 if (k0 > l0) {
624 DECLARE(vkj, k, j);
625 DEF_DM(k, i);
626 DEF_DM(l, i);
627 LOCATE(vlj, l, j);
628 int i, j, k, l, ijkl, icomp;
629 double tmp0, tmp1;
630
631 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
632 for (l = 0; l < dl; l++) {
633 for (k = 0; k < dk; k++) {
634 for (j = 0; j < dj; j++) {
635 tmp0 = vkj[k*dj+j];
636 tmp1 = vlj[l*dj+j];
637 for (i = 0; i < di; i++, ijkl++) {
638 tmp0 += eri[ijkl] * dmli[l*di+i];
639 tmp1 += eri[ijkl] * dmki[k*di+i];
640 }
641 vkj[k*dj+j] = tmp0;
642 vlj[l*dj+j] = tmp1;
643 } } }
644 vkj += dkj;
645 vlj += dlj;
646 }
647 } else {
648 nrs1_li_s1kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
649 }
650 }
651 ADD_JKOP(nrs2kl_li_s1kj, L, I, K, J, s2kl);
652
nrs2kl_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)653 static void nrs2kl_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
654 int i0, int i1, int j0, int j1,
655 int k0, int k1, int l0, int l1)
656 {
657 if (l0 >= j0) {
658 nrs2kl_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
659 } else {
660 nrs1_li_s2kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
661 }
662 }
663 ADD_JKOP(nrs2kl_li_s2kj, L, I, K, J, s2kl);
664
665
nrs4_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)666 static void nrs4_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
667 int i0, int i1, int j0, int j1,
668 int k0, int k1, int l0, int l1)
669 {
670 if (i0 == j0) {
671 nrs2kl_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
672 } else if (k0 == l0) {
673 nrs2ij_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
674 } else {
675 DECLARE(vkl, k, l);
676 LOCATE(vlk, l, k);
677 int dij = di * dj;
678 DEF_DM(i, j);
679 DEF_DM(j, i);
680 int i, j, k, l, ij, icomp;
681 double tdm[MAXCGTO*MAXCGTO];
682 double tmp;
683
684 for (ij = 0, j = 0; j < dj; j++) {
685 for (i = 0; i < di; i++, ij++) {
686 tdm[ij] = dmij[i*dj+j] + dmji[j*di+i];
687 }
688 }
689
690 for (icomp = 0; icomp < ncomp; icomp++) {
691 for (l = 0; l < dl; l++) {
692 for (k = 0; k < dk; k++) {
693 tmp = 0;
694 for (ij = 0; ij < dij; ij++) {
695 tmp += eri[ij] * tdm[ij];
696 }
697 vkl[k*dl+l] += tmp;
698 vlk[l*dk+k] += tmp;
699 eri += dij;
700 } }
701 vkl += dkl;
702 vlk += dkl;
703 }
704 }
705 }
706 ADD_JKOP(nrs4_ji_s1kl, J, I, K, L, s4);
707
nrs4_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)708 static void nrs4_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
709 int i0, int i1, int j0, int j1,
710 int k0, int k1, int l0, int l1)
711 {
712 nrs2ij_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
713 }
714 ADD_JKOP(nrs4_ji_s2kl, J, I, K, L, s4);
715
716
nrs4_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)717 static void nrs4_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
718 int i0, int i1, int j0, int j1,
719 int k0, int k1, int l0, int l1)
720 {
721 if (i0 == j0) {
722 nrs2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
723 } else if (k0 == l0) {
724 nrs2ij_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
725 } else {
726 DECLARE(vij, i, j);
727 LOCATE(vji, j, i);
728 DEF_DM(k, l);
729 DEF_DM(l, k);
730 int i, j, k, l, ij, icomp;
731 double buf[MAXCGTO*MAXCGTO];
732 double tdm;
733
734 for (icomp = 0; icomp < ncomp; icomp++) {
735
736 for (i = 0; i < dij; i++) { buf[i] = 0; }
737 for (l = 0; l < dl; l++) {
738 for (k = 0; k < dk; k++) {
739 tdm = dmlk[l*dk+k] + dmkl[k*dl+l];
740 for (ij = 0; ij < dij; ij++) {
741 buf[ij] += eri[ij] * tdm;
742 }
743 eri += dij;
744 } }
745
746 for (ij = 0, j = 0; j < dj; j++) {
747 for (i = 0; i < di; i++, ij++) {
748 vij[i*dj+j] += buf[ij];
749 vji[ij] += buf[ij];
750 } }
751 vij += dij;
752 vji += dij;
753 }
754 }
755 }
756 ADD_JKOP(nrs4_lk_s1ij, L, K, I, J, s4);
757
nrs4_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)758 static void nrs4_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
759 int i0, int i1, int j0, int j1,
760 int k0, int k1, int l0, int l1)
761 {
762 nrs2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
763 }
764 ADD_JKOP(nrs4_lk_s2ij, L, K, I, J, s4);
765
nrs4_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)766 static void nrs4_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
767 int i0, int i1, int j0, int j1,
768 int k0, int k1, int l0, int l1)
769 {
770 if (i0 == j0) {
771 nrs2kl_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
772 } else if (k0 == l0) {
773 nrs2ij_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
774 } else {
775 DECLARE(vik, i, k);
776 LOCATE(vil, i, l);
777 LOCATE(vjk, j, k);
778 LOCATE(vjl, j, l);
779 DEF_DM(i, l);
780 DEF_DM(i, k);
781 DEF_DM(j, l);
782 DEF_DM(j, k);
783 int i, j, k, l, ijkl, icomp;
784 double s0, s1, tmp0, tmp1;
785
786 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
787 for (l = 0; l < dl; l++) {
788 for (k = 0; k < dk; k++) {
789 for (j = 0; j < dj; j++) {
790 s0 = dmjl[j*dl+l];
791 s1 = dmjk[j*dk+k];
792 tmp0 = vjk[j*dk+k];
793 tmp1 = vjl[j*dl+l];
794 for (i = 0; i < di; i++, ijkl++) {
795 tmp0 += eri[ijkl] * dmil[i*dl+l];
796 tmp1 += eri[ijkl] * dmik[i*dk+k];
797 vik[i*dk+k] += eri[ijkl] * s0;
798 vil[i*dl+l] += eri[ijkl] * s1;
799 }
800 vjk[j*dk+k] = tmp0;
801 vjl[j*dl+l] = tmp1;
802 } } }
803 vjk += djk;
804 vjl += djl;
805 vik += dik;
806 vil += dil;
807 }
808 }
809 }
810 ADD_JKOP(nrs4_jk_s1il, J, K, I, L, s4);
811
nrs4_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)812 static void nrs4_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
813 int i0, int i1, int j0, int j1,
814 int k0, int k1, int l0, int l1)
815 {
816 if (i0 == j0) {
817 nrs2kl_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
818 } else if (k0 == l0) {
819 nrs2ij_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
820 } else if (i0 < l0) {
821 } else if (i0 < k0) {
822 if (j0 < l0) { // j < l <= i < k
823 DECLARE(v, i, l);
824 DEF_DM(j, k);
825 int i, j, k, l, ijkl, icomp;
826 double s;
827 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
828 for (l = 0; l < dl; l++) {
829 for (k = 0; k < dk; k++) {
830 for (j = 0; j < dj; j++) {
831 s = dmjk[j*dk+k];
832 for (i = 0; i < di; i++, ijkl++) {
833 v[i*dl+l] += eri[ijkl] * s;
834 }
835 } } }
836 v += dil;
837 }
838 } else { // l <= j < i < k
839 DECLARE(vil, i, l);
840 LOCATE(vjl, j, l);
841 DEF_DM(i, k);
842 DEF_DM(j, k);
843 int i, j, k, l, ijkl, icomp;
844 double s, tmp;
845 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
846 for (l = 0; l < dl; l++) {
847 for (k = 0; k < dk; k++) {
848 for (j = 0; j < dj; j++) {
849 s = dmjk[j*dk+k];
850 tmp = vjl[j*dl+l];
851 for (i = 0; i < di; i++, ijkl++) {
852 tmp += eri[ijkl] * dmik[i*dk+k];
853 vil[i*dl+l] += eri[ijkl] * s;
854 }
855 vjl[j*dl+l] = tmp;
856 } } }
857 vjl += djl;
858 vil += dil;
859 }
860 }
861 } else if (j0 < l0) { // j < l < k <= i
862 DECLARE(vil, i, l);
863 LOCATE(vik, i, k);
864 DEF_DM(j, k);
865 DEF_DM(j, l);
866 int i, j, k, l, ijkl, icomp;
867 double s0, s1;
868 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
869 for (l = 0; l < dl; l++) {
870 for (k = 0; k < dk; k++) {
871 for (j = 0; j < dj; j++) {
872 s0 = dmjk[j*dk+k];
873 s1 = dmjl[j*dl+l];
874 for (i = 0; i < di; i++, ijkl++) {
875 vil[i*dl+l] += eri[ijkl] * s0;
876 vik[i*dk+k] += eri[ijkl] * s1;
877 }
878 } } }
879 vil += dil;
880 vik += dik;
881 }
882 } else if (j0 < k0) { // l <= j < k <= i
883 DECLARE(vjl, j, l);
884 LOCATE(vik, i, k);
885 LOCATE(vil, i, l);
886 DEF_DM(i, k);
887 DEF_DM(j, k);
888 DEF_DM(j, l);
889 int i, j, k, l, ijkl, icomp;
890 double s0, s1, tmp;
891 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
892 for (l = 0; l < dl; l++) {
893 for (k = 0; k < dk; k++) {
894 for (j = 0; j < dj; j++) {
895 s0 = dmjk[j*dk+k];
896 s1 = dmjl[j*dl+l];
897 tmp = vjl[j*dl+l];
898 for (i = 0; i < di; i++, ijkl++) {
899 tmp += eri[ijkl] * dmik[i*dk+k];
900 vil[i*dl+l] += eri[ijkl] * s0;
901 vik[i*dk+k] += eri[ijkl] * s1;
902 }
903 vjl[j*dl+l] = tmp;
904 } } }
905 vjl += djl;
906 vil += dil;
907 vik += dik;
908 }
909 } else { // l < k <= j < i
910 DECLARE(vjl, j, l);
911 LOCATE(vik, i, k);
912 LOCATE(vjk, j, k);
913 LOCATE(vil, i, l);
914 DEF_DM(i, l);
915 DEF_DM(i, k);
916 DEF_DM(j, l);
917 DEF_DM(j, k);
918 int i, j, k, l, ijkl, icomp;
919 double s0, s1, tmp0, tmp1;
920 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
921 for (l = 0; l < dl; l++) {
922 for (k = 0; k < dk; k++) {
923 for (j = 0; j < dj; j++) {
924 s0 = dmjl[j*dl+l];
925 s1 = dmjk[j*dk+k];
926 tmp0 = vjk[j*dk+k];
927 tmp1 = vjl[j*dl+l];
928 for (i = 0; i < di; i++, ijkl++) {
929 tmp0 += eri[ijkl] * dmil[i*dl+l];
930 tmp1 += eri[ijkl] * dmik[i*dk+k];
931 vik[i*dk+k] += eri[ijkl] * s0;
932 vil[i*dl+l] += eri[ijkl] * s1;
933 }
934 vjk[j*dk+k] = tmp0;
935 vjl[j*dl+l] = tmp1;
936 } } }
937 vjk += djk;
938 vjl += djl;
939 vik += dik;
940 vil += dil;
941 }
942 }
943 }
944 ADD_JKOP(nrs4_jk_s2il, J, K, I, L, s4);
945
946
nrs4_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)947 static void nrs4_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
948 int i0, int i1, int j0, int j1,
949 int k0, int k1, int l0, int l1)
950 {
951 if (i0 == j0) {
952 nrs2kl_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
953 } else if (k0 == l0) {
954 nrs2ij_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
955 } else {
956 DECLARE(vki, k, i);
957 LOCATE(vkj, k, j);
958 LOCATE(vli, l, i);
959 LOCATE(vlj, l, j);
960 DEF_DM(l, i);
961 DEF_DM(l, j);
962 DEF_DM(k, i);
963 DEF_DM(k, j);
964 int i, j, k, l, ijkl, icomp;
965 double s0, s1, tmp0, tmp1;
966 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
967 for (l = 0; l < dl; l++) {
968 for (k = 0; k < dk; k++) {
969 for (j = 0; j < dj; j++) {
970 s0 = dmlj[l*dj+j];
971 s1 = dmkj[k*dj+j];
972 tmp0 = vkj[k*dj+j];
973 tmp1 = vlj[l*dj+j];
974 for (i = 0; i < di; i++, ijkl++) {
975 tmp0 += eri[ijkl] * dmli[l*di+i];
976 tmp1 += eri[ijkl] * dmki[k*di+i];
977 vki[k*di+i] += eri[ijkl] * s0;
978 vli[l*di+i] += eri[ijkl] * s1;
979 }
980 vkj[k*dj+j] = tmp0;
981 vlj[l*dj+j] = tmp1;
982 } } }
983 vkj += dkj;
984 vki += dki;
985 vlj += dlj;
986 vli += dli;
987 }
988 }
989 }
990 ADD_JKOP(nrs4_li_s1kj, L, I, K, J, s4);
991
nrs4_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)992 static void nrs4_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
993 int i0, int i1, int j0, int j1,
994 int k0, int k1, int l0, int l1)
995 {
996 if (i0 == j0) {
997 nrs2kl_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
998 } else if (k0 == l0) {
999 nrs2ij_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1000 } else if (k0 < j0) {
1001 } else if (k0 < i0) {
1002 if (l0 < j0) { // l < j < k < i
1003 DECLARE(v, k, j);
1004 DEF_DM(l, i);
1005 int i, j, k, l, ijkl, icomp;
1006 double tmp;
1007 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1008 for (l = 0; l < dl; l++) {
1009 for (k = 0; k < dk; k++) {
1010 for (j = 0; j < dj; j++) {
1011 tmp = v[k*dj+j];
1012 for (i = 0; i < di; i++, ijkl++) {
1013 tmp += eri[ijkl] * dmli[l*di+i];
1014 }
1015 v[k*dj+j] = tmp;
1016 } } }
1017 v += dkj;
1018 }
1019 } else { // j <= l < k < i
1020 DECLARE(vkj, k, j);
1021 LOCATE(vlj, l, j);
1022 DEF_DM(l, i);
1023 DEF_DM(k, i);
1024 int i, j, k, l, ijkl, icomp;
1025 double tmp0, tmp1;
1026 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1027 for (l = 0; l < dl; l++) {
1028 for (k = 0; k < dk; k++) {
1029 for (j = 0; j < dj; j++) {
1030 tmp0 = vkj[k*dj+j];
1031 tmp1 = vlj[l*dj+j];
1032 for (i = 0; i < di; i++, ijkl++) {
1033 tmp0 += eri[ijkl] * dmli[l*di+i];
1034 tmp1 += eri[ijkl] * dmki[k*di+i];
1035 }
1036 vkj[k*dj+j] = tmp0;
1037 vlj[l*dj+j] = tmp1;
1038 } } }
1039 vkj += dkj;
1040 vlj += dlj;
1041 }
1042 }
1043 } else if (l0 < j0) { // l < j < i <= k
1044 DECLARE(vki, k, i);
1045 LOCATE(vkj, k, j);
1046 DEF_DM(l, i);
1047 DEF_DM(l, j);
1048 int i, j, k, l, ijkl, icomp;
1049 double s, tmp;
1050 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1051 for (l = 0; l < dl; l++) {
1052 for (k = 0; k < dk; k++) {
1053 for (j = 0; j < dj; j++) {
1054 s = dmlj[l*dj+j];
1055 tmp = vkj[k*dj+j];
1056 for (i = 0; i < di; i++, ijkl++) {
1057 tmp += eri[ijkl] * dmli[l*di+i];
1058 vki[k*di+i] += eri[ijkl] * s;
1059 }
1060 vkj[k*dj+j] = tmp;
1061 } } }
1062 vkj += dkj;
1063 vki += dki;
1064 }
1065 } else if (l0 < i0) { // j <= l < i <= k
1066 DECLARE(vki, k, i);
1067 LOCATE(vkj, k, j);
1068 LOCATE(vlj, l, j);
1069 DEF_DM(l, i);
1070 DEF_DM(l, j);
1071 DEF_DM(k, i);
1072 int i, j, k, l, ijkl, icomp;
1073 double s, tmp0, tmp1;
1074 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1075 for (l = 0; l < dl; l++) {
1076 for (k = 0; k < dk; k++) {
1077 for (j = 0; j < dj; j++) {
1078 s = dmlj[l*dj+j];
1079 tmp0 = vkj[k*dj+j];
1080 tmp1 = vlj[l*dj+j];
1081 for (i = 0; i < di; i++, ijkl++) {
1082 vki[k*di+i] += eri[ijkl] * s;
1083 tmp0 += eri[ijkl] * dmli[l*di+i];
1084 tmp1 += eri[ijkl] * dmki[k*di+i];
1085 }
1086 vkj[k*dj+j] = tmp0;
1087 vlj[l*dj+j] = tmp1;
1088 } } }
1089 vkj += dkj;
1090 vki += dki;
1091 vlj += dlj;
1092 }
1093 } else { // j < i <= l < k
1094 DECLARE(vki, k, i);
1095 LOCATE(vkj, k, j);
1096 LOCATE(vli, l, i);
1097 LOCATE(vlj, l, j);
1098 DEF_DM(l, i);
1099 DEF_DM(l, j);
1100 DEF_DM(k, i);
1101 DEF_DM(k, j);
1102 int i, j, k, l, ijkl, icomp;
1103 double s0, s1, tmp0, tmp1;
1104 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1105 for (l = 0; l < dl; l++) {
1106 for (k = 0; k < dk; k++) {
1107 for (j = 0; j < dj; j++) {
1108 s0 = dmlj[l*dj+j];
1109 s1 = dmkj[k*dj+j];
1110 tmp0 = vkj[k*dj+j];
1111 tmp1 = vlj[l*dj+j];
1112 for (i = 0; i < di; i++, ijkl++) {
1113 tmp0 += eri[ijkl] * dmli[l*di+i];
1114 tmp1 += eri[ijkl] * dmki[k*di+i];
1115 vki[k*di+i] += eri[ijkl] * s0;
1116 vli[l*di+i] += eri[ijkl] * s1;
1117 }
1118 vkj[k*dj+j] = tmp0;
1119 vlj[l*dj+j] = tmp1;
1120 } } }
1121 vkj += dkj;
1122 vki += dki;
1123 vlj += dlj;
1124 vli += dli;
1125 }
1126 }
1127 }
1128 ADD_JKOP(nrs4_li_s2kj, L, I, K, J, s4);
1129
1130
nrs8_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1131 static void nrs8_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
1132 int i0, int i1, int j0, int j1,
1133 int k0, int k1, int l0, int l1)
1134 {
1135 if (i0 == k0 && j0 == l0) {
1136 nrs4_ji_s1kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1137 } else if (i0 == j0 || k0 == l0) {
1138 nrs4_ji_s1kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1139 nrs4_lk_s1ij (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1140 } else {
1141 DECLARE(vij, i, j);
1142 LOCATE(vji, j, i);
1143 LOCATE(vkl, k, l);
1144 LOCATE(vlk, l, k);
1145 DEF_DM(i, j);
1146 DEF_DM(j, i);
1147 DEF_DM(k, l);
1148 DEF_DM(l, k);
1149 int i, j, k, l, ij, icomp;
1150 double tdm[MAXCGTO*MAXCGTO];
1151 double buf[MAXCGTO*MAXCGTO];
1152 double tdm2, tmp;
1153
1154 for (icomp = 0; icomp < ncomp; icomp++) {
1155 for (ij = 0, j = 0; j < dj; j++) {
1156 for (i = 0; i < di; i++, ij++) {
1157 tdm[ij] = dmij[i*dj+j] + dmji[j*di+i];
1158 buf[ij] = 0;
1159 }
1160 }
1161
1162 for (l = 0; l < dl; l++) {
1163 for (k = 0; k < dk; k++) {
1164 tmp = 0;
1165 tdm2 = dmkl[k*dl+l] + dmlk[l*dk+k];
1166 for (ij = 0; ij < dij; ij++) {
1167 tmp += eri[ij] * tdm[ij];
1168 buf[ij] += eri[ij] * tdm2;
1169 }
1170 vkl[k*dl+l] += tmp;
1171 vlk[l*dk+k] += tmp;
1172 eri += dij;
1173 } }
1174
1175 for (ij = 0, j = 0; j < dj; j++) {
1176 for (i = 0; i < di; i++, ij++) {
1177 vij[i*dj+j] += buf[ij];
1178 vji[ij] += buf[ij];
1179 } }
1180 vij += dij;
1181 vji += dji;
1182 vkl += dkl;
1183 vlk += dlk;
1184 }
1185 }
1186 }
1187 ADD_JKOP(nrs8_ji_s1kl, J, I, K, L, s8);
1188
nrs8_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1189 static void nrs8_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
1190 int i0, int i1, int j0, int j1,
1191 int k0, int k1, int l0, int l1)
1192 {
1193 if (i0 == k0 && j0 == l0) {
1194 nrs4_ji_s2kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1195 } else if (i0 == j0 || k0 == l0) {
1196 nrs4_ji_s2kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1197 nrs4_lk_s2ij (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1198 } else {
1199 DECLARE(vij, i, j);
1200 LOCATE(vkl, k, l);
1201 DEF_DM(i, j);
1202 DEF_DM(j, i);
1203 DEF_DM(k, l);
1204 DEF_DM(l, k);
1205 int i, j, k, l, ij, icomp;
1206 double tdm[MAXCGTO*MAXCGTO];
1207 double buf[MAXCGTO*MAXCGTO];
1208 double tmp, tdm2;
1209
1210 for (icomp = 0; icomp < ncomp; icomp++) {
1211 for (ij = 0, j = 0; j < dj; j++) {
1212 for (i = 0; i < di; i++, ij++) {
1213 tdm[ij] = dmij[i*dj+j] + dmji[j*di+i];
1214 buf[ij] = 0;
1215 }
1216 }
1217
1218 for (l = 0; l < dl; l++) {
1219 for (k = 0; k < dk; k++) {
1220 tdm2 = dmkl[k*dl+l] + dmlk[l*dk+k];
1221 tmp = 0;
1222 for (ij = 0; ij < dij; ij++) {
1223 buf[ij] += eri[ij] * tdm2;
1224 tmp += eri[ij] * tdm[ij];
1225 }
1226 vkl[k*dl+l] += tmp;
1227 eri += dij;
1228 } }
1229
1230 for (ij = 0, i = 0; i < di; i++) {
1231 for (j = 0; j < dj; j++, ij++) {
1232 vij[ij] += buf[j*di+i];
1233 } }
1234 vij += dij;
1235 vkl += dkl;
1236 }
1237 }
1238 }
1239 ADD_JKOP(nrs8_ji_s2kl, J, I, K, L, s8);
1240
1241
nrs8_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1242 static void nrs8_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
1243 int i0, int i1, int j0, int j1,
1244 int k0, int k1, int l0, int l1)
1245 {
1246 nrs8_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1247 }
1248 ADD_JKOP(nrs8_lk_s1ij, L, K, I, J, s8);
1249
nrs8_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1250 static void nrs8_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
1251 int i0, int i1, int j0, int j1,
1252 int k0, int k1, int l0, int l1)
1253 {
1254 nrs8_ji_s2kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1255 }
1256 ADD_JKOP(nrs8_lk_s2ij, L, K, I, J, s8);
1257
1258
nrs8_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1259 static void nrs8_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
1260 int i0, int i1, int j0, int j1,
1261 int k0, int k1, int l0, int l1)
1262 {
1263 if (i0 == k0 && j0 == l0) {
1264 nrs4_li_s1kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1265 } else if (i0 == j0 || k0 == l0) { // i0==l0 => i0==k0==l0
1266 nrs4_li_s1kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1267 nrs4_jk_s1il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1268 } else {
1269 DECLARE(vkj, k, j);
1270 LOCATE(vki, k, i);
1271 LOCATE(vlj, l, j);
1272 LOCATE(vli, l, i);
1273 LOCATE(vik, i, k);
1274 LOCATE(vil, i, l);
1275 LOCATE(vjk, j, k);
1276 LOCATE(vjl, j, l);
1277 DEF_DM(l, i);
1278 DEF_DM(l, j);
1279 DEF_DM(k, i);
1280 DEF_DM(k, j);
1281 DEF_DM(j, l);
1282 DEF_DM(j, k);
1283 DEF_DM(i, l);
1284 DEF_DM(i, k);
1285 int i, j, k, l, ijkl, icomp;
1286 double s, s0, s1, s2, s3, tmp0, tmp1, tmp2, tmp3;
1287 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1288 for (l = 0; l < dl; l++) {
1289 for (k = 0; k < dk; k++) {
1290 for (j = 0; j < dj; j++) {
1291 s0 = dmlj[l*dj+j];
1292 s1 = dmkj[k*dj+j];
1293 s2 = dmjl[j*dl+l];
1294 s3 = dmjk[j*dk+k];
1295 tmp0 = 0;
1296 tmp1 = 0;
1297 tmp2 = 0;
1298 tmp3 = 0;
1299 for (i = 0; i < di; i++, ijkl++) {
1300 s = eri[ijkl];
1301 vki[k*di+i] += s * s0;
1302 vli[l*di+i] += s * s1;
1303 vik[i*dk+k] += s * s2;
1304 vil[i*dl+l] += s * s3;
1305 tmp0 += s * dmli[l*di+i];
1306 tmp1 += s * dmki[k*di+i];
1307 tmp2 += s * dmil[i*dl+l];
1308 tmp3 += s * dmik[i*dk+k];
1309 }
1310 vkj[k*dj+j] += tmp0;
1311 vlj[l*dj+j] += tmp1;
1312 vjk[j*dk+k] += tmp2; // vkj, vjk may share memory
1313 vjl[j*dl+l] += tmp3; // vlj, vjl may share memory
1314 } } }
1315 vkj += dkj;
1316 vki += dki;
1317 vlj += dlj;
1318 vli += dli;
1319 vik += dik;
1320 vil += dil;
1321 vjk += djk;
1322 vjl += djl;
1323 }
1324 }
1325 }
1326 ADD_JKOP(nrs8_li_s1kj, L, I, K, J, s8);
1327
nrs8_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1328 static void nrs8_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
1329 int i0, int i1, int j0, int j1,
1330 int k0, int k1, int l0, int l1)
1331 {
1332 if (i0 == k0) {
1333 nrs4_li_s2kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1334 if (j0 != l0) {
1335 nrs4_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1336 }
1337 } else if (i0 == j0 || k0 == l0) { // i0==l0 => i0==k0==l0
1338 nrs4_li_s2kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1339 nrs4_jk_s2il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1340 } else {
1341 int i, j, k, l, ijkl, icomp;
1342 double s, tjl, tjk, sjk, sjl, skj, slj;
1343 if (j0 < l0) { // j <= l < k < i
1344 DECLARE(vkj, k, j);
1345 LOCATE(vlj, l, j);
1346 LOCATE(vik, i, k);
1347 LOCATE(vil, i, l);
1348 DEF_DM(l, i);
1349 DEF_DM(k, i);
1350 DEF_DM(j, l);
1351 DEF_DM(j, k);
1352 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1353 for (l = 0; l < dl; l++) {
1354 for (k = 0; k < dk; k++) {
1355 for (j = 0; j < dj; j++) {
1356 tjl = dmjl[j*dl+l];
1357 tjk = dmjk[j*dk+k];
1358 skj = 0;
1359 slj = 0;
1360 for (i = 0; i < di; i++, ijkl++) {
1361 //vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
1362 //vlj[l*dj+j] += eri[ijkl] * dmki[k*di+i];
1363 //vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
1364 //vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
1365 s = eri[ijkl];
1366 skj += s * dmli[l*di+i];
1367 slj += s * dmki[k*di+i];
1368 vik[i*dk+k] += s * tjl;
1369 vil[i*dl+l] += s * tjk;
1370 }
1371 vkj[k*dj+j] += skj;
1372 vlj[l*dj+j] += slj;
1373 } } }
1374 vkj += dkj;
1375 vlj += dlj;
1376 vik += dik;
1377 vil += dil;
1378 }
1379 } else if (j0 == l0) { // j == l < k < i
1380 DECLARE(vkj, k, j);
1381 LOCATE(vlj, l, j);
1382 LOCATE(vik, i, k);
1383 LOCATE(vil, i, l);
1384 LOCATE(vjl, j, l);
1385 DEF_DM(l, i);
1386 DEF_DM(k, i);
1387 DEF_DM(j, l);
1388 DEF_DM(j, k);
1389 DEF_DM(i, k);
1390 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1391 for (l = 0; l < dl; l++) {
1392 for (k = 0; k < dk; k++) {
1393 for (j = 0; j < dj; j++) {
1394 tjl = dmjl[j*dl+l];
1395 tjk = dmjk[j*dk+k];
1396 skj = 0;
1397 slj = 0;
1398 sjl = 0;
1399 for (i = 0; i < di; i++, ijkl++) {
1400 //vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
1401 //vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
1402 //vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
1403 //vlj[l*dj+j] += eri[ijkl] * dmki[k*di+i];
1404 //vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
1405 s = eri[ijkl];
1406 vik[i*dk+k] += s * tjl;
1407 vil[i*dl+l] += s * tjk;
1408 skj += s * dmli[l*di+i];
1409 slj += s * dmki[k*di+i];
1410 sjl += s * dmik[i*dk+k];
1411 }
1412 vlj[l*dj+j] += slj;
1413 vkj[k*dj+j] += skj;
1414 vjl[j*dl+l] += sjl; // vjl, vlj may share memory
1415 } } }
1416 vkj += dkj;
1417 vlj += dlj;
1418 vik += dik;
1419 vil += dil;
1420 vjl += djl;
1421 }
1422 } else if (j0 < k0) { // l < j < k < i
1423 DECLARE(vkj, k, j);
1424 LOCATE(vik, i, k);
1425 LOCATE(vil, i, l);
1426 LOCATE(vjl, j, l);
1427 DEF_DM(l, i);
1428 DEF_DM(j, l);
1429 DEF_DM(j, k);
1430 DEF_DM(i, k);
1431 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1432 for (l = 0; l < dl; l++) {
1433 for (k = 0; k < dk; k++) {
1434 for (j = 0; j < dj; j++) {
1435 tjl = dmjl[j*dl+l];
1436 tjk = dmjk[j*dk+k];
1437 skj = 0;
1438 sjl = 0;
1439 for (i = 0; i < di; i++, ijkl++) {
1440 //vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
1441 //vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
1442 //vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
1443 //vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
1444 s = eri[ijkl];
1445 vik[i*dk+k] += s * tjl;
1446 vil[i*dl+l] += s * tjk;
1447 skj += s * dmli[l*di+i];
1448 sjl += s * dmik[i*dk+k];
1449 }
1450 vkj[k*dj+j] += skj;
1451 vjl[j*dl+l] += sjl;
1452 } } }
1453 vkj += dkj;
1454 vik += dik;
1455 vil += dil;
1456 vjl += djl;
1457 }
1458 } else if (j0 == k0) { // l < j == k < i
1459 DECLARE(vkj, k, j);
1460 LOCATE(vik, i, k);
1461 LOCATE(vil, i, l);
1462 LOCATE(vjk, j, k);
1463 LOCATE(vjl, j, l);
1464 DEF_DM(l, i);
1465 DEF_DM(i, l);
1466 DEF_DM(j, l);
1467 DEF_DM(j, k);
1468 DEF_DM(i, k);
1469 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1470 for (l = 0; l < dl; l++) {
1471 for (k = 0; k < dk; k++) {
1472 for (j = 0; j < dj; j++) {
1473 tjl = dmjl[j*dl+l];
1474 tjk = dmjk[j*dk+k];
1475 sjk = 0;
1476 sjl = 0;
1477 skj = 0;
1478 for (i = 0; i < di; i++, ijkl++) {
1479 //vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
1480 //vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
1481 //vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
1482 //vjk[j*dk+k] += eri[ijkl] * dmil[i*dl+l];
1483 //vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
1484 s = eri[ijkl];
1485 vik[i*dk+k] += s * tjl;
1486 vil[i*dl+l] += s * tjk;
1487 skj += s * dmli[l*di+i];
1488 sjk += s * dmil[i*dl+l];
1489 sjl += s * dmik[i*dk+k];
1490 }
1491 vjk[j*dk+k] += sjk;
1492 vjl[j*dl+l] += sjl;
1493 vkj[k*dj+j] += skj; // vjk, vkj may share memory
1494 } } }
1495 vkj += dkj;
1496 vjk += djk;
1497 vik += dik;
1498 vil += dil;
1499 vjl += djl;
1500 }
1501 } else { // l < k < j < i
1502 DECLARE(vik, i, k);
1503 LOCATE(vil, i, l);
1504 LOCATE(vjk, j, k);
1505 LOCATE(vjl, j, l);
1506 DEF_DM(j, l);
1507 DEF_DM(j, k);
1508 DEF_DM(i, l);
1509 DEF_DM(i, k);
1510 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1511 for (l = 0; l < dl; l++) {
1512 for (k = 0; k < dk; k++) {
1513 for (j = 0; j < dj; j++) {
1514 tjl = dmjl[j*dl+l];
1515 tjk = dmjk[j*dk+k];
1516 sjk = 0;
1517 sjl = 0;
1518 for (i = 0; i < di; i++, ijkl++) {
1519 //vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
1520 //vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
1521 //vjk[j*dk+k] += eri[ijkl] * dmil[i*dl+l];
1522 //vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
1523 s = eri[ijkl];
1524 vik[i*dk+k] += s * tjl;
1525 vil[i*dl+l] += s * tjk;
1526 sjk += s * dmil[i*dl+l];
1527 sjl += s * dmik[i*dk+k];
1528 }
1529 vjk[j*dk+k] += sjk;
1530 vjl[j*dl+l] += sjl;
1531 } } }
1532 vik += dik;
1533 vil += dil;
1534 vjk += djk;
1535 vjl += djl;
1536 }
1537 }
1538 }
1539 }
1540 ADD_JKOP(nrs8_li_s2kj, L, I, K, J, s8);
1541
1542
nrs8_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1543 static void nrs8_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
1544 int i0, int i1, int j0, int j1,
1545 int k0, int k1, int l0, int l1)
1546 {
1547 nrs8_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1548 }
1549 ADD_JKOP(nrs8_jk_s1il, J, K, I, L, s8);
1550
nrs8_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1551 static void nrs8_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
1552 int i0, int i1, int j0, int j1,
1553 int k0, int k1, int l0, int l1)
1554 {
1555 nrs8_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1556 }
1557 ADD_JKOP(nrs8_jk_s2il, J, K, I, L, s8);
1558
1559
1560 /*************************************************
1561 * For anti symmetrized integrals
1562 *************************************************/
nra2ij_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1563 static void nra2ij_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
1564 int i0, int i1, int j0, int j1,
1565 int k0, int k1, int l0, int l1)
1566 {
1567 if (i0 > j0) {
1568 DECLARE(v, k, l);
1569 DEF_DM(i, j);
1570 DEF_DM(j, i);
1571 int dij = di * dj;
1572 int i, j, k, l, ij, icomp;
1573 double tdm[MAXCGTO*MAXCGTO];
1574 double tmp;
1575
1576 for (ij = 0, j = 0; j < dj; j++) {
1577 for (i = 0; i < di; i++, ij++) {
1578 tdm[ij] = dmji[j*di+i] - dmij[i*dj+j];
1579 }
1580 }
1581
1582 for (icomp = 0; icomp < ncomp; icomp++) {
1583 for (l = 0; l < dl; l++) {
1584 for (k = 0; k < dk; k++) {
1585 tmp = 0;
1586 for (ij = 0; ij < dij; ij++) {
1587 tmp += eri[ij] * tdm[ij];
1588 }
1589 v[k*dl+l] += tmp;
1590 eri += dij;
1591 } }
1592 v += dkl;
1593 }
1594 } else {
1595 nrs1_ji_s1kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1596 }
1597 }
1598 ADD_JKOP(nra2ij_ji_s1kl, J, I, K, L, s2ij);
1599
nra2ij_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1600 static void nra2ij_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
1601 int i0, int i1, int j0, int j1,
1602 int k0, int k1, int l0, int l1)
1603 {
1604 if (k0 >= l0) {
1605 nra2ij_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1606 }
1607 }
1608 ADD_JKOP(nra2ij_ji_s2kl, J, I, K, L, s2ij);
1609
nra2ij_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1610 static void nra2ij_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
1611 int i0, int i1, int j0, int j1,
1612 int k0, int k1, int l0, int l1)
1613 {
1614 if (i0 > j0) {
1615 DECLARE(vij, i, j);
1616 LOCATE(vji, j, i);
1617 DEF_DM(l, k);
1618 int i, j, k, l, ij, icomp;
1619 double buf[MAXCGTO*MAXCGTO];
1620
1621 for (icomp = 0; icomp < ncomp; icomp++) {
1622
1623 for (i = 0; i < dij; i++) { buf[i] = 0; }
1624 for (l = 0; l < dl; l++) {
1625 for (k = 0; k < dk; k++) {
1626 for (ij = 0; ij < dij; ij++) {
1627 buf[ij] += eri[ij] * dmlk[l*dk+k];
1628 }
1629 eri += dij;
1630 } }
1631
1632 for (ij = 0, j = 0; j < dj; j++) {
1633 for (i = 0; i < di; i++, ij++) {
1634 vij[i*dj+j] += buf[ij];
1635 vji[ij] -= buf[ij];
1636 } }
1637 vij += dij;
1638 vji += dij;
1639 }
1640 } else {
1641 nrs1_lk_s1ij (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1642 }
1643 }
1644 ADD_JKOP(nra2ij_lk_s1ij, L, K, I, J, s2ij);
1645
nra2ij_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1646 static void nra2ij_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
1647 int i0, int i1, int j0, int j1,
1648 int k0, int k1, int l0, int l1)
1649 {
1650 nra2ij_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1651 }
1652 ADD_JKOP(nra2ij_lk_s2ij, L, K, I, J, s2ij);
1653
nra2ij_lk_a2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1654 static void nra2ij_lk_a2ij(double *eri, double *dm, JKArray *out, int *shls,
1655 int i0, int i1, int j0, int j1,
1656 int k0, int k1, int l0, int l1)
1657 {
1658 nra2ij_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1659 }
1660 ADD_JKOP(nra2ij_lk_a2ij, L, K, I, J, s2ij);
1661
nra2ij_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1662 static void nra2ij_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
1663 int i0, int i1, int j0, int j1,
1664 int k0, int k1, int l0, int l1)
1665 {
1666 if (i0 > j0) {
1667 DECLARE(vil, i, l);
1668 DEF_DM(i, k);
1669 DEF_DM(j, k);
1670 LOCATE(vjl, j, l);
1671 int i, j, k, l, ijkl, icomp;
1672
1673 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1674 for (l = 0; l < dl; l++) {
1675 for (k = 0; k < dk; k++) {
1676 for (j = 0; j < dj; j++) {
1677 for (i = 0; i < di; i++, ijkl++) {
1678 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
1679 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
1680 } } } }
1681 vil += dil;
1682 vjl += djl;
1683 }
1684 } else {
1685 nrs1_jk_s1il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1686 }
1687 }
1688 ADD_JKOP(nra2ij_jk_s1il, J, K, I, L, s2ij);
1689
nra2ij_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1690 static void nra2ij_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
1691 int i0, int i1, int j0, int j1,
1692 int k0, int k1, int l0, int l1)
1693 {
1694 if (j0 >= l0) {
1695 nra2ij_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1696 } else {
1697 nrs1_jk_s2il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1698 }
1699 }
1700 ADD_JKOP(nra2ij_jk_s2il, J, K, I, L, s2ij);
1701
nra2ij_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1702 static void nra2ij_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
1703 int i0, int i1, int j0, int j1,
1704 int k0, int k1, int l0, int l1)
1705 {
1706 if (i0 > j0) {
1707 DECLARE(vkj, k, j);
1708 DEF_DM(l, i);
1709 DEF_DM(l, j);
1710 LOCATE(vki, k, i);
1711 int i, j, k, l, ijkl, icomp;
1712
1713 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1714 for (l = 0; l < dl; l++) {
1715 for (k = 0; k < dk; k++) {
1716 for (j = 0; j < dj; j++) {
1717 for (i = 0; i < di; i++, ijkl++) {
1718 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
1719 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
1720 } } } }
1721 vkj += dkj;
1722 vki += dki;
1723 }
1724 } else {
1725 nrs1_li_s1kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1726 }
1727 }
1728 ADD_JKOP(nra2ij_li_s1kj, L, I, K, J, s2ij);
1729
nra2ij_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1730 static void nra2ij_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
1731 int i0, int i1, int j0, int j1,
1732 int k0, int k1, int l0, int l1)
1733 {
1734 if (k0 >= i0) {
1735 nra2ij_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1736 } else {
1737 nrs1_li_s2kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1738 }
1739 }
1740 ADD_JKOP(nra2ij_li_s2kj, L, I, K, J, s2ij);
1741
nra2kl_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1742 static void nra2kl_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
1743 int i0, int i1, int j0, int j1,
1744 int k0, int k1, int l0, int l1)
1745 {
1746 if (k0 > l0) {
1747 DECLARE(vkl, k, l);
1748 LOCATE(vlk, l, k);
1749 DEF_DM(j, i);
1750 int dij = di * dj;
1751 int k, l, ij, icomp;
1752 double tmp;
1753
1754 for (icomp = 0; icomp < ncomp; icomp++) {
1755 for (l = 0; l < dl; l++) {
1756 for (k = 0; k < dk; k++) {
1757 tmp = 0;
1758 for (ij = 0; ij < dij; ij++) {
1759 tmp += eri[ij] * dmji[ij];
1760 }
1761 vkl[k*dl+l] += tmp;
1762 vlk[l*dk+k] -= tmp;
1763 eri += dij;
1764 } }
1765 vkl += dkl;
1766 vlk += dkl;
1767 }
1768 } else {
1769 nrs1_ji_s1kl (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1770 }
1771 }
1772 ADD_JKOP(nra2kl_ji_s1kl, J, I, K, L, s2kl);
1773
nra2kl_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1774 static void nra2kl_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
1775 int i0, int i1, int j0, int j1,
1776 int k0, int k1, int l0, int l1)
1777 {
1778 nrs1_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1779 }
1780 ADD_JKOP(nra2kl_ji_s2kl, J, I, K, L, s2kl);
1781
nra2kl_ji_a2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1782 static void nra2kl_ji_a2kl(double *eri, double *dm, JKArray *out, int *shls,
1783 int i0, int i1, int j0, int j1,
1784 int k0, int k1, int l0, int l1)
1785 {
1786 nrs1_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1787 }
1788 ADD_JKOP(nra2kl_ji_a2kl, J, I, K, L, s2kl);
1789
nra2kl_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1790 static void nra2kl_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
1791 int i0, int i1, int j0, int j1,
1792 int k0, int k1, int l0, int l1)
1793 {
1794 if (k0 > l0) {
1795 DECLARE(v, i, j);
1796 DEF_DM(k, l);
1797 DEF_DM(l, k);
1798 int i, j, k, l, ij, icomp;
1799 double buf[MAXCGTO*MAXCGTO];
1800 double tdm;
1801
1802 for (icomp = 0; icomp < ncomp; icomp++) {
1803
1804 for (i = 0; i < dij; i++) { buf[i] = 0; }
1805 for (l = 0; l < dl; l++) {
1806 for (k = 0; k < dk; k++) {
1807 tdm = dmlk[l*dk+k] - dmkl[k*dl+l];
1808 for (ij = 0; ij < dij; ij++) {
1809 buf[ij] += eri[ij] * tdm;
1810 }
1811 eri += dij;
1812 } }
1813
1814 for (ij = 0, j = 0; j < dj; j++) {
1815 for (i = 0; i < di; i++, ij++) {
1816 v[i*dj+j] += buf[ij];
1817 } }
1818 v += dij;
1819 }
1820 } else {
1821 nrs1_lk_s1ij (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1822 }
1823 }
1824 ADD_JKOP(nra2kl_lk_s1ij, L, K, I, J, s2kl);
1825
nra2kl_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1826 static void nra2kl_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
1827 int i0, int i1, int j0, int j1,
1828 int k0, int k1, int l0, int l1)
1829 {
1830 if (i0 >= j0) {
1831 nra2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1832 }
1833 }
1834 ADD_JKOP(nra2kl_lk_s2ij, L, K, I, J, s2kl);
1835
nra2kl_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1836 static void nra2kl_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
1837 int i0, int i1, int j0, int j1,
1838 int k0, int k1, int l0, int l1)
1839 {
1840 if (k0 > l0) {
1841 DECLARE(vil, i, l);
1842 DEF_DM(j, k);
1843 DEF_DM(j, l);
1844 LOCATE(vik, i, k);
1845 int i, j, k, l, ijkl, icomp;
1846
1847 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1848 for (l = 0; l < dl; l++) {
1849 for (k = 0; k < dk; k++) {
1850 for (j = 0; j < dj; j++) {
1851 for (i = 0; i < di; i++, ijkl++) {
1852 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
1853 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
1854 } } } }
1855 vik += dik;
1856 vil += dil;
1857 }
1858 } else {
1859 nrs1_jk_s1il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1860 }
1861 }
1862 ADD_JKOP(nra2kl_jk_s1il, J, K, I, L, s2kl);
1863
nra2kl_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1864 static void nra2kl_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
1865 int i0, int i1, int j0, int j1,
1866 int k0, int k1, int l0, int l1)
1867 {
1868 if (i0 >= k0) {
1869 nra2kl_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1870 } else {
1871 nrs1_jk_s2il (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1872 }
1873 }
1874 ADD_JKOP(nra2kl_jk_s2il, J, K, I, L, s2kl);
1875
nra2kl_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1876 static void nra2kl_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
1877 int i0, int i1, int j0, int j1,
1878 int k0, int k1, int l0, int l1)
1879 {
1880 if (k0 > l0) {
1881 DECLARE(vkj, k, j);
1882 DEF_DM(l, i);
1883 DEF_DM(k, i);
1884 LOCATE(vlj, l, j);
1885 int i, j, k, l, ijkl, icomp;
1886
1887 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
1888 for (l = 0; l < dl; l++) {
1889 for (k = 0; k < dk; k++) {
1890 for (j = 0; j < dj; j++) {
1891 for (i = 0; i < di; i++, ijkl++) {
1892 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
1893 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
1894 } } } }
1895 vkj += dkj;
1896 vlj += dlj;
1897 }
1898 } else {
1899 nrs1_li_s1kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1900 }
1901 }
1902 ADD_JKOP(nra2kl_li_s1kj, L, I, K, J, s2kl);
1903
nra2kl_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1904 static void nra2kl_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
1905 int i0, int i1, int j0, int j1,
1906 int k0, int k1, int l0, int l1)
1907 {
1908 if (l0 >= j0) {
1909 nra2kl_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1910 } else {
1911 nrs1_li_s2kj (eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1912 }
1913 }
1914 ADD_JKOP(nra2kl_li_s2kj, L, I, K, J, s2kl);
1915
nra4ij_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1916 static void nra4ij_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
1917 int i0, int i1, int j0, int j1,
1918 int k0, int k1, int l0, int l1)
1919 {
1920 if (i0 == j0) {
1921 nrs2kl_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1922 } else if (k0 == l0) {
1923 nra2ij_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1924 } else {
1925 DECLARE(vkl, k, l);
1926 LOCATE(vlk, l, k);
1927 DEF_DM(i, j);
1928 DEF_DM(j, i);
1929 int dij = di * dj;
1930 int i, j, k, l, ij, icomp;
1931 double tdm[MAXCGTO*MAXCGTO];
1932 double tmp;
1933
1934 for (ij = 0, j = 0; j < dj; j++) {
1935 for (i = 0; i < di; i++, ij++) {
1936 tdm[ij] = dmji[j*di+i] - dmij[i*dj+j];
1937 }
1938 }
1939
1940 for (icomp = 0; icomp < ncomp; icomp++) {
1941 for (l = 0; l < dl; l++) {
1942 for (k = 0; k < dk; k++) {
1943 tmp = 0;
1944 for (ij = 0; ij < dij; ij++) {
1945 tmp += eri[ij] * tdm[ij];
1946 }
1947 vkl[k*dl+l] += tmp;
1948 vlk[l*dk+k] += tmp;
1949 eri += dij;
1950 } }
1951 vkl += dkl;
1952 vlk += dkl;
1953 }
1954 }
1955 }
1956 ADD_JKOP(nra4ij_ji_s1kl, J, I, K, L, s4);
1957
nra4ij_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1958 static void nra4ij_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
1959 int i0, int i1, int j0, int j1,
1960 int k0, int k1, int l0, int l1)
1961 {
1962 nra2ij_ji_s2kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1963 }
1964 ADD_JKOP(nra4ij_ji_s2kl, J, I, K, L, s4);
1965
nra4ij_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)1966 static void nra4ij_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
1967 int i0, int i1, int j0, int j1,
1968 int k0, int k1, int l0, int l1)
1969 {
1970 if (i0 == j0) {
1971 nrs2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1972 } else if (k0 == l0) {
1973 nra2ij_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
1974 } else {
1975 DECLARE(vij, i, j);
1976 LOCATE(vji, j, i);
1977 DEF_DM(k, l);
1978 DEF_DM(l, k);
1979 int i, j, k, l, ij, icomp;
1980 double buf[MAXCGTO*MAXCGTO];
1981 double tdm;
1982
1983 for (icomp = 0; icomp < ncomp; icomp++) {
1984
1985 for (i = 0; i < dij; i++) { buf[i] = 0; }
1986 for (l = 0; l < dl; l++) {
1987 for (k = 0; k < dk; k++) {
1988 tdm = dmlk[l*dk+k] + dmkl[k*dl+l];
1989 for (ij = 0; ij < dij; ij++) {
1990 buf[ij] += eri[ij] * tdm;
1991 }
1992 eri += dij;
1993 } }
1994
1995 for (ij = 0, j = 0; j < dj; j++) {
1996 for (i = 0; i < di; i++, ij++) {
1997 vij[i*dj+j] += buf[ij];
1998 vji[ij] -= buf[ij];
1999 } }
2000 vij += dij;
2001 vji += dij;
2002 }
2003 }
2004 }
2005 ADD_JKOP(nra4ij_lk_s1ij, L, K, I, J, s4);
2006
nra4ij_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2007 static void nra4ij_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
2008 int i0, int i1, int j0, int j1,
2009 int k0, int k1, int l0, int l1)
2010 {
2011 nrs2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2012 }
2013 ADD_JKOP(nra4ij_lk_s2ij, L, K, I, J, s4);
2014
nra4ij_lk_a2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2015 static void nra4ij_lk_a2ij(double *eri, double *dm, JKArray *out, int *shls,
2016 int i0, int i1, int j0, int j1,
2017 int k0, int k1, int l0, int l1)
2018 {
2019 nrs2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2020 }
2021 ADD_JKOP(nra4ij_lk_a2ij, L, K, I, J, s4);
2022
nra4ij_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2023 static void nra4ij_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
2024 int i0, int i1, int j0, int j1,
2025 int k0, int k1, int l0, int l1)
2026 {
2027 if (i0 == j0) {
2028 nrs2kl_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2029 } else if (k0 == l0) {
2030 nra2ij_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2031 } else {
2032 DECLARE(vik, i, k);
2033 LOCATE(vil, i, l);
2034 LOCATE(vjk, j, k);
2035 LOCATE(vjl, j, l);
2036 DEF_DM(i, l);
2037 DEF_DM(i, k);
2038 DEF_DM(j, l);
2039 DEF_DM(j, k);
2040 int i, j, k, l, ijkl, icomp;
2041
2042 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2043 for (l = 0; l < dl; l++) {
2044 for (k = 0; k < dk; k++) {
2045 for (j = 0; j < dj; j++) {
2046 for (i = 0; i < di; i++, ijkl++) {
2047 vjk[j*dk+k] -= eri[ijkl] * dmil[i*dl+l];
2048 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
2049 vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
2050 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2051 } } } }
2052 vjk += djk;
2053 vjl += djl;
2054 vik += dik;
2055 vil += dil;
2056 }
2057 }
2058 }
2059 ADD_JKOP(nra4ij_jk_s1il, J, K, I, L, s4);
2060
nra4ij_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2061 static void nra4ij_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
2062 int i0, int i1, int j0, int j1,
2063 int k0, int k1, int l0, int l1)
2064 {
2065 if (i0 == j0) {
2066 nrs2kl_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2067 } else if (k0 == l0) {
2068 nra2ij_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2069 } else if (i0 < l0) {
2070 } else if (i0 < k0) {
2071 if (j0 < l0) { // j < l <= i < k
2072 DECLARE(v, i, l);
2073 DEF_DM(j, k);
2074 int i, j, k, l, ijkl, icomp;
2075 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2076 for (l = 0; l < dl; l++) {
2077 for (k = 0; k < dk; k++) {
2078 for (j = 0; j < dj; j++) {
2079 for (i = 0; i < di; i++, ijkl++) {
2080 v[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2081 } } } }
2082 v += dil;
2083 }
2084 } else { // l <= j < i < k
2085 DECLARE(vil, i, l);
2086 DEF_DM(i, k);
2087 DEF_DM(j, k);
2088 LOCATE(vjl, j, l);
2089 int i, j, k, l, ijkl, icomp;
2090 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2091 for (l = 0; l < dl; l++) {
2092 for (k = 0; k < dk; k++) {
2093 for (j = 0; j < dj; j++) {
2094 for (i = 0; i < di; i++, ijkl++) {
2095 vjl[j*dl+l] -= eri[ijkl] *dmik[i*dk+k];
2096 vil[i*dl+l] += eri[ijkl] *dmjk[j*dk+k];
2097 } } } }
2098 vjl += djl;
2099 vil += dil;
2100 }
2101 }
2102 } else if (j0 < l0) { // j < l < k <= i
2103 DECLARE(vil, i, l);
2104 DEF_DM(j, k);
2105 DEF_DM(j, l);
2106 LOCATE(vik, i, k);
2107 int i, j, k, l, ijkl, icomp;
2108 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2109 for (l = 0; l < dl; l++) {
2110 for (k = 0; k < dk; k++) {
2111 for (j = 0; j < dj; j++) {
2112 for (i = 0; i < di; i++, ijkl++) {
2113 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2114 vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
2115 } } } }
2116 vil += dil;
2117 vik += dik;
2118 }
2119 } else if (j0 < k0) { // l <= j < k <= i
2120 DECLARE(vjl, j, l);
2121 LOCATE(vik, i, k);
2122 LOCATE(vil, i, l);
2123 DEF_DM(i, k);
2124 DEF_DM(j, k);
2125 DEF_DM(j, l);
2126 int i, j, k, l, ijkl, icomp;
2127 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2128 for (l = 0; l < dl; l++) {
2129 for (k = 0; k < dk; k++) {
2130 for (j = 0; j < dj; j++) {
2131 for (i = 0; i < di; i++, ijkl++) {
2132 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
2133 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2134 vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
2135 } } } }
2136 vjl += djl;
2137 vil += dil;
2138 vik += dik;
2139 }
2140 } else { // l < k <= j < i
2141 DECLARE(vjl, j, l);
2142 LOCATE(vik, i, k);
2143 LOCATE(vjk, j, k);
2144 LOCATE(vil, i, l);
2145 DEF_DM(i, l);
2146 DEF_DM(i, k);
2147 DEF_DM(j, l);
2148 DEF_DM(j, k);
2149 int i, j, k, l, ijkl, icomp;
2150 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2151 for (l = 0; l < dl; l++) {
2152 for (k = 0; k < dk; k++) {
2153 for (j = 0; j < dj; j++) {
2154 for (i = 0; i < di; i++, ijkl++) {
2155 vjk[j*dk+k] -= eri[ijkl] * dmil[i*dl+l];
2156 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
2157 vik[i*dk+k] += eri[ijkl] * dmjl[j*dl+l];
2158 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2159 } } } }
2160 vjk += djk;
2161 vjl += djl;
2162 vik += dik;
2163 vil += dil;
2164 }
2165 }
2166 }
2167 ADD_JKOP(nra4ij_jk_s2il, J, K, I, L, s4);
2168
nra4ij_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2169 static void nra4ij_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
2170 int i0, int i1, int j0, int j1,
2171 int k0, int k1, int l0, int l1)
2172 {
2173 if (i0 == j0) {
2174 nrs2kl_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2175 } else if (k0 == l0) {
2176 nra2ij_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2177 } else {
2178 DECLARE(vki, k, i);
2179 LOCATE(vkj, k, j);
2180 LOCATE(vli, l, i);
2181 LOCATE(vlj, l, j);
2182 DEF_DM(l, i);
2183 DEF_DM(l, j);
2184 DEF_DM(k, i);
2185 DEF_DM(k, j);
2186 int i, j, k, l, ijkl, icomp;
2187 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2188 for (l = 0; l < dl; l++) {
2189 for (k = 0; k < dk; k++) {
2190 for (j = 0; j < dj; j++) {
2191 for (i = 0; i < di; i++, ijkl++) {
2192 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2193 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
2194 vlj[l*dj+j] += eri[ijkl] * dmki[k*di+i];
2195 vli[l*di+i] -= eri[ijkl] * dmkj[k*dj+j];
2196 } } } }
2197 vkj += dkj;
2198 vki += dki;
2199 vlj += dlj;
2200 vli += dli;
2201 }
2202 }
2203 }
2204 ADD_JKOP(nra4ij_li_s1kj, L, I, K, J, s4);
2205
nra4ij_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2206 static void nra4ij_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
2207 int i0, int i1, int j0, int j1,
2208 int k0, int k1, int l0, int l1)
2209 {
2210 if (i0 == j0) {
2211 nrs2kl_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2212 } else if (k0 == l0) {
2213 nra2ij_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2214 } else if (k0 < j0) {
2215 } else if (k0 < i0) {
2216 if (l0 < j0) { // l < j < k < i
2217 DECLARE(v, k, j);
2218 DEF_DM(l, i);
2219 int i, j, k, l, ijkl, icomp;
2220 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2221 for (l = 0; l < dl; l++) {
2222 for (k = 0; k < dk; k++) {
2223 for (j = 0; j < dj; j++) {
2224 for (i = 0; i < di; i++, ijkl++) {
2225 v[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2226 } } } }
2227 v += dkj;
2228 }
2229 } else { // j <= l < k < i
2230 DECLARE(vkj, k, j);
2231 LOCATE(vlj, l, j);
2232 DEF_DM(l, i);
2233 DEF_DM(k, i);
2234 int i, j, k, l, ijkl, icomp;
2235 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2236 for (l = 0; l < dl; l++) {
2237 for (k = 0; k < dk; k++) {
2238 for (j = 0; j < dj; j++) {
2239 for (i = 0; i < di; i++, ijkl++) {
2240 vkj[k*dj+j] += eri[ijkl] *dmli[l*di+i];
2241 vlj[l*dj+j] += eri[ijkl] *dmki[k*di+i];
2242 } } } }
2243 vkj += dkj;
2244 vlj += dlj;
2245 }
2246 }
2247 } else if (l0 < j0) { // l < j < i <= k
2248 DECLARE(vki, k, i);
2249 LOCATE(vkj, k, j);
2250 DEF_DM(l, i);
2251 DEF_DM(l, j);
2252 int i, j, k, l, ijkl, icomp;
2253 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2254 for (l = 0; l < dl; l++) {
2255 for (k = 0; k < dk; k++) {
2256 for (j = 0; j < dj; j++) {
2257 for (i = 0; i < di; i++, ijkl++) {
2258 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2259 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
2260 } } } }
2261 vkj += dkj;
2262 vki += dki;
2263 }
2264 } else if (l0 < i0) { // j <= l < i <= k
2265 DECLARE(vki, k, i);
2266 LOCATE(vkj, k, j);
2267 LOCATE(vlj, l, j);
2268 DEF_DM(l, i);
2269 DEF_DM(l, j);
2270 DEF_DM(k, i);
2271 int i, j, k, l, ijkl, icomp;
2272 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2273 for (l = 0; l < dl; l++) {
2274 for (k = 0; k < dk; k++) {
2275 for (j = 0; j < dj; j++) {
2276 for (i = 0; i < di; i++, ijkl++) {
2277 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2278 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
2279 vlj[l*dj+j] += eri[ijkl] * dmki[k*di+i];
2280 } } } }
2281 vkj += dkj;
2282 vki += dki;
2283 vlj += dlj;
2284 }
2285 } else { // j < i <= l < k
2286 DECLARE(vki, k, i);
2287 LOCATE(vkj, k, j);
2288 LOCATE(vli, l, i);
2289 LOCATE(vlj, l, j);
2290 DEF_DM(l, i);
2291 DEF_DM(l, j);
2292 DEF_DM(k, i);
2293 DEF_DM(k, j);
2294 int i, j, k, l, ijkl, icomp;
2295 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2296 for (l = 0; l < dl; l++) {
2297 for (k = 0; k < dk; k++) {
2298 for (j = 0; j < dj; j++) {
2299 for (i = 0; i < di; i++, ijkl++) {
2300 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2301 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
2302 vlj[l*dj+j] += eri[ijkl] * dmki[k*di+i];
2303 vli[l*di+i] -= eri[ijkl] * dmkj[k*dj+j];
2304 } } } }
2305 vkj += dkj;
2306 vki += dki;
2307 vlj += dlj;
2308 vli += dli;
2309 }
2310 }
2311 }
2312 ADD_JKOP(nra4ij_li_s2kj, L, I, K, J, s4);
2313
nra4kl_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2314 static void nra4kl_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
2315 int i0, int i1, int j0, int j1,
2316 int k0, int k1, int l0, int l1)
2317 {
2318 if (i0 == j0) {
2319 nra2kl_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2320 } else if (k0 == l0) {
2321 nrs2ij_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2322 } else {
2323 DECLARE(vkl, k, l);
2324 LOCATE(vlk, l, k);
2325 DEF_DM(i, j);
2326 DEF_DM(j, i);
2327 int dij = di * dj;
2328 int i, j, k, l, ij, icomp;
2329 double tdm[MAXCGTO*MAXCGTO];
2330 double tmp;
2331
2332 for (ij = 0, j = 0; j < dj; j++) {
2333 for (i = 0; i < di; i++, ij++) {
2334 tdm[ij] = dmij[i*dj+j] + dmji[j*di+i];
2335 }
2336 }
2337
2338 for (icomp = 0; icomp < ncomp; icomp++) {
2339 for (l = 0; l < dl; l++) {
2340 for (k = 0; k < dk; k++) {
2341 tmp = 0;
2342 for (ij = 0; ij < dij; ij++) {
2343 tmp += eri[ij] * tdm[ij];
2344 }
2345 vkl[k*dl+l] += tmp;
2346 vlk[l*dk+k] -= tmp;
2347 eri += dij;
2348 } }
2349 vkl += dkl;
2350 vlk += dkl;
2351 }
2352 }
2353 }
2354 ADD_JKOP(nra4kl_ji_s1kl, J, I, K, L, s4);
2355
nra4kl_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2356 static void nra4kl_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
2357 int i0, int i1, int j0, int j1,
2358 int k0, int k1, int l0, int l1)
2359 {
2360 nra2kl_ji_s2kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2361 }
2362 ADD_JKOP(nra4kl_ji_s2kl, J, I, K, L, s4);
2363
nra4kl_ji_a2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2364 static void nra4kl_ji_a2kl(double *eri, double *dm, JKArray *out, int *shls,
2365 int i0, int i1, int j0, int j1,
2366 int k0, int k1, int l0, int l1)
2367 {
2368 nra2kl_ji_s2kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2369 }
2370 ADD_JKOP(nra4kl_ji_a2kl, J, I, K, L, s4);
2371
nra4kl_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2372 static void nra4kl_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
2373 int i0, int i1, int j0, int j1,
2374 int k0, int k1, int l0, int l1)
2375 {
2376 if (i0 == j0) {
2377 nra2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2378 } else if (k0 == l0) {
2379 nrs2ij_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2380 } else {
2381 DECLARE(vij, i, j);
2382 LOCATE(vji, j, i);
2383 DEF_DM(k, l);
2384 DEF_DM(l, k);
2385 int i, j, k, l, ij, icomp;
2386 double buf[MAXCGTO*MAXCGTO];
2387 double tdm;
2388
2389 for (icomp = 0; icomp < ncomp; icomp++) {
2390
2391 for (i = 0; i < dij; i++) { buf[i] = 0; }
2392 for (l = 0; l < dl; l++) {
2393 for (k = 0; k < dk; k++) {
2394 tdm = dmlk[l*dk+k] - dmkl[k*dl+l];
2395 for (ij = 0; ij < dij; ij++) {
2396 buf[ij] += eri[ij] * tdm;
2397 }
2398 eri += dij;
2399 } }
2400
2401 for (ij = 0, j = 0; j < dj; j++) {
2402 for (i = 0; i < di; i++, ij++) {
2403 vij[i*dj+j] += buf[ij];
2404 vji[ij] += buf[ij];
2405 } }
2406 vij += dij;
2407 vji += dij;
2408 }
2409 }
2410 }
2411 ADD_JKOP(nra4kl_lk_s1ij, L, K, I, J, s4);
2412
nra4kl_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2413 static void nra4kl_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
2414 int i0, int i1, int j0, int j1,
2415 int k0, int k1, int l0, int l1)
2416 {
2417 nra2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2418 }
2419 ADD_JKOP(nra4kl_lk_s2ij, L, K, I, J, s4);
2420
nra4kl_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2421 static void nra4kl_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
2422 int i0, int i1, int j0, int j1,
2423 int k0, int k1, int l0, int l1)
2424 {
2425 if (i0 == j0) {
2426 nra2kl_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2427 } else if (k0 == l0) {
2428 nrs2ij_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2429 } else {
2430 DECLARE(vik, i, k);
2431 LOCATE(vil, i, l);
2432 LOCATE(vjk, j, k);
2433 LOCATE(vjl, j, l);
2434 DEF_DM(i, l);
2435 DEF_DM(i, k);
2436 DEF_DM(j, l);
2437 DEF_DM(j, k);
2438 int i, j, k, l, ijkl, icomp;
2439
2440 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2441 for (l = 0; l < dl; l++) {
2442 for (k = 0; k < dk; k++) {
2443 for (j = 0; j < dj; j++) {
2444 for (i = 0; i < di; i++, ijkl++) {
2445 vjk[j*dk+k] -= eri[ijkl] * dmil[i*dl+l];
2446 vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
2447 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2448 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2449 } } } }
2450 vjk += djk;
2451 vjl += djl;
2452 vik += dik;
2453 vil += dil;
2454 }
2455 }
2456 }
2457 ADD_JKOP(nra4kl_jk_s1il, J, K, I, L, s4);
2458
nra4kl_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2459 static void nra4kl_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
2460 int i0, int i1, int j0, int j1,
2461 int k0, int k1, int l0, int l1)
2462 {
2463 if (i0 == j0) {
2464 nra2kl_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2465 } else if (k0 == l0) {
2466 nrs2ij_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2467 } else if (i0 < l0) {
2468 } else if (i0 < k0) {
2469 if (j0 < l0) { // j < l <= i < k
2470 DECLARE(v, i, l);
2471 DEF_DM(j, k);
2472 int i, j, k, l, ijkl, icomp;
2473 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2474 for (l = 0; l < dl; l++) {
2475 for (k = 0; k < dk; k++) {
2476 for (j = 0; j < dj; j++) {
2477 for (i = 0; i < di; i++, ijkl++) {
2478 v[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2479 } } } }
2480 v += dil;
2481 }
2482 } else { // l <= j < i < k
2483 DECLARE(vil, i, l);
2484 LOCATE(vjl, j, l);
2485 DEF_DM(i, k);
2486 DEF_DM(j, k);
2487 int i, j, k, l, ijkl, icomp;
2488 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2489 for (l = 0; l < dl; l++) {
2490 for (k = 0; k < dk; k++) {
2491 for (j = 0; j < dj; j++) {
2492 for (i = 0; i < di; i++, ijkl++) {
2493 vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
2494 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2495 } } } }
2496 vjl += djl;
2497 vil += dil;
2498 }
2499 }
2500 } else if (j0 < l0) { // j < l < k <= i
2501 DECLARE(vil, i, l);
2502 LOCATE(vik, i, k);
2503 DEF_DM(j, k);
2504 DEF_DM(j, l);
2505 int i, j, k, l, ijkl, icomp;
2506 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2507 for (l = 0; l < dl; l++) {
2508 for (k = 0; k < dk; k++) {
2509 for (j = 0; j < dj; j++) {
2510 for (i = 0; i < di; i++, ijkl++) {
2511 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2512 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2513 } } } }
2514 vil += dil;
2515 vik += dik;
2516 }
2517 } else if (j0 < k0) { // l <= j < k <= i
2518 DECLARE(vjl, j, l);
2519 LOCATE(vik, i, k);
2520 LOCATE(vil, i, l);
2521 DEF_DM(i, k);
2522 DEF_DM(j, k);
2523 DEF_DM(j, l);
2524 int i, j, k, l, ijkl, icomp;
2525 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2526 for (l = 0; l < dl; l++) {
2527 for (k = 0; k < dk; k++) {
2528 for (j = 0; j < dj; j++) {
2529 for (i = 0; i < di; i++, ijkl++) {
2530 vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
2531 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2532 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2533 } } } }
2534 vjl += djl;
2535 vil += dil;
2536 vik += dik;
2537 }
2538 } else { // l < k <= j < i
2539 DECLARE(vjl, j, l);
2540 LOCATE(vik, i, k);
2541 LOCATE(vjk, j, k);
2542 LOCATE(vil, i, l);
2543 DEF_DM(i, l);
2544 DEF_DM(i, k);
2545 DEF_DM(j, l);
2546 DEF_DM(j, k);
2547 int i, j, k, l, ijkl, icomp;
2548 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2549 for (l = 0; l < dl; l++) {
2550 for (k = 0; k < dk; k++) {
2551 for (j = 0; j < dj; j++) {
2552 for (i = 0; i < di; i++, ijkl++) {
2553 vjk[j*dk+k] -= eri[ijkl] * dmil[i*dl+l];
2554 vjl[j*dl+l] += eri[ijkl] * dmik[i*dk+k];
2555 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2556 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2557 } } } }
2558 vjk += djk;
2559 vjl += djl;
2560 vik += dik;
2561 vil += dil;
2562 }
2563 }
2564 }
2565 ADD_JKOP(nra4kl_jk_s2il, J, K, I, L, s4);
2566
nra4kl_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2567 static void nra4kl_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
2568 int i0, int i1, int j0, int j1,
2569 int k0, int k1, int l0, int l1)
2570 {
2571 if (i0 == j0) {
2572 nra2kl_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2573 } else if (k0 == l0) {
2574 nrs2ij_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2575 } else {
2576 DECLARE(vki, k, i);
2577 LOCATE(vkj, k, j);
2578 LOCATE(vli, l, i);
2579 LOCATE(vlj, l, j);
2580 DEF_DM(l, i);
2581 DEF_DM(l, j);
2582 DEF_DM(k, i);
2583 DEF_DM(k, j);
2584 int i, j, k, l, ijkl, icomp;
2585 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2586 for (l = 0; l < dl; l++) {
2587 for (k = 0; k < dk; k++) {
2588 for (j = 0; j < dj; j++) {
2589 for (i = 0; i < di; i++, ijkl++) {
2590 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2591 vki[k*di+i] += eri[ijkl] * dmlj[l*dj+j];
2592 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
2593 vli[l*di+i] -= eri[ijkl] * dmkj[k*dj+j];
2594 } } } }
2595 vkj += dkj;
2596 vki += dki;
2597 vlj += dlj;
2598 vli += dli;
2599 }
2600 }
2601 }
2602 ADD_JKOP(nra4kl_li_s1kj, L, I, K, J, s4);
2603
nra4kl_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2604 static void nra4kl_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
2605 int i0, int i1, int j0, int j1,
2606 int k0, int k1, int l0, int l1)
2607 {
2608 if (i0 == j0) {
2609 nrs2kl_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2610 } else if (k0 == l0) {
2611 nrs2ij_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2612 } else if (k0 < j0) {
2613 } else if (k0 < i0) {
2614 if (l0 < j0) { // l < j < k < i
2615 DECLARE(v, k, j);
2616 DEF_DM(l, i);
2617 int i, j, k, l, ijkl, icomp;
2618 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2619 for (l = 0; l < dl; l++) {
2620 for (k = 0; k < dk; k++) {
2621 for (j = 0; j < dj; j++) {
2622 for (i = 0; i < di; i++, ijkl++) {
2623 v[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2624 } } } }
2625 v += dkj;
2626 }
2627 } else { // j <= l < k < i
2628 DECLARE(vkj, k, j);
2629 LOCATE(vlj, l, j);
2630 DEF_DM(l, i);
2631 DEF_DM(k, i);
2632 int i, j, k, l, ijkl, icomp;
2633 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2634 for (l = 0; l < dl; l++) {
2635 for (k = 0; k < dk; k++) {
2636 for (j = 0; j < dj; j++) {
2637 for (i = 0; i < di; i++, ijkl++) {
2638 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2639 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
2640 } } } }
2641 vkj += dkj;
2642 vlj += dlj;
2643 }
2644 }
2645 } else if (l0 < j0) { // l < j < i <= k
2646 DECLARE(vki, k, i);
2647 LOCATE(vkj, k, j);
2648 DEF_DM(l, i);
2649 DEF_DM(l, j);
2650 int i, j, k, l, ijkl, icomp;
2651 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2652 for (l = 0; l < dl; l++) {
2653 for (k = 0; k < dk; k++) {
2654 for (j = 0; j < dj; j++) {
2655 for (i = 0; i < di; i++, ijkl++) {
2656 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2657 vki[k*di+i] += eri[ijkl] * dmlj[l*dj+j];
2658 } } } }
2659 vkj += dkj;
2660 vki += dki;
2661 }
2662 } else if (l0 < i0) { // j <= l < i <= k
2663 DECLARE(vki, k, i);
2664 LOCATE(vkj, k, j);
2665 LOCATE(vlj, l, j);
2666 DEF_DM(l, i);
2667 DEF_DM(l, j);
2668 DEF_DM(k, i);
2669 int i, j, k, l, ijkl, icomp;
2670 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2671 for (l = 0; l < dl; l++) {
2672 for (k = 0; k < dk; k++) {
2673 for (j = 0; j < dj; j++) {
2674 for (i = 0; i < di; i++, ijkl++) {
2675 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2676 vki[k*di+i] += eri[ijkl] * dmlj[l*dj+j];
2677 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
2678 } } } }
2679 vkj += dkj;
2680 vki += dki;
2681 vlj += dlj;
2682 }
2683 } else { // j < i <= l < k
2684 DECLARE(vki, k, i);
2685 LOCATE(vkj, k, j);
2686 LOCATE(vli, l, i);
2687 LOCATE(vlj, l, j);
2688 DEF_DM(l, i);
2689 DEF_DM(l, j);
2690 DEF_DM(k, i);
2691 DEF_DM(k, j);
2692 int i, j, k, l, ijkl, icomp;
2693 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2694 for (l = 0; l < dl; l++) {
2695 for (k = 0; k < dk; k++) {
2696 for (j = 0; j < dj; j++) {
2697 for (i = 0; i < di; i++, ijkl++) {
2698 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2699 vki[k*di+i] += eri[ijkl] * dmlj[l*dj+j];
2700 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
2701 vli[l*di+i] -= eri[ijkl] * dmkj[k*dj+j];
2702 } } } }
2703 vkj += dkj;
2704 vki += dki;
2705 vlj += dlj;
2706 vli += dli;
2707 }
2708 }
2709 }
2710 ADD_JKOP(nra4kl_li_s2kj, L, I, K, J, s4);
2711
2712 /*
2713 * aa4: 4-fold permutation symmetry with anti-symm for ij and anti-symm for kl
2714 */
nraa4_ji_s1kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2715 static void nraa4_ji_s1kl(double *eri, double *dm, JKArray *out, int *shls,
2716 int i0, int i1, int j0, int j1,
2717 int k0, int k1, int l0, int l1)
2718 {
2719 if (i0 == j0) {
2720 nra2kl_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2721 } else if (k0 == l0) {
2722 nra2ij_ji_s1kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2723 } else {
2724 DECLARE(vkl, k, l);
2725 LOCATE(vlk, l, k);
2726 DEF_DM(i, j);
2727 DEF_DM(j, i);
2728 int dij = di * dj;
2729 int i, j, k, l, ij, icomp;
2730 double tdm[MAXCGTO*MAXCGTO];
2731 double tmp;
2732
2733 for (ij = 0, j = 0; j < dj; j++) {
2734 for (i = 0; i < di; i++, ij++) {
2735 tdm[ij] = dmji[j*di+i] - dmij[i*dj+j];
2736 }
2737 }
2738
2739 for (icomp = 0; icomp < ncomp; icomp++) {
2740 for (l = 0; l < dl; l++) {
2741 for (k = 0; k < dk; k++) {
2742 tmp = 0;
2743 for (ij = 0; ij < dij; ij++) {
2744 tmp += eri[ij] * tdm[ij];
2745 }
2746 vkl[k*dl+l] += tmp;
2747 vlk[l*dk+k] -= tmp;
2748 eri += dij;
2749 } }
2750 vkl += dkl;
2751 vlk += dkl;
2752 }
2753 }
2754 }
2755 ADD_JKOP(nraa4_ji_s1kl, J, I, K, L, s4);
2756
nraa4_ji_s2kl(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2757 static void nraa4_ji_s2kl(double *eri, double *dm, JKArray *out, int *shls,
2758 int i0, int i1, int j0, int j1,
2759 int k0, int k1, int l0, int l1)
2760 {
2761 nra2ij_ji_s2kl(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2762 }
2763 ADD_JKOP(nraa4_ji_s2kl, J, I, K, L, s4);
2764
nraa4_lk_s1ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2765 static void nraa4_lk_s1ij(double *eri, double *dm, JKArray *out, int *shls,
2766 int i0, int i1, int j0, int j1,
2767 int k0, int k1, int l0, int l1)
2768 {
2769 if (i0 == j0) {
2770 nra2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2771 } else if (k0 == l0) {
2772 nra2ij_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2773 } else {
2774 DECLARE(vij, i, j);
2775 LOCATE(vji, j, i);
2776 DEF_DM(k, l);
2777 DEF_DM(l, k);
2778 int i, j, k, l, ij, icomp;
2779 double buf[MAXCGTO*MAXCGTO];
2780 double tdm;
2781
2782 for (icomp = 0; icomp < ncomp; icomp++) {
2783
2784 for (i = 0; i < dij; i++) { buf[i] = 0; }
2785 for (l = 0; l < dl; l++) {
2786 for (k = 0; k < dk; k++) {
2787 tdm = dmlk[l*dk+k] - dmkl[k*dl+l];
2788 for (ij = 0; ij < dij; ij++) {
2789 buf[ij] += eri[ij] * tdm;
2790 }
2791 eri += dij;
2792 } }
2793
2794 for (ij = 0, j = 0; j < dj; j++) {
2795 for (i = 0; i < di; i++, ij++) {
2796 vij[i*dj+j] += buf[ij];
2797 vji[ij] -= buf[ij];
2798 } }
2799 vij += dij;
2800 vji += dij;
2801 }
2802 }
2803 }
2804 ADD_JKOP(nraa4_lk_s1ij, L, K, I, J, s4);
2805
nraa4_lk_s2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2806 static void nraa4_lk_s2ij(double *eri, double *dm, JKArray *out, int *shls,
2807 int i0, int i1, int j0, int j1,
2808 int k0, int k1, int l0, int l1)
2809 {
2810 nra2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2811 }
2812 ADD_JKOP(nraa4_lk_s2ij, L, K, I, J, s4);
2813
nraa4_lk_a2ij(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2814 static void nraa4_lk_a2ij(double *eri, double *dm, JKArray *out, int *shls,
2815 int i0, int i1, int j0, int j1,
2816 int k0, int k1, int l0, int l1)
2817 {
2818 nra2kl_lk_s1ij(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2819 }
2820 ADD_JKOP(nraa4_lk_a2ij, L, K, I, J, s4);
2821
nraa4_jk_s1il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2822 static void nraa4_jk_s1il(double *eri, double *dm, JKArray *out, int *shls,
2823 int i0, int i1, int j0, int j1,
2824 int k0, int k1, int l0, int l1)
2825 {
2826 if (i0 == j0) {
2827 nra2kl_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2828 } else if (k0 == l0) {
2829 nra2ij_jk_s1il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2830 } else {
2831 DECLARE(vik, i, k);
2832 LOCATE(vil, i, l);
2833 LOCATE(vjk, j, k);
2834 LOCATE(vjl, j, l);
2835 DEF_DM(i, l);
2836 DEF_DM(i, k);
2837 DEF_DM(j, l);
2838 DEF_DM(j, k);
2839 int i, j, k, l, ijkl, icomp;
2840
2841 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2842 for (l = 0; l < dl; l++) {
2843 for (k = 0; k < dk; k++) {
2844 for (j = 0; j < dj; j++) {
2845 for (i = 0; i < di; i++, ijkl++) {
2846 vjk[j*dk+k] += eri[ijkl] * dmil[i*dl+l];
2847 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
2848 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2849 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2850 } } } }
2851 vjk += djk;
2852 vjl += djl;
2853 vik += dik;
2854 vil += dil;
2855 }
2856 }
2857 }
2858 ADD_JKOP(nraa4_jk_s1il, J, K, I, L, s4);
2859
nraa4_jk_s2il(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2860 static void nraa4_jk_s2il(double *eri, double *dm, JKArray *out, int *shls,
2861 int i0, int i1, int j0, int j1,
2862 int k0, int k1, int l0, int l1)
2863 {
2864 if (i0 == j0) {
2865 nra2kl_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2866 } else if (k0 == l0) {
2867 nra2ij_jk_s2il(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2868 } else if (i0 < l0) {
2869 } else if (i0 < k0) {
2870 if (j0 < l0) { // j < l <= i < k
2871 DECLARE(v, i, l);
2872 DEF_DM(j, k);
2873 int i, j, k, l, ijkl, icomp;
2874 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2875 for (l = 0; l < dl; l++) {
2876 for (k = 0; k < dk; k++) {
2877 for (j = 0; j < dj; j++) {
2878 for (i = 0; i < di; i++, ijkl++) {
2879 v[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2880 } } } }
2881 v += dil;
2882 }
2883 } else { // l <= j < i < k
2884 DECLARE(vil, i, l);
2885 LOCATE(vjl, j, l);
2886 DEF_DM(i, k);
2887 DEF_DM(j, k);
2888 int i, j, k, l, ijkl, icomp;
2889 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2890 for (l = 0; l < dl; l++) {
2891 for (k = 0; k < dk; k++) {
2892 for (j = 0; j < dj; j++) {
2893 for (i = 0; i < di; i++, ijkl++) {
2894 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
2895 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2896 } } } }
2897 vjl += djl;
2898 vil += dil;
2899 }
2900 }
2901 } else if (j0 < l0) { // j < l < k <= i
2902 DECLARE(vil, i, l);
2903 LOCATE(vik, i, k);
2904 DEF_DM(j, k);
2905 DEF_DM(j, l);
2906 int i, j, k, l, ijkl, icomp;
2907 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2908 for (l = 0; l < dl; l++) {
2909 for (k = 0; k < dk; k++) {
2910 for (j = 0; j < dj; j++) {
2911 for (i = 0; i < di; i++, ijkl++) {
2912 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2913 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2914 } } } }
2915 vil += dil;
2916 vik += dik;
2917 }
2918 } else if (j0 < k0) { // l <= j < k <= i
2919 DECLARE(vjl, j, l);
2920 LOCATE(vik, i, k);
2921 LOCATE(vil, i, l);
2922 DEF_DM(i, k);
2923 DEF_DM(j, k);
2924 DEF_DM(j, l);
2925 int i, j, k, l, ijkl, icomp;
2926 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2927 for (l = 0; l < dl; l++) {
2928 for (k = 0; k < dk; k++) {
2929 for (j = 0; j < dj; j++) {
2930 for (i = 0; i < di; i++, ijkl++) {
2931 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
2932 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2933 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2934 } } } }
2935 vjl += djl;
2936 vil += dil;
2937 vik += dik;
2938 }
2939 } else { // l < k <= j < i
2940 DECLARE(vjl, j, l);
2941 LOCATE(vik, i, k);
2942 LOCATE(vjk, j, k);
2943 LOCATE(vil, i, l);
2944 DEF_DM(i, l);
2945 DEF_DM(i, k);
2946 DEF_DM(j, l);
2947 DEF_DM(j, k);
2948 int i, j, k, l, ijkl, icomp;
2949 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2950 for (l = 0; l < dl; l++) {
2951 for (k = 0; k < dk; k++) {
2952 for (j = 0; j < dj; j++) {
2953 for (i = 0; i < di; i++, ijkl++) {
2954 vjk[j*dk+k] += eri[ijkl] * dmil[i*dl+l];
2955 vjl[j*dl+l] -= eri[ijkl] * dmik[i*dk+k];
2956 vik[i*dk+k] -= eri[ijkl] * dmjl[j*dl+l];
2957 vil[i*dl+l] += eri[ijkl] * dmjk[j*dk+k];
2958 } } } }
2959 vjk += djk;
2960 vjl += djl;
2961 vik += dik;
2962 vil += dil;
2963 }
2964 }
2965 }
2966 ADD_JKOP(nraa4_jk_s2il, J, K, I, L, s4);
2967
nraa4_li_s1kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)2968 static void nraa4_li_s1kj(double *eri, double *dm, JKArray *out, int *shls,
2969 int i0, int i1, int j0, int j1,
2970 int k0, int k1, int l0, int l1)
2971 {
2972 if (i0 == j0) {
2973 nra2kl_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2974 } else if (k0 == l0) {
2975 nra2ij_li_s1kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
2976 } else {
2977 DECLARE(vki, k, i);
2978 LOCATE(vkj, k, j);
2979 LOCATE(vli, l, i);
2980 LOCATE(vlj, l, j);
2981 DEF_DM(l, i);
2982 DEF_DM(l, j);
2983 DEF_DM(k, i);
2984 DEF_DM(k, j);
2985 int i, j, k, l, ijkl, icomp;
2986 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
2987 for (l = 0; l < dl; l++) {
2988 for (k = 0; k < dk; k++) {
2989 for (j = 0; j < dj; j++) {
2990 for (i = 0; i < di; i++, ijkl++) {
2991 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
2992 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
2993 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
2994 vli[l*di+i] += eri[ijkl] * dmkj[k*dj+j];
2995 } } } }
2996 vkj += dkj;
2997 vki += dki;
2998 vlj += dlj;
2999 vli += dli;
3000 }
3001 }
3002 }
3003 ADD_JKOP(nraa4_li_s1kj, L, I, K, J, s4);
3004
nraa4_li_s2kj(double * eri,double * dm,JKArray * out,int * shls,int i0,int i1,int j0,int j1,int k0,int k1,int l0,int l1)3005 static void nraa4_li_s2kj(double *eri, double *dm, JKArray *out, int *shls,
3006 int i0, int i1, int j0, int j1,
3007 int k0, int k1, int l0, int l1)
3008 {
3009 if (i0 == j0) {
3010 nra2kl_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
3011 } else if (k0 == l0) {
3012 nra2ij_li_s2kj(eri, dm, out, shls, i0, i1, j0, j1, k0, k1, l0, l1);
3013 } else if (k0 < j0) {
3014 } else if (k0 < i0) {
3015 if (l0 < j0) { // l < j < k < i
3016 DECLARE(v, k, j);
3017 DEF_DM(l, i);
3018 int i, j, k, l, ijkl, icomp;
3019 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
3020 for (l = 0; l < dl; l++) {
3021 for (k = 0; k < dk; k++) {
3022 for (j = 0; j < dj; j++) {
3023 for (i = 0; i < di; i++, ijkl++) {
3024 v[k*dj+j] += eri[ijkl] * dmli[l*di+i];
3025 } } } }
3026 v += dkj;
3027 }
3028 } else { // j <= l < k < i
3029 DECLARE(vkj, k, j);
3030 LOCATE(vlj, l, j);
3031 DEF_DM(l, i);
3032 DEF_DM(k, i);
3033 int i, j, k, l, ijkl, icomp;
3034 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
3035 for (l = 0; l < dl; l++) {
3036 for (k = 0; k < dk; k++) {
3037 for (j = 0; j < dj; j++) {
3038 for (i = 0; i < di; i++, ijkl++) {
3039 vkj[k*dj+j] += eri[ijkl] *dmli[l*di+i];
3040 vlj[l*dj+j] -= eri[ijkl] *dmki[k*di+i];
3041 } } } }
3042 vkj += dkj;
3043 vlj += dlj;
3044 }
3045 }
3046 } else if (l0 < j0) { // l < j < i <= k
3047 DECLARE(vki, k, i);
3048 LOCATE(vkj, k, j);
3049 DEF_DM(l, i);
3050 DEF_DM(l, j);
3051 int i, j, k, l, ijkl, icomp;
3052 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
3053 for (l = 0; l < dl; l++) {
3054 for (k = 0; k < dk; k++) {
3055 for (j = 0; j < dj; j++) {
3056 for (i = 0; i < di; i++, ijkl++) {
3057 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
3058 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
3059 } } } }
3060 vkj += dkj;
3061 vki += dki;
3062 }
3063 } else if (l0 < i0) { // j <= l < i <= k
3064 DECLARE(vki, k, i);
3065 LOCATE(vkj, k, j);
3066 LOCATE(vlj, l, j);
3067 DEF_DM(l, i);
3068 DEF_DM(l, j);
3069 DEF_DM(k, i);
3070 int i, j, k, l, ijkl, icomp;
3071 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
3072 for (l = 0; l < dl; l++) {
3073 for (k = 0; k < dk; k++) {
3074 for (j = 0; j < dj; j++) {
3075 for (i = 0; i < di; i++, ijkl++) {
3076 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
3077 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
3078 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
3079 } } } }
3080 vkj += dkj;
3081 vki += dki;
3082 vlj += dlj;
3083 }
3084 } else { // j < i <= l < k
3085 DECLARE(vki, k, i);
3086 LOCATE(vkj, k, j);
3087 LOCATE(vli, l, i);
3088 LOCATE(vlj, l, j);
3089 DEF_DM(l, i);
3090 DEF_DM(l, j);
3091 DEF_DM(k, i);
3092 DEF_DM(k, j);
3093 int i, j, k, l, ijkl, icomp;
3094 for (ijkl = 0, icomp = 0; icomp < ncomp; icomp++) {
3095 for (l = 0; l < dl; l++) {
3096 for (k = 0; k < dk; k++) {
3097 for (j = 0; j < dj; j++) {
3098 for (i = 0; i < di; i++, ijkl++) {
3099 vkj[k*dj+j] += eri[ijkl] * dmli[l*di+i];
3100 vki[k*di+i] -= eri[ijkl] * dmlj[l*dj+j];
3101 vlj[l*dj+j] -= eri[ijkl] * dmki[k*di+i];
3102 vli[l*di+i] += eri[ijkl] * dmkj[k*dj+j];
3103 } } } }
3104 vkj += dkj;
3105 vki += dki;
3106 vlj += dlj;
3107 vli += dli;
3108 }
3109 }
3110 }
3111 ADD_JKOP(nraa4_li_s2kj, L, I, K, J, s4);
3112
3113