MPQC 3.0.0-alpha
Loading...
Searching...
No Matches
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
35namespace 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 --------------------------------------------------------------------------------*/
45inline 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 Wed Sep 25 2024 02:45:29 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.12.0.