1 /*
2 * Author: Qiming Sun <osirpt.sun@gmail.com>
3 *
4 * Breit = Gaunt + gauge
5 * Gaunt ~ - \sigma1\dot\sigma2/r12
6 * gauge ~ 1/2 \sigma1\dot\sigma2/r12 - 1/2 (\sigma1\dot r12) (\sigma2\dot r12)/r12^3
7 * Breit ~ -1/2 \sigma1\dot\sigma2/r12 - 1/2 (\sigma1\dot r12) (\sigma2\dot r12)/r12^3
8 */
9
10 #include <stdlib.h>
11 #include <complex.h>
12 #include "cint_bas.h"
13 #include "cart2sph.h"
14 #include "g1e.h"
15 #include "g2e.h"
16 #include "optimizer.h"
17 #include "cint1e.h"
18 #include "cint2e.h"
19 #include "misc.h"
20 #include "c2f.h"
21
22 #define DECLARE(X) FINT X(double complex *out, FINT *dims, FINT *shls, \
23 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, \
24 CINTOpt *opt, double *cache)
25
26 #define BREIT0(X, ncomp_tensor) \
27 DECLARE(int2e_##X##_spinor); \
28 DECLARE(int2e_gauge_r1_##X##_spinor); \
29 DECLARE(int2e_gauge_r2_##X##_spinor); \
30 void int2e_breit_##X##_optimizer(CINTOpt **opt, FINT *atm, FINT natm, \
31 FINT *bas, FINT nbas, double *env) \
32 { \
33 *opt = NULL; \
34 } \
35 CACHE_SIZE_T int2e_breit_##X##_spinor(double complex *out, FINT *dims, FINT *shls, \
36 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, \
37 CINTOpt *opt, double *cache) \
38 { \
39 return _int2e_breit_drv(out, dims, shls, atm, natm, bas, nbas, env, opt, cache, \
40 ncomp_tensor, &int2e_##X##_spinor, \
41 &int2e_gauge_r1_##X##_spinor, &int2e_gauge_r2_##X##_spinor); \
42 } \
43 FINT cint2e_breit_##X##_spinor(double complex *out, FINT *shls, \
44 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, \
45 CINTOpt *opt) \
46 { \
47 return int2e_breit_##X##_spinor(out, NULL, shls, atm, natm, bas, nbas, env, opt, NULL); \
48 }
49
_copy_to_out(double complex * out,double complex * in,FINT * dims,FINT * counts)50 static void _copy_to_out(double complex *out, double complex *in, FINT *dims, FINT *counts)
51 {
52 if (out == in) {
53 return;
54 }
55 FINT ni = dims[0];
56 FINT nj = dims[1];
57 FINT nk = dims[2];
58 FINT nij = ni * nj;
59 FINT nijk = nij * nk;
60 FINT di = counts[0];
61 FINT dj = counts[1];
62 FINT dk = counts[2];
63 FINT dl = counts[3];
64 FINT dij = di * dj;
65 FINT dijk = dij * dk;
66 FINT i, j, k, l;
67 double complex *pin, *pout;
68 for (l = 0; l < dl; l++) {
69 for (k = 0; k < dk; k++) {
70 pin = in + k * dij;
71 pout = out + k * nij;
72 for (j = 0; j < dj; j++) {
73 for (i = 0; i < di; i++) {
74 pout[j*ni+i] = pin[j*di+i];
75 } }
76 }
77 in += dijk;
78 out += nijk;
79 }
80 }
81
_int2e_breit_drv(double complex * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache,FINT ncomp_tensor,FINT (* f_gaunt)(),FINT (* f_gauge_r1)(),FINT (* f_gauge_r2)())82 static CACHE_SIZE_T _int2e_breit_drv(double complex *out, FINT *dims, FINT *shls,
83 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env,
84 CINTOpt *opt, double *cache, FINT ncomp_tensor,
85 FINT (*f_gaunt)(), FINT (*f_gauge_r1)(), FINT (*f_gauge_r2)())
86 {
87 if (out == NULL) {
88 CACHE_SIZE_T cache_size1 = (*f_gauge_r1)(NULL, NULL, shls,
89 atm, natm, bas, nbas, env, NULL, cache);
90 CACHE_SIZE_T cache_size2 = (*f_gauge_r2)(NULL, NULL, shls,
91 atm, natm, bas, nbas, env, NULL, cache);
92 return MAX(cache_size1, cache_size2);
93 }
94
95 FINT counts[4];
96 counts[0] = CINTcgto_spinor(shls[0], bas);
97 counts[1] = CINTcgto_spinor(shls[1], bas);
98 counts[2] = CINTcgto_spinor(shls[2], bas);
99 counts[3] = CINTcgto_spinor(shls[3], bas);
100 FINT nop = counts[0] * counts[1] * counts[2] * counts[3] * ncomp_tensor;
101 double complex *buf = malloc(sizeof(double complex) * nop*2);
102 double complex *buf1;
103 if (dims == NULL) {
104 dims = counts;
105 buf1 = out;
106 } else {
107 buf1 = buf + nop;
108 }
109
110 FINT has_value = (*f_gaunt)(buf1, NULL, shls, atm, natm, bas, nbas, env, NULL, cache);
111
112 FINT i;
113 has_value = ((*f_gauge_r1)(buf, NULL, shls, atm, natm, bas, nbas, env, NULL, cache) ||
114 has_value);
115 /* [1/2 gaunt] - [1/2 xxx*\sigma1\dot r1] */
116 if (has_value) {
117 for (i = 0; i < nop; i++) {
118 buf1[i] = -buf1[i] - buf[i];
119 }
120 }
121 /* ... [- 1/2 xxx*\sigma1\dot(-r2)] */
122 has_value = ((*f_gauge_r2)(buf, NULL, shls, atm, natm, bas, nbas, env, NULL, cache) ||
123 has_value);
124 if (has_value) {
125 for (i = 0; i < nop; i++) {
126 buf1[i] = (buf1[i] + buf[i]) * .5;
127 }
128 }
129 _copy_to_out(out, buf1, dims, counts);
130 free(buf);
131 return has_value;
132 }
133
134
135 BREIT0(ssp1ssp2, 1);
136 BREIT0(ssp1sps2, 1);
137 BREIT0(sps1ssp2, 1);
138 BREIT0(sps1sps2, 1);
139
140 /* based on
141 * '("int2e_breit_r1p2" ( nabla \, r0 \| dot nabla-r12 \| \, nabla ))
142 */
CINTgout2e_int2e_breit_r1p2(double * gout,double * g,FINT * idx,CINTEnvVars * envs,FINT gout_empty)143 static void CINTgout2e_int2e_breit_r1p2(double *gout, double *g,
144 FINT *idx, CINTEnvVars *envs, FINT gout_empty) {
145 FINT nf = envs->nf;
146 FINT nrys_roots = envs->nrys_roots;
147 FINT ix, iy, iz, i, n;
148 double *g0 = g;
149 double *g1 = g0 + envs->g_size * 3;
150 double *g2 = g1 + envs->g_size * 3;
151 double *g3 = g2 + envs->g_size * 3;
152 double *g4 = g3 + envs->g_size * 3;
153 double *g5 = g4 + envs->g_size * 3;
154 double *g6 = g5 + envs->g_size * 3;
155 double *g7 = g6 + envs->g_size * 3;
156 double *g8 = g7 + envs->g_size * 3;
157 double *g9 = g8 + envs->g_size * 3;
158 double *g10 = g9 + envs->g_size * 3;
159 double *g11 = g10 + envs->g_size * 3;
160 double *g12 = g11 + envs->g_size * 3;
161 double *g13 = g12 + envs->g_size * 3;
162 double *g14 = g13 + envs->g_size * 3;
163 double *g15 = g14 + envs->g_size * 3;
164 G2E_D_L(g1, g0, envs->i_l+2, envs->j_l+2, envs->k_l+0, envs->l_l+0);
165 G2E_R0J(g3, g1, envs->i_l+1, envs->j_l+0, envs->k_l, envs->l_l);
166 G2E_D_J(g4, g0, envs->i_l+1, envs->j_l+1, envs->k_l, envs->l_l);
167 G2E_D_I(g5, g0, envs->i_l+1, envs->j_l+1, envs->k_l, envs->l_l);
168 for (ix = 0; ix < envs->g_size * 3; ix++) {g4[ix] += g5[ix];}
169 G2E_D_J(g5, g1, envs->i_l+1, envs->j_l+1, envs->k_l, envs->l_l);
170 G2E_D_I(g6, g1, envs->i_l+1, envs->j_l+1, envs->k_l, envs->l_l);
171 for (ix = 0; ix < envs->g_size * 3; ix++) {g5[ix] += g6[ix];}
172 G2E_R0J(g7, g5, envs->i_l+1, envs->j_l+0, envs->k_l, envs->l_l);
173 G2E_D_I(g12, g4, envs->i_l+0, envs->j_l, envs->k_l, envs->l_l);
174 G2E_D_I(g15, g7, envs->i_l+0, envs->j_l, envs->k_l, envs->l_l);
175 double s;
176 for (n = 0; n < nf; n++) {
177 ix = idx[0+n*3];
178 iy = idx[1+n*3];
179 iz = idx[2+n*3];
180 s = 0;
181 for (i = 0; i < nrys_roots; i++) {
182 s += g15[ix+i] * g0[iy+i] * g0[iz+i];
183 s += g12[ix+i] * g3[iy+i] * g0[iz+i];
184 s += g12[ix+i] * g0[iy+i] * g3[iz+i];
185 s += g3[ix+i] * g12[iy+i] * g0[iz+i];
186 s += g0[ix+i] * g15[iy+i] * g0[iz+i];
187 s += g0[ix+i] * g12[iy+i] * g3[iz+i];
188 s += g3[ix+i] * g0[iy+i] * g12[iz+i];
189 s += g0[ix+i] * g3[iy+i] * g12[iz+i];
190 s += g0[ix+i] * g0[iy+i] * g15[iz+i];
191 }
192 if (gout_empty) {
193 gout[n] = s;
194 } else {
195 gout[n] += s;
196 }
197 }
198 }
int2e_breit_r1p2_optimizer(CINTOpt ** opt,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env)199 void int2e_breit_r1p2_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {
200 FINT ng[] = {2, 2, 0, 1, 4, 1, 1, 1};
201 CINTall_2e_optimizer(opt, ng, atm, natm, bas, nbas, env);
202 }
int2e_breit_r1p2_cart(double * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)203 CACHE_SIZE_T int2e_breit_r1p2_cart(double *out, FINT *dims, FINT *shls,
204 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
205 FINT ng[] = {2, 2, 0, 1, 4, 1, 1, 1};
206 CINTEnvVars envs;
207 CINTinit_int2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
208 envs.f_gout = &CINTgout2e_int2e_breit_r1p2;
209 return CINT2e_drv(out, dims, &envs, opt, cache, &c2s_cart_2e1);
210 } // int2e_breit_r1p2_cart
int2e_breit_r1p2_sph(double * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)211 CACHE_SIZE_T int2e_breit_r1p2_sph(double *out, FINT *dims, FINT *shls,
212 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
213 FINT ng[] = {2, 2, 0, 1, 4, 1, 1, 1};
214 CINTEnvVars envs;
215 CINTinit_int2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
216 envs.f_gout = &CINTgout2e_int2e_breit_r1p2;
217 return CINT2e_drv(out, dims, &envs, opt, cache, &c2s_sph_2e1);
218 } // int2e_breit_r1p2_sph
int2e_breit_r1p2_spinor(double complex * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)219 CACHE_SIZE_T int2e_breit_r1p2_spinor(double complex *out, FINT *dims, FINT *shls,
220 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
221 FINT ng[] = {2, 2, 0, 1, 4, 1, 1, 1};
222 CINTEnvVars envs;
223 CINTinit_int2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
224 envs.f_gout = &CINTgout2e_int2e_breit_r1p2;
225 return CINT2e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_2e1i, &c2s_sf_2e2i);
226 } // int2e_breit_r1p2_spinor
227 ALL_CINT(int2e_breit_r1p2)
ALL_CINT_FORTRAN_(int2e_breit_r1p2)228 ALL_CINT_FORTRAN_(int2e_breit_r1p2)
229
230 /* based on
231 * '("int2e_breit_r2p2" ( nabla \, r0 \| dot nabla-r12 \| \, nabla ))
232 */
233 static void CINTgout2e_int2e_breit_r2p2(double *gout,
234 double *g, FINT *idx, CINTEnvVars *envs, FINT gout_empty) {
235 FINT nf = envs->nf;
236 FINT nrys_roots = envs->nrys_roots;
237 FINT ix, iy, iz, i, n;
238 double *g0 = g;
239 double *g1 = g0 + envs->g_size * 3;
240 double *g2 = g1 + envs->g_size * 3;
241 double *g3 = g2 + envs->g_size * 3;
242 double *g4 = g3 + envs->g_size * 3;
243 double *g5 = g4 + envs->g_size * 3;
244 double *g6 = g5 + envs->g_size * 3;
245 double *g7 = g6 + envs->g_size * 3;
246 double *g8 = g7 + envs->g_size * 3;
247 double *g9 = g8 + envs->g_size * 3;
248 double *g10 = g9 + envs->g_size * 3;
249 double *g11 = g10 + envs->g_size * 3;
250 double *g12 = g11 + envs->g_size * 3;
251 double *g13 = g12 + envs->g_size * 3;
252 double *g14 = g13 + envs->g_size * 3;
253 double *g15 = g14 + envs->g_size * 3;
254 G2E_R0L(g2, g0, envs->i_l+2, envs->j_l+1, envs->k_l+0, envs->l_l+1);
255 G2E_D_L(g3, g2, envs->i_l+2, envs->j_l+1, envs->k_l+0, envs->l_l+0);
256 G2E_D_J(g4, g0, envs->i_l+1, envs->j_l+0, envs->k_l, envs->l_l);
257 G2E_D_I(g5, g0, envs->i_l+1, envs->j_l+0, envs->k_l, envs->l_l);
258 for (ix = 0; ix < envs->g_size * 3; ix++) {g4[ix] += g5[ix];}
259 G2E_D_J(g7, g3, envs->i_l+1, envs->j_l+0, envs->k_l, envs->l_l);
260 G2E_D_I(g8, g3, envs->i_l+1, envs->j_l+0, envs->k_l, envs->l_l);
261 for (ix = 0; ix < envs->g_size * 3; ix++) {g7[ix] += g8[ix];}
262 G2E_D_I(g12, g4, envs->i_l+0, envs->j_l, envs->k_l, envs->l_l);
263 G2E_D_I(g15, g7, envs->i_l+0, envs->j_l, envs->k_l, envs->l_l);
264 double s;
265 for (n = 0; n < nf; n++) {
266 ix = idx[0+n*3];
267 iy = idx[1+n*3];
268 iz = idx[2+n*3];
269 s = 0;
270 for (i = 0; i < nrys_roots; i++) {
271 s += g15[ix+i] * g0[iy+i] * g0[iz+i];
272 s += g12[ix+i] * g3[iy+i] * g0[iz+i];
273 s += g12[ix+i] * g0[iy+i] * g3[iz+i];
274 s += g3[ix+i] * g12[iy+i] * g0[iz+i];
275 s += g0[ix+i] * g15[iy+i] * g0[iz+i];
276 s += g0[ix+i] * g12[iy+i] * g3[iz+i];
277 s += g3[ix+i] * g0[iy+i] * g12[iz+i];
278 s += g0[ix+i] * g3[iy+i] * g12[iz+i];
279 s += g0[ix+i] * g0[iy+i] * g15[iz+i];
280 }
281 if (gout_empty) {
282 gout[n] = s;
283 } else {
284 gout[n] += s;
285 }
286 }
287 }
int2e_breit_r2p2_optimizer(CINTOpt ** opt,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env)288 void int2e_breit_r2p2_optimizer(CINTOpt **opt, FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env) {
289 FINT ng[] = {2, 1, 0, 2, 4, 1, 1, 1};
290 CINTall_2e_optimizer(opt, ng, atm, natm, bas, nbas, env);
291 }
int2e_breit_r2p2_cart(double * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)292 CACHE_SIZE_T int2e_breit_r2p2_cart(double *out, FINT *dims, FINT *shls,
293 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
294 FINT ng[] = {2, 1, 0, 2, 4, 1, 1, 1};
295 CINTEnvVars envs;
296 CINTinit_int2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
297 envs.f_gout = &CINTgout2e_int2e_breit_r2p2;
298 return CINT2e_drv(out, dims, &envs, opt, cache, &c2s_cart_2e1);
299 } // int2e_breit_r2p2_cart
int2e_breit_r2p2_sph(double * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)300 CACHE_SIZE_T int2e_breit_r2p2_sph(double *out, FINT *dims, FINT *shls,
301 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
302 FINT ng[] = {2, 1, 0, 2, 4, 1, 1, 1};
303 CINTEnvVars envs;
304 CINTinit_int2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
305 envs.f_gout = &CINTgout2e_int2e_breit_r2p2;
306 return CINT2e_drv(out, dims, &envs, opt, cache, &c2s_sph_2e1);
307 } // int2e_breit_r2p2_sph
int2e_breit_r2p2_spinor(double complex * out,FINT * dims,FINT * shls,FINT * atm,FINT natm,FINT * bas,FINT nbas,double * env,CINTOpt * opt,double * cache)308 CACHE_SIZE_T int2e_breit_r2p2_spinor(double complex *out, FINT *dims, FINT *shls,
309 FINT *atm, FINT natm, FINT *bas, FINT nbas, double *env, CINTOpt *opt, double *cache) {
310 FINT ng[] = {2, 1, 0, 2, 4, 1, 1, 1};
311 CINTEnvVars envs;
312 CINTinit_int2e_EnvVars(&envs, ng, shls, atm, natm, bas, nbas, env);
313 envs.f_gout = &CINTgout2e_int2e_breit_r2p2;
314 return CINT2e_spinor_drv(out, dims, &envs, opt, cache, &c2s_sf_2e1i, &c2s_sf_2e2i);
315 } // int2e_breit_r2p2_spinor
316 ALL_CINT(int2e_breit_r2p2)
317 ALL_CINT_FORTRAN_(int2e_breit_r2p2)
318