source: ThirdParty/mpqc_open/src/lib/chemistry/qc/cints/hcore.cc

Candidate_v1.6.1
Last change on this file was 860145, checked in by Frederik Heber <heber@…>, 8 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 10.3 KB
Line 
1//
2// hcore.cc
3//
4// Copyright (C) 2001 Edward Valeev
5//
6// Author: Edward Valeev <edward.valeev@chemistry.gatech.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#ifdef __GNUG__
29#pragma implementation
30#endif
31
32#include <util/misc/math.h>
33
34#include <chemistry/qc/cints/int1e.h>
35#include <chemistry/qc/cints/macros.h>
36
37using namespace sc;
38
39void Int1eCints::hcore(int sh1, int sh2)
40{
41 zero_buffers_();
42 compute_doublet_info_(sh1, sh2);
43
44 int maxam1 = int_shell1_->max_am();
45 int minam1 = int_shell1_->min_am();
46 int maxam2 = int_shell2_->max_am();
47 int minam2 = int_shell2_->min_am();
48
49 if (maxam1 != minam1 || maxam2 != minam2) {
50 // fail();
51 hcore_full_general_();
52 }
53 else {
54 hcore_full_general_();
55 }
56}
57
58
59void Int1eCints::hcore_full_general_()
60{
61 int maxam1 = int_shell1_->max_am();
62 int maxam2 = int_shell2_->max_am();
63 int z1weight = 1;
64 int y1weight = maxam1 + 1;
65 int x1weight = y1weight * y1weight;
66 int z2weight = 1;
67 int y2weight = maxam2 + 1;
68 int x2weight = y2weight * y2weight;
69
70 /* See if need to transform to spherical harmonics */
71 bool need_cart2sph_transform = false;
72 if (int_shell1_->has_pure() ||
73 int_shell2_->has_pure())
74 need_cart2sph_transform = true;
75
76 /* See if contraction quartets need to be resorted into a shell quartet */
77 bool need_sort_to_shell_doublet = false;
78 int num_gen_shells = 0;
79 if (int_shell1_->ncontraction() > 1)
80 num_gen_shells++;
81 if (int_shell2_->ncontraction() > 1)
82 num_gen_shells++;
83 if (maxam1 + maxam2 && num_gen_shells >= 1)
84 need_sort_to_shell_doublet = true;
85
86 /* Determine where integrals need to go at each stage */
87 if (need_sort_to_shell_doublet) {
88 prim_ints_ = cart_ints_;
89 if (need_cart2sph_transform)
90 contr_doublets_ = sphharm_ints_;
91 else
92 contr_doublets_ = cart_ints_;
93 shell_doublet_ = target_ints_buffer_;
94 }
95 else {
96 if (need_cart2sph_transform) {
97 prim_ints_ = cart_ints_;
98 contr_doublets_ = target_ints_buffer_;
99 shell_doublet_ = target_ints_buffer_;
100 }
101 else {
102 prim_ints_ = target_ints_buffer_;
103 shell_doublet_ = target_ints_buffer_;
104 }
105 }
106
107 /* Begin loops over primitives. */
108 for (int p1=0; p1<int_shell1_->nprimitive(); p1++) {
109 double a1 = int_shell1_->exponent(p1);
110 for (int p2=0; p2<int_shell2_->nprimitive(); p2++) {
111 double a2 = int_shell2_->exponent(p2);
112
113 double gamma = a1+a2;
114 double oog = 1.0/gamma;
115 double over_pf = exp(-a1*a2*doublet_info_.AB2*oog)*sqrt(M_PI*oog)*M_PI*oog;
116
117 double P[3], PA[3], PB[3], PC[3];
118 for(int xyz=0; xyz<3; xyz++) {
119 P[xyz] = (a1*doublet_info_.A[xyz] + a2*doublet_info_.B[xyz])*oog;
120 PA[xyz] = P[xyz] - doublet_info_.A[xyz];
121 PB[xyz] = P[xyz] - doublet_info_.B[xyz];
122 }
123
124 OI_OSrecurs_(OIX_,OIY_,OIZ_,PA,PB,gamma,maxam1+2,maxam2+2);
125
126 /*--- contract each buffer into appropriate location ---*/
127 double *ints_buf = prim_ints_;
128 for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
129 double norm1 = int_shell1_->coefficient_unnorm(gc1,p1);
130 int am1 = int_shell1_->am(gc1);
131 for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
132 double norm2 = int_shell2_->coefficient_unnorm(gc2,p2);
133 int am2 = int_shell2_->am(gc2);
134 double total_pf = over_pf * norm1 * norm2;
135
136 int k1,l1,m1,k2,l2,m2;
137 FOR_CART(k1,l1,m1,am1)
138 FOR_CART(k2,l2,m2,am2)
139 double x0 = OIX_[k1][k2];
140 double y0 = OIY_[l1][l2];
141 double z0 = OIZ_[m1][m2];
142 double tx = a2*(2*k2+1)*OIX_[k1][k2] - 2*a2*a2*OIX_[k1][k2+2];
143 if (k2 >= 2)
144 tx -= 0.5*k2*(k2-1)*OIX_[k1][k2-2];
145 double ty = a2*(2*l2+1)*OIY_[l1][l2] - 2*a2*a2*OIY_[l1][l2+2];
146 if (l2 >= 2)
147 ty -= 0.5*l2*(l2-1)*OIY_[l1][l2-2];
148 double tz = a2*(2*m2+1)*OIZ_[m1][m2] - 2*a2*a2*OIZ_[m1][m2+2];
149 if (m2 >= 2)
150 tz -= 0.5*m2*(m2-1)*OIZ_[m1][m2-2];
151 *(ints_buf++) += total_pf*(tx*y0*z0 + x0*ty*z0 + x0*y0*tz);
152 END_FOR_CART
153 END_FOR_CART
154
155 }
156 }
157
158 if (bs1_->molecule() != bs2_->molecule()) {
159 // fail();
160 }
161
162 int natom = bs1_->ncenter();
163 for(int atom=0; atom<natom; atom++) {
164 // if charge is 0 - skip to the next one
165 int Z = bs1_->molecule()->Z(atom);
166 if (Z == 0.0)
167 continue;
168 PC[0] = P[0] - bs1_->r(atom,0);
169 PC[1] = P[1] - bs1_->r(atom,1);
170 PC[2] = P[2] - bs1_->r(atom,2);
171 AI_OSrecurs_(AI0_,PA,PB,PC,gamma,maxam1,maxam2);
172
173 /*--- contract each buffer into appropriate location ---*/
174 double *ints_buf = prim_ints_;
175 for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
176 double norm1 = int_shell1_->coefficient_unnorm(gc1,p1);
177 int am1 = int_shell1_->am(gc1);
178 for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
179 double norm2 = int_shell2_->coefficient_unnorm(gc2,p2);
180 int am2 = int_shell2_->am(gc2);
181 double total_pf = over_pf * norm1 * norm2 * Z;
182
183 int k1,l1,m1,k2,l2,m2;
184 FOR_CART(k1,l1,m1,am1)
185 int ind1 = k1*x1weight + l1*y1weight + m1*z1weight;
186 FOR_CART(k2,l2,m2,am2)
187 int ind2 = k2*x2weight + l2*y2weight + m2*z2weight;
188 *ints_buf -= AI0_[ind1][ind2][0] * total_pf;
189 ints_buf++;
190 END_FOR_CART
191 END_FOR_CART
192
193 }
194 }
195 }
196
197 }
198 }
199
200 if (need_cart2sph_transform)
201 transform_contrquartets_(prim_ints_,contr_doublets_);
202
203 if (need_sort_to_shell_doublet)
204 sort_contrdoublets_to_shelldoublet_(contr_doublets_,shell_doublet_);
205}
206
207
208void Int1eCints::hcore_sameam_general_()
209{
210 int tam1 = int_shell1_->am(0);
211 int tam2 = int_shell2_->am(0);
212 int z1weight = 1;
213 int y1weight = tam1 + 1;
214 int x1weight = y1weight * y1weight;
215 int z2weight = 1;
216 int y2weight = tam2 + 1;
217 int x2weight = y2weight * y2weight;
218
219 /* Begin loops over primitives. */
220 for (int p1=0; p1<int_shell1_->nprimitive(); p1++) {
221 double a1 = int_shell1_->exponent(p1);
222 for (int p2=0; p2<int_shell2_->nprimitive(); p2++) {
223 double a2 = int_shell2_->exponent(p2);
224
225 double gamma = a1+a2;
226 double oog = 1.0/gamma;
227 double over_pf = exp(-a1*a2*doublet_info_.AB2*oog)*sqrt(M_PI*oog)*M_PI*oog;
228
229 double P[3], PA[3], PB[3], PC[3];
230 for(int xyz=0; xyz<3; xyz++) {
231 P[xyz] = (a1*doublet_info_.A[xyz] + a2*doublet_info_.B[xyz])*oog;
232 PA[xyz] = P[xyz] - doublet_info_.A[xyz];
233 PB[xyz] = P[xyz] - doublet_info_.B[xyz];
234 }
235
236 OI_OSrecurs_(OIX_,OIY_,OIZ_,PA,PB,gamma,tam1+2,tam2+2);
237
238 /*--- contract each buffer into appropriate location ---*/
239 double *ints_buf = cart_ints_;
240 for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
241 double norm1 = int_shell1_->coefficient_unnorm(gc1,p1);
242 for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
243 double norm2 = int_shell2_->coefficient_unnorm(gc2,p2);
244 double total_pf = over_pf * norm1 * norm2;
245
246 int k1,l1,m1,k2,l2,m2;
247 FOR_CART(k1,l1,m1,tam1)
248 FOR_CART(k2,l2,m2,tam2)
249 double x0 = OIX_[k1][k2];
250 double y0 = OIY_[l1][l2];
251 double z0 = OIZ_[m1][m2];
252 double tx = a2*(2*k2+1)*OIX_[k1][k2] - 2*a2*a2*OIX_[k1][k2+2];
253 if (k2 >= 2)
254 tx -= 0.5*k2*(k2-1)*OIX_[k1][k2-2];
255 double ty = a2*(2*l2+1)*OIY_[l1][l2] - 2*a2*a2*OIY_[l1][l2+2];
256 if (l2 >= 2)
257 ty -= 0.5*l2*(l2-1)*OIY_[l1][l2-2];
258 double tz = a2*(2*m2+1)*OIZ_[m1][m2] - 2*a2*a2*OIZ_[m1][m2+2];
259 if (m2 >= 2)
260 tz -= 0.5*m2*(m2-1)*OIZ_[m1][m2-2];
261 *(ints_buf++) += total_pf*(tx*y0*z0 + x0*ty*z0 + x0*y0*tz);
262 END_FOR_CART
263 END_FOR_CART
264
265 }
266 }
267
268 if (bs1_->molecule() != bs2_->molecule()) {
269 // fail();
270 }
271
272 int natom = bs1_->ncenter();
273 for(int atom=0; atom<natom; atom++) {
274 // if charge is 0 - skip to the next one
275 int Z = bs1_->molecule()->Z(atom);
276 if (Z == 0.0)
277 continue;
278 PC[0] = P[0] - bs1_->r(atom,0);
279 PC[1] = P[1] - bs1_->r(atom,1);
280 PC[2] = P[2] - bs1_->r(atom,2);
281 AI_OSrecurs_(AI0_,PA,PB,PC,gamma,tam1,tam2);
282
283 /*--- contract each buffer into appropriate location ---*/
284 double *ints_buf = cart_ints_;
285 for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
286 double norm1 = int_shell1_->coefficient_unnorm(gc1,p1);
287 for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
288 double norm2 = int_shell2_->coefficient_unnorm(gc2,p2);
289 double total_pf = over_pf * norm1 * norm2 * Z;
290
291 int k1,l1,m1,k2,l2,m2;
292 FOR_CART(k1,l1,m1,tam1)
293 int ind1 = k1*x1weight + l1*y1weight + m1*z1weight;
294 FOR_CART(k2,l2,m2,tam2)
295 int ind2 = k2*x2weight + l2*y2weight + m2*z2weight;
296 *ints_buf -= AI0_[ind1][ind2][0] * total_pf;
297 ints_buf++;
298 END_FOR_CART
299 END_FOR_CART
300
301 }
302 }
303 }
304
305 }
306 }
307
308 /*----------------------------------------------------------------------
309 transform to spherical harmonics and/or resort to the target ordering
310 ----------------------------------------------------------------------*/
311
312 /*--- sort to the target ordering ---*/
313 double *source_ints_buf = cart_ints_;
314 double *target_ints_buf = target_ints_buffer_;
315 int target_bf1_offset = 0;
316 int target_bf2_offset = 0;
317 int nbf2 = int_shell2_->nfunction();
318 for (int gc1=0; gc1<int_shell1_->ncontraction(); gc1++) {
319 int tsize1 = INT_NCART_NN(tam1);
320 for (int gc2=0; gc2<int_shell2_->ncontraction(); gc2++) {
321 int tsize2 = INT_NCART_NN(tam2);
322
323 int k1,l1,m1,k2,l2,m2;
324 int bf1 = 0;
325 FOR_CART(k1,l1,m1,tam1)
326 double *target_ints_buf = target_ints_buffer_ + (target_bf1_offset+bf1)*nbf2 +
327 target_bf2_offset;
328 FOR_CART(k2,l2,m2,tam2)
329 *(target_ints_buf++) = *(source_ints_buf++);
330 END_FOR_CART
331 bf1++;
332 END_FOR_CART
333 target_bf2_offset += tsize2;
334 }
335 target_bf1_offset += tsize1;
336 }
337}
338
339/////////////////////////////////////////////////////////////////////////////
340
341// Local Variables:
342// mode: c++
343// c-file-style: "CLJ"
344// End:
Note: See TracBrowser for help on using the repository browser.