MPQC  3.0.0-alpha
g12nc_quartet_data.h
1 //
2 // g12nc_quartet_data.h
3 //
4 // Copyright (C) 2005 Edward Valeev
5 //
6 // Author: Edward Valeev <evaleev@vt.edu>
7 // Maintainer: EV
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27 
28 #include <util/misc/math.h>
29 #include <chemistry/qc/libint2/static.h>
30 #include <chemistry/qc/libint2/libint2_utils.h>
31 
32 #ifndef _chemistry_qc_libint2_g12ncquartetdata_h
33 #define _chemistry_qc_libint2_g12ncquartetdata_h
34 
35 namespace sc {
36 
37 /*--------------------------------------------------------------------------------
38  This function computes constants used in OSRR for a given quartet of primitives
39 
40  gamma is the exponents of the Gaussian geminal in the integral
41 
42  r12_2_g12_pfac is the prefactor by which the r12^2 * g12 integrals are scaled
43  it's necessary to produce [g12,[T1,g12]]
44  --------------------------------------------------------------------------------*/
45 inline void G12NCLibint2::g12nc_quartet_data_(prim_data *Data, double scale, OperType otype,
46  const ContractedGeminal* gbra, const ContractedGeminal* gket)
47 {
48 #define STATIC_OO2NP1
49 #include "static.h"
50 
51  /*----------------
52  Local variables
53  ----------------*/
54  double P[3], Q[3], PQ[3], W[3];
55  const double small_T = 1E-15; /*--- Use only one term in Taylor expansion of Fj(T) if T < small_T ---*/
56 
57  const int p1 = quartet_info_.p1;
58  const int p2 = quartet_info_.p2;
59  const int p3 = quartet_info_.p3;
60  const int p4 = quartet_info_.p4;
61 
62  const double a1 = int_shell1_->exponent(quartet_info_.p1);
63  const double a2 = int_shell2_->exponent(quartet_info_.p2);
64  const double a3 = int_shell3_->exponent(quartet_info_.p3);
65  const double a4 = int_shell4_->exponent(quartet_info_.p4);
66 
67  prim_pair_t* pair12;
68  prim_pair_t* pair34;
69  if (!quartet_info_.p13p24) {
70  pair12 = quartet_info_.shell_pair12->prim_pair(*quartet_info_.op1,*quartet_info_.op2);
71  pair34 = quartet_info_.shell_pair34->prim_pair(*quartet_info_.op3,*quartet_info_.op4);
72  }
73  else {
74  pair12 = quartet_info_.shell_pair34->prim_pair(*quartet_info_.op3,*quartet_info_.op4);
75  pair34 = quartet_info_.shell_pair12->prim_pair(*quartet_info_.op1,*quartet_info_.op2);
76  }
77 
78  //
79  // prefactors for OSRR do not depend on the integral kernel
80  //
81  const double zeta = pair12->gamma;
82  const double eta = pair34->gamma;
83  const double ooz = 1.0/zeta;
84  const double ooe = 1.0/eta;
85  const double ooze = 1.0/(zeta+eta);
86 #if LIBINT2_DEFINED(eri,roz)
87  Data->roz[0] = eta*ooze;
88  double rho = zeta*Data->roz[0];
89 #else
90  double rho = zeta * eta * ooze;
91 #endif
92 
93  const double pfac_norm = int_shell1_->coefficient_unnorm(quartet_info_.gc1,p1)*
94  int_shell2_->coefficient_unnorm(quartet_info_.gc2,p2)*
95  int_shell3_->coefficient_unnorm(quartet_info_.gc3,p3)*
96  int_shell4_->coefficient_unnorm(quartet_info_.gc4,p4);
97  const double pfac_normovlp = pfac_norm * pair12->ovlp * pair34->ovlp * scale;
98 
99 
100  P[0] = pair12->P[0];
101  P[1] = pair12->P[1];
102  P[2] = pair12->P[2];
103  Q[0] = pair34->P[0];
104  Q[1] = pair34->P[1];
105  Q[2] = pair34->P[2];
106  PQ[0] = P[0] - Q[0];
107  PQ[1] = P[1] - Q[1];
108  PQ[2] = P[2] - Q[2];
109  const double PQ2 = PQ[0]*PQ[0] + PQ[1]*PQ[1] + PQ[2]*PQ[2];
110  const double T = rho*PQ2;
111 
112  // Coulomb integral
113  if (otype == coulomb) {
114  double pfac = 2.0*sqrt(rho*M_1_PI)*pfac_normovlp;
115  if(T < small_T){
116  assign_FjT(Data,quartet_info_.am,oo2np1,pfac);
117  }
118  else {
119  Fm_Eval_->eval(Fm_table_,T,quartet_info_.am);
120  assign_FjT(Data,quartet_info_.am,Fm_table_,pfac);
121  }
122  }
123  else {
124 
125  // this stores (ss|Oper|ss)^(m) integrals
126  double ss_oper_ss[4*LIBINT2_MAX_AM_eri+1];
127  std::fill(ss_oper_ss, ss_oper_ss + quartet_info_.am + 1, 0.0);
128 
129  // f12_coulomb and f12 integrals
130  if (otype == f12 || otype == f12_coulomb) {
131 
132  MPQC_ASSERT(gbra != 0);
133  const size_t ngbra = gbra->size();
134  for(size_t ig=0; ig<ngbra; ++ig) {
135  const PrimitiveGeminal& gbra_i = gbra->operator[](ig);
136  const double gamma = gbra_i.first;
137  const double gcoef = gbra_i.second;
138  const double rhog = rho + gamma;
139  const double oorhog = 1.0/rhog;
140 
141  const double gorg = gamma * oorhog;
142  const double rorg = rho * oorhog;
143  const double sqrt_rorg = sqrt(rorg);
145  const double SS_K0G12_SS = rorg * sqrt_rorg * exp(-gorg*T) * pfac_normovlp;
146 
147  if (otype == f12_coulomb) {
148  const double rorgT = rorg * T;
149  double pfac = 2.0 * sqrt(rhog*M_1_PI) * SS_K0G12_SS;
150 
151  const double *F;
152  if(rorgT < small_T){
153  F = oo2np1;
154  }
155  else {
156  Fm_Eval_->eval(Fm_table_,rorgT,quartet_info_.am);
157  F = Fm_table_;
158  }
159 
160  double g_i[4*LIBINT2_MAX_AM_eri+1];
161  double r_i[4*LIBINT2_MAX_AM_eri+1];
162  double oorhog_i[4*LIBINT2_MAX_AM_eri+1];
163  g_i[0] = 1.0;
164  r_i[0] = 1.0;
165  oorhog_i[0] = 1.0;
166  for(int i=1; i<=quartet_info_.am; i++) {
167  g_i[i] = g_i[i-1] * gamma;
168  r_i[i] = r_i[i-1] * rho;
169  oorhog_i[i] = oorhog_i[i-1] * oorhog;
170  }
171  for(int m=0; m<=quartet_info_.am; m++) {
172  double ssss = 0.0;
173  for(int k=0; k<=m; k++) {
174  ssss += ExpMath_.bc[m][k] * r_i[k] * g_i[m-k] * F[k];
175  }
176  ss_oper_ss[m] += gcoef * pfac * ssss * oorhog_i[m];
177  }
178 
179  }
180 
181  if (otype == f12) {
182 
183  double ss_oper_ss_m = SS_K0G12_SS * gcoef;
184  ss_oper_ss[0] += ss_oper_ss_m;
185  for(int m=1; m<=quartet_info_.am; ++m) {
186  ss_oper_ss_m *= gorg;
187  ss_oper_ss[m] += ss_oper_ss_m;
188  }
189 
190  }
191 
192  } // loop over gaussian geminals
193  } // one correlation factor involved
194 
195  // f12_2 and f12_T1_f12 integrals
196  if (otype == f12_2 || otype == f12_T1_f12) {
197 
198  MPQC_ASSERT(gbra != 0);
199  MPQC_ASSERT(gket != 0);
200  const size_t ngbra = gbra->size();
201  const size_t ngket = gket->size();
202  for(size_t igbra=0; igbra<ngbra; ++igbra) {
203  const PrimitiveGeminal& gbra_i = gbra->operator[](igbra);
204  const double gamma_bra = gbra_i.first;
205  const double gcoef_bra = gbra_i.second;
206  for(size_t igket=0; igket<ngket; ++igket) {
207  const PrimitiveGeminal& gket_i = gket->operator[](igket);
208  const double gamma_ket = gket_i.first;
209  const double gcoef_ket = gket_i.second;
210 
211  const double gamma = gamma_bra + gamma_ket;
212  const double gcoef = gcoef_bra * gcoef_ket;
213 
214  const double rhog = rho + gamma;
215  const double oorhog = 1.0/rhog;
216 
217  const double gorg = gamma * oorhog;
218  const double rorg = rho * oorhog;
219  const double sqrt_rorg = sqrt(rorg);
221  const double SS_K0G12_SS = rorg * sqrt_rorg * exp(-gorg*T) * pfac_normovlp;
222  const double rorgT = rorg * T;
224  const double SS_K2G12_SS_0 = (1.5 + rorgT) * (SS_K0G12_SS * oorhog);
225  const double SS_K2G12_SS_m1 = rorg * (SS_K0G12_SS * oorhog);
226 
227  if (otype == f12_2) {
228 
229  double ss_oper_ss_m = SS_K0G12_SS * gcoef;
230  ss_oper_ss[0] += ss_oper_ss_m;
231  for(int m=1; m<=quartet_info_.am; ++m) {
232  ss_oper_ss_m *= gorg;
233  ss_oper_ss[m] += ss_oper_ss_m;
234  }
235 
236  }
237 
238  if (otype == f12_T1_f12) {
239 
240  double SS_K2G12_SS_gorg_m = SS_K2G12_SS_0 * (gcoef * 4.0 * gamma_bra * gamma_ket);
241  double SS_K2G12_SS_gorg_m1 = SS_K2G12_SS_m1 * (gcoef * 4.0 * gamma_bra * gamma_ket);
242  ss_oper_ss[0] += SS_K2G12_SS_gorg_m;
243  for(int m=1; m<=quartet_info_.am; ++m) {
244  SS_K2G12_SS_gorg_m *= gorg;
245  ss_oper_ss[m] += SS_K2G12_SS_gorg_m - m * SS_K2G12_SS_gorg_m1;
246  SS_K2G12_SS_gorg_m1 *= gorg;
247  }
248 
249  }
250 
251  }
252 
253  } // loop over gaussian geminals
254  } // two correlation factors involved
255 
256  assign_FjT(Data,quartet_info_.am,ss_oper_ss,1.0);
257 
258  }
259 
260  // these prefactors only necessary if angular momenta != 0
261  if (quartet_info_.am != 0) {
262 #if LIBINT2_DEFINED(eri,oo2ze)
263  Data->oo2ze[0] = 0.5 * ooze;
264 #endif
265 #if LIBINT2_DEFINED(eri,roe)
266  Data->roe[0] = zeta * ooze;
267 #endif
268 #if LIBINT2_DEFINED(eri,oo2z)
269  Data->oo2z[0] = 0.5 * ooz;
270 #endif
271 #if LIBINT2_DEFINED(eri,oo2e)
272  Data->oo2e[0] = 0.5 * ooe;
273 #endif
274  W[0] = (zeta * P[0] + eta * Q[0]) * ooze;
275  W[1] = (zeta * P[1] + eta * Q[1]) * ooze;
276  W[2] = (zeta * P[2] + eta * Q[2]) * ooze;
277 
278  /* PA */
279 #if LIBINT2_DEFINED(eri,PA_x)
280  Data->PA_x[0] = P[0] - quartet_info_.A[0];
281 #endif
282 #if LIBINT2_DEFINED(eri,PA_y)
283  Data->PA_y[0] = P[1] - quartet_info_.A[1];
284 #endif
285 #if LIBINT2_DEFINED(eri,PA_z)
286  Data->PA_z[0] = P[2] - quartet_info_.A[2];
287 #endif
288  /* QC */
289 #if LIBINT2_DEFINED(eri,QC_x)
290  Data->QC_x[0] = Q[0] - quartet_info_.C[0];
291 #endif
292 #if LIBINT2_DEFINED(eri,QC_y)
293  Data->QC_y[0] = Q[1] - quartet_info_.C[1];
294 #endif
295 #if LIBINT2_DEFINED(eri,QC_z)
296  Data->QC_z[0] = Q[2] - quartet_info_.C[2];
297 #endif
298  /* WP */
299 #if LIBINT2_DEFINED(eri,WP_x)
300  Data->WP_x[0] = W[0] - P[0];
301 #endif
302 #if LIBINT2_DEFINED(eri,WP_y)
303  Data->WP_y[0] = W[1] - P[1];
304 #endif
305 #if LIBINT2_DEFINED(eri,WP_z)
306  Data->WP_z[0] = W[2] - P[2];
307 #endif
308  /* WQ */
309 #if LIBINT2_DEFINED(eri,WQ_x)
310  Data->WQ_x[0] = W[0] - Q[0];
311 #endif
312 #if LIBINT2_DEFINED(eri,WQ_y)
313  Data->WQ_y[0] = W[1] - Q[1];
314 #endif
315 #if LIBINT2_DEFINED(eri,WQ_z)
316  Data->WQ_z[0] = W[2] - Q[2];
317 #endif
318 
319  // using ITR?
320 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x)
321  Data->TwoPRepITR_pfac0_0_0_x[0] = - (a2*(quartet_info_.A[0]-quartet_info_.B[0]) + a4*(quartet_info_.C[0]-quartet_info_.D[0]))/zeta;
322 #endif
323 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y)
324  Data->TwoPRepITR_pfac0_0_0_y[0] = - (a2*(quartet_info_.A[1]-quartet_info_.B[1]) + a4*(quartet_info_.C[1]-quartet_info_.D[1]))/zeta;
325 #endif
326 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z)
327  Data->TwoPRepITR_pfac0_0_0_z[0] = - (a2*(quartet_info_.A[2]-quartet_info_.B[2]) + a4*(quartet_info_.C[2]-quartet_info_.D[2]))/zeta;
328 #endif
329 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x)
330  Data->TwoPRepITR_pfac0_1_0_x[0] = - (a2*(quartet_info_.A[0]-quartet_info_.B[0]) + a4*(quartet_info_.C[0]-quartet_info_.D[0]))/eta;
331 #endif
332 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y)
333  Data->TwoPRepITR_pfac0_1_0_y[0] = - (a2*(quartet_info_.A[1]-quartet_info_.B[1]) + a4*(quartet_info_.C[1]-quartet_info_.D[1]))/eta;
334 #endif
335 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z)
336  Data->TwoPRepITR_pfac0_1_0_z[0] = - (a2*(quartet_info_.A[2]-quartet_info_.B[2]) + a4*(quartet_info_.C[2]-quartet_info_.D[2]))/eta;
337 #endif
338 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x)
339  Data->TwoPRepITR_pfac0_0_1_x[0] = (a1*(quartet_info_.A[0]-quartet_info_.B[0]) + a3*(quartet_info_.C[0]-quartet_info_.D[0]))/zeta;
340 #endif
341 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y)
342  Data->TwoPRepITR_pfac0_0_1_y[0] = (a1*(quartet_info_.A[1]-quartet_info_.B[1]) + a3*(quartet_info_.C[1]-quartet_info_.D[1]))/zeta;
343 #endif
344 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z)
345  Data->TwoPRepITR_pfac0_0_1_z[0] = (a1*(quartet_info_.A[2]-quartet_info_.B[2]) + a3*(quartet_info_.C[2]-quartet_info_.D[2]))/zeta;
346 #endif
347 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_x)
348  Data->TwoPRepITR_pfac0_1_1_x[0] = (a1*(quartet_info_.A[0]-quartet_info_.B[0]) + a3*(quartet_info_.C[0]-quartet_info_.D[0]))/eta;
349 #endif
350 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y)
351  Data->TwoPRepITR_pfac0_1_1_y[0] = (a1*(quartet_info_.A[1]-quartet_info_.B[1]) + a3*(quartet_info_.C[1]-quartet_info_.D[1]))/eta;
352 #endif
353 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z)
354  Data->TwoPRepITR_pfac0_1_1_z[0] = (a1*(quartet_info_.A[2]-quartet_info_.B[2]) + a3*(quartet_info_.C[2]-quartet_info_.D[2]))/eta;
355 #endif
356 #if LIBINT2_DEFINED(eri,eoz)
357  Data->eoz[0] = eta * ooz;
358 #endif
359 #if LIBINT2_DEFINED(eri,zoe)
360  Data->zoe[0] = zeta * ooe;
361 #endif
362 
363  }
364 
365  return;
366 }
367 
368 }
369 
370 #endif
371 
372 // Local Variables:
373 // mode: c++
374 // c-file-style: "CLJ"
375 // End:
sc::GaussianShell::coefficient_unnorm
double coefficient_unnorm(int con, int prim) const
Returns the contraction coef for unnormalized primitives.
Definition: gaussshell.h:232
sc::GaussianShell::exponent
double exponent(int iprim) const
Returns the exponents of the given primitive.
Definition: gaussshell.h:238
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
Generated at Sun Jan 26 2020 23:23:58 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.