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