Context Navigation

source: ThirdParty/mpqc_open/src/lib/chemistry/qc/intv3/comp2e.cc

Visit:

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: 65.4 KB

Line
1	//
2	// comp2e.cc
3	//
4	// Copyright (C) 1996 Limit Point Systems, Inc.
5	//
6	// Author: Curtis Janssen <cljanss@limitpt.com>
7	// Maintainer: LPS
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 <stdarg.h>
29
30	#include <util/misc/formio.h>
31	#include <chemistry/qc/intv3/macros.h>
32	#include <chemistry/qc/intv3/flags.h>
33	#include <chemistry/qc/intv3/types.h>
34	#include <chemistry/qc/intv3/int2e.h>
35	#include <chemistry/qc/intv3/utils.h>
36	#include <chemistry/qc/intv3/tformv3.h>
37
38	using namespace std;
39	using namespace sc;
40
41	#undef DER_TIMING
42	#undef EREP_TIMING
43
44	#if defined(DER_TIMING)\|\|defined(EREP_TIMING)
45	# include <util/misc/timer.h>
46	#endif
47
48	static inline void
49	swtch(GaussianBasisSet* &i,GaussianBasisSet* &j)
50	{
51	GaussianBasisSet *tmp;
52	tmp = i;
53	i = j;
54	j = tmp;
55	}
56
57	static inline void
58	sswtch(GaussianShelli,GaussianShellj)
59	{
60	GaussianShell*tmp;
61	tmp = *i;
62	i = j;
63	*j = tmp;
64	}
65
66	static inline void
67	iswtch(int i,int j)
68	{
69	int tmp;
70	tmp = *i;
71	i = j;
72	*j = tmp;
73	}
74
75	static void
76	fail()
77	{
78	ExEnv::errn() << scprintf("failing module:\n%s",__FILE__) << endl;
79	abort();
80	}
81
82	/* This computes the 2erep integrals for a shell quartet
83	* specified by psh1, psh2, psh3, psh4.
84	* The routine int_initialize_2erep must be called before
85	* any integrals can be computed.
86	* This routine may decide to change the shell ordering.
87	* The new ordering is placed in *psh1,4 on exit.
88	* for the derivatives.
89	*/
90	void
91	Int2eV3::erep(int &psh1, int &psh2, int &psh3, int &psh4)
92	{
93	compute_erep(0,&psh1,&psh2,&psh3,&psh4,0,0,0,0);
94	}
95
96	/* This is an alternate interface to int_erep. It takes
97	* as arguments the flags, an integer vector of shell numbers
98	* and an integer vector which will be filled in with size
99	* information, if it is non-NULL. */
100	void
101	Int2eV3::erep(int shells, int sizes)
102	{
103	erep(shells[0],shells[1],shells[2],shells[3]);
104	if (sizes) {
105	sizes[0] = bs1_->shell(shells[0]).nfunction();
106	sizes[1] = bs2_->shell(shells[1]).nfunction();
107	sizes[2] = bs3_->shell(shells[2]).nfunction();
108	sizes[3] = bs4_->shell(shells[3]).nfunction();
109	}
110	}
111
112	/* If we need a computation with adjusted angular momentum, then
113	* this lower level routine can be called instead of int_erep.
114	* The dam{1,2,3,4} arguments given the amount by which the
115	* angular momentum is to adjusted. This differs from libint version
116	* 1 in that it used total angular momentum here.
117	*/
118	void
119	Int2eV3::compute_erep(int flags, int psh1, int psh2, int psh3, int psh4,
120	int dam1, int dam2, int dam3, int dam4)
121	{
122	#ifdef EREP_TIMING
123	char section[30];
124	#endif
125	GaussianBasisSet *pbs1=bs1_.pointer();
126	GaussianBasisSet *pbs2=bs2_.pointer();
127	GaussianBasisSet *pbs3=bs3_.pointer();
128	GaussianBasisSet *pbs4=bs4_.pointer();
129	int size;
130	int ii;
131	int size1, size2, size3, size4;
132	int tam1,tam2,tam3,tam4;
133	int i,j,k,l;
134	int ogc1,ogc2,ogc3,ogc4;
135	int sh1,sh2,sh3,sh4;
136	int am1,am2,am3,am4,am12, am34;
137	int minam1,minam2,minam3,minam4;
138	int redundant_index;
139	int e12,e13e24,e34;
140	int p12,p34,p13p24;
141	int eAB;
142
143	/* Compute the offset shell numbers. */
144	osh1 = *psh1 + bs1_shell_offset_;
145	osh2 = *psh2 + bs2_shell_offset_;
146	osh3 = *psh3 + bs3_shell_offset_;
147	osh4 = *psh4 + bs4_shell_offset_;
148
149	sh1 = *psh1;
150	sh2 = *psh2;
151	sh3 = *psh3;
152	sh4 = *psh4;
153
154	/* Test the arguments to make sure that they are sensible. */
155	if ( sh1 < 0 \|\| sh1 >= bs1_->nbasis()
156	\|\|( !int_unit2 && (sh2 < 0 \|\| sh2 >= bs2_->nbasis()))
157	\|\| sh3 < 0 \|\| sh3 >= bs3_->nbasis()
158	\|\|( !int_unit4 && (sh4 < 0 \|\| sh4 >= bs4_->nbasis()))) {
159	ExEnv::errn() << scprintf("compute_erep has been incorrectly used\n");
160	ExEnv::errn() << scprintf("shells (bounds): %d (%d), %d (%d), %d (%d), %d (%d)\n",
161	sh1,bs1_->nbasis()-1,
162	sh2,(bs2_.null()?0:bs2_->nbasis())-1,
163	sh3,bs3_->nbasis()-1,
164	sh4,(bs4_.null()?0:bs4_->nbasis())-1);
165	fail();
166	}
167
168	/* Set up pointers to the current shells. */
169	int_shell1 = &bs1_->shell(sh1);
170	if (!int_unit2) int_shell2 = &bs2_->shell(sh2);
171	else int_shell2 = int_unit_shell;
172	int_shell3 = &bs3_->shell(sh3);
173	if (!int_unit4) int_shell4 = &bs4_->shell(sh4);
174	else int_shell4 = int_unit_shell;
175
176
177	/* Compute the maximum angular momentum on each centers to
178	* determine the most efficient way to invoke the building and shifting
179	* routines. The minimum angular momentum will be computed at the
180	* same time. */
181	minam1 = int_shell1->min_am();
182	minam2 = int_shell2->min_am();
183	minam3 = int_shell3->min_am();
184	minam4 = int_shell4->min_am();
185	am1 = int_shell1->max_am();
186	am2 = int_shell2->max_am();
187	am3 = int_shell3->max_am();
188	am4 = int_shell4->max_am();
189
190	am1 += dam1; minam1 += dam1;
191	am2 += dam2; minam2 += dam2;
192	am3 += dam3; minam3 += dam3;
193	am4 += dam4; minam4 += dam4;
194	am12 = am1 + am2;
195	am34 = am3 + am4;
196
197	/* if no angular momentum remains on one of the centers return */
198	if (am1 < 0 \|\| am2 < 0 \|\| am3 < 0 \|\| am4 < 0) {
199	return;
200	}
201
202	#ifdef EREP_TIMING
203	sprintf(section,"erep am=%02d",am12+am34);
204	tim_enter(section);
205	tim_enter("setup");
206	#endif
207
208	/* Convert the integral to the most efficient form. */
209	p12 = 0;
210	p34 = 0;
211	p13p24 = 0;
212
213	if (am2 > am1) {
214	p12 = 1;
215	iswtch(&am1,&am2);iswtch(&sh1,&sh2);iswtch(psh1,psh2);iswtch(&osh1,&osh2);
216	iswtch(&dam1,&dam2);
217	iswtch(&minam1,&minam2);
218	sswtch(&int_shell1,&int_shell2);
219	swtch(pbs1,pbs2);
220	}
221	if (am4 > am3) {
222	p34 = 1;
223	iswtch(&am3,&am4);iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
224	iswtch(&dam3,&dam4);
225	iswtch(&minam3,&minam4);
226	sswtch(&int_shell3,&int_shell4);
227	swtch(pbs3,pbs4);
228	}
229	if ((osh1 == osh4) && (osh2 == osh3) && (osh1 != osh2)) {
230	/* Don't make the permutation unless we won't override what was
231	* decided above about p34. */
232	if (am4 == am3) {
233	p34 = 1;
234	iswtch(&am3,&am4);iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
235	iswtch(&dam3,&dam4);
236	iswtch(&minam3,&minam4);
237	sswtch(&int_shell3,&int_shell4);
238	swtch(pbs3,pbs4);
239	}
240	}
241	if ((am34 > am12)\|\|((am34 == am12)&&(minam1 > minam3))) {
242	p13p24 = 1;
243	iswtch(&am1,&am3);iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
244	iswtch(&am2,&am4);iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
245	iswtch(&int_unit2,&int_unit4);
246	iswtch(&am12,&am34);
247	iswtch(&dam1,&dam3);
248	iswtch(&minam1,&minam3);
249	sswtch(&int_shell1,&int_shell3);
250	swtch(pbs1,pbs3);
251	iswtch(&dam2,&dam4);
252	iswtch(&minam2,&minam4);
253	sswtch(&int_shell2,&int_shell4);
254	swtch(pbs2,pbs4);
255	}
256	/* This tries to make centers A and B equivalent, if possible. */
257	else if ( (am3 == am1)
258	&&(am4 == am2)
259	&& !int_unit2
260	&& !int_unit4
261	&&(minam1 == minam3)
262	&&(!( (bs1_ == bs2_)
263	&&(bs1_->shell_to_center(sh1)==bs2_->shell_to_center(sh2))))
264	&&( (bs3_ == bs4_)
265	&&(bs3_->shell_to_center(sh3)==bs4_->shell_to_center(sh4)))) {
266	p13p24 = 1;
267	iswtch(&am1,&am3);iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
268	iswtch(&am2,&am4);iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
269	iswtch(&am12,&am34);
270	iswtch(&dam1,&dam3);
271	iswtch(&minam1,&minam3);
272	sswtch(&int_shell1,&int_shell3);
273	swtch(pbs1,pbs3);
274	iswtch(&dam2,&dam4);
275	iswtch(&minam2,&minam4);
276	sswtch(&int_shell2,&int_shell4);
277	swtch(pbs2,pbs4);
278	}
279
280	if ((pbs1 == pbs2)
281	&&(pbs1->shell_to_center(sh1)==pbs2->shell_to_center(sh2))) {
282	eAB = 1;
283	}
284	else {
285	eAB = 0;
286	}
287
288	/* If the centers were permuted, then the int_expweighted variable may
289	* need to be changed. */
290	if (p12) {
291	iswtch(&int_expweight1,&int_expweight2);
292	}
293	if (p34) {
294	iswtch(&int_expweight3,&int_expweight4);
295	}
296	if (p13p24) {
297	iswtch(&int_expweight1,&int_expweight3);
298	iswtch(&int_expweight2,&int_expweight4);
299	}
300
301	pbs1_ = pbs1;
302	pbs2_ = pbs2;
303	pbs3_ = pbs3;
304	pbs4_ = pbs4;
305
306	int nc1 = int_shell1->ncontraction();
307	int nc2 = int_shell2->ncontraction();
308	int nc3 = int_shell3->ncontraction();
309	int nc4 = int_shell4->ncontraction();
310
311	/* Compute the shell sizes. */
312	for (ii=size1=0; ii<nc1; ii++)
313	size1 += INT_NCART(int_shell1->am(ii)+dam1);
314	for (ii=size2=0; ii<nc2; ii++)
315	size2 += INT_NCART(int_shell2->am(ii)+dam2);
316	for (ii=size3=0; ii<nc3; ii++)
317	size3 += INT_NCART(int_shell3->am(ii)+dam3);
318	for (ii=size4=0; ii<nc4; ii++)
319	size4 += INT_NCART(int_shell4->am(ii)+dam4);
320	size = size1size2size3*size4;
321
322	if (int_integral_storage) {
323	#ifdef EREP_TIMING
324	tim_change("check storage");
325	#endif
326	if (dam1 \|\| dam2 \|\| dam3 \|\| dam4) {
327	ExEnv::errn() << scprintf("cannot use integral storage and dam\n");
328	fail();
329	}
330	if ( !int_unit2
331	&& !int_unit4
332	&& int_have_stored_integral(sh1,sh2,sh3,sh4,p12,p34,p13p24)) {
333	goto post_computation;
334	}
335	}
336
337	/* Buildam up on center 1 and 3. */
338	#ifdef EREP_TIMING
339	tim_change("build");
340	#endif
341	int_buildgcam(minam1,minam2,minam3,minam4,
342	am1,am2,am3,am4,
343	dam1,dam2,dam3,dam4,
344	sh1,sh2,sh3,sh4,
345	eAB);
346	#ifdef EREP_TIMING
347	tim_change("cleanup");
348	#endif
349
350	/* Begin loop over generalized contractions. */
351	ogc1 = 0;
352	for (i=0; i<nc1; i++) {
353	tam1 = int_shell1->am(i) + dam1;
354	if (tam1 < 0) continue;
355	int tsize1 = INT_NCART_NN(tam1);
356	ogc2 = 0;
357	for (j=0; j<nc2; j++) {
358	tam2 = int_shell2->am(j) + dam2;
359	if (tam2 < 0) continue;
360	int tsize2 = INT_NCART_NN(tam2);
361	ogc3 = 0;
362	for (k=0; k<nc3; k++) {
363	tam3 = int_shell3->am(k) + dam3;
364	if (tam3 < 0) continue;
365	int tsize3 = INT_NCART_NN(tam3);
366	ogc4 = 0;
367	for (l=0; l<nc4; l++) {
368	tam4 = int_shell4->am(l) + dam4;
369	if (tam4 < 0) continue;
370	int tsize4 = INT_NCART_NN(tam4);
371
372	#ifdef EREP_TIMING
373	tim_change("shift");
374	#endif
375	/* Shift angular momentum from 1 to 2 and from 3 to 4. */
376	double *shiftbuffer = int_shiftgcam(i,j,k,l,tam1,tam2,tam3,tam4, eAB);
377	#ifdef EREP_TIMING
378	tim_change("cleanup");
379	#endif
380
381	/* Place the integrals in the integral buffer. */
382	/* If permute_ is not set, then repack the integrals while copying. */
383	if ((!permute_)&&(p12\|\|p34\|\|p13p24)) {
384	int pam1,pam2,pam3,pam4;
385	int psize234,psize34;
386	int pogc1,pogc2,pogc3,pogc4;
387	int psize1,psize2,psize3,psize4;
388	pam1 = tam1;
389	pam2 = tam2;
390	pam3 = tam3;
391	pam4 = tam4;
392	pogc1 = ogc1;
393	pogc2 = ogc2;
394	pogc3 = ogc3;
395	pogc4 = ogc4;
396	psize1 = size1;
397	psize2 = size2;
398	psize3 = size3;
399	psize4 = size4;
400	if (p13p24) {
401	iswtch(&pam1,&pam3);
402	iswtch(&pam2,&pam4);
403	iswtch(&pogc1,&pogc3);
404	iswtch(&pogc2,&pogc4);
405	iswtch(&psize1,&psize3);
406	iswtch(&psize2,&psize4);
407	}
408	if (p34) {
409	iswtch(&pam3,&pam4);
410	iswtch(&pogc3,&pogc4);
411	iswtch(&psize3,&psize4);
412	}
413	if (p12) {
414	iswtch(&pam1,&pam2);
415	iswtch(&pogc1,&pogc2);
416	iswtch(&psize1,&psize2);
417	}
418	psize34 = psize4 * psize3;
419	psize234 = psize34 * psize2;
420	redundant_index = 0;
421	int newindexoffset = pogc1psize234 + pogc2psize34 + pogc3*psize4 + pogc4;
422	if (p13p24\|\|p34) {
423	int stride1=psize234;
424	int stride2=psize34;
425	int stride3=psize4;
426	int stride4=1;
427	int tmp;
428	if (p12) {
429	tmp=stride1; stride1=stride2; stride2=tmp;
430	}
431	if (p34) {
432	tmp=stride3; stride3=stride4; stride4=tmp;
433	}
434	if (p13p24) {
435	tmp=stride1; stride1=stride3; stride3=tmp;
436	tmp=stride2; stride2=stride4; stride4=tmp;
437	}
438	int newindex1 = newindexoffset;
439	for (int ci1=0; ci1<tsize1; ci1++) {
440	int newindex12 = newindex1;
441	for (int ci2=0; ci2<tsize2; ci2++) {
442	int newindex123 = newindex12;
443	for (int ci3=0; ci3<tsize3; ci3++) {
444	double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
445	int newindex1234 = newindex123;
446	for (int ci4=0; ci4<tsize4; ci4++) {
447	int_buffer[newindex1234] = tmp_shiftbuffer[ci4];
448	newindex1234 += stride4;
449	}
450	redundant_index+=tsize4;
451	newindex123 += stride3;
452	}
453	newindex12 += stride2;
454	}
455	newindex1 += stride1;
456	}
457	}
458	else if (nc3 == 1 && nc4 == 1) {
459	// this is the p12 only case w/o gen contractions on 3 & 4
460	// this special case collapses the 3rd and 4th indices together
461	for (int ci1=0; ci1<tsize1; ci1++) {
462	for (int ci2=0; ci2<tsize2; ci2++) {
463	int pci1=ci1;
464	int pci2=ci2;
465	if (p12) {
466	int tmp=pci1; pci1=pci2; pci2=tmp;
467	}
468	int newindex123 = newindexoffset + pci1psize234 + pci2psize34;
469	double *tmp_int_buffer = &int_buffer[newindex123];
470	double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
471	for (int ci34=0; ci34<psize34; ci34++)
472	tmp_int_buffer[ci34] = tmp_shiftbuffer[ci34];
473	redundant_index += psize34;
474	}
475	}
476	}
477	else {
478	// this is the p12 only case w/gen. contr. on 3 & 4
479	for (int ci1=0; ci1<tsize1; ci1++) {
480	for (int ci2=0; ci2<tsize2; ci2++) {
481	int pci1=ci1;
482	int pci2=ci2;
483	if (p12) {
484	int tmp=pci1; pci1=pci2; pci2=tmp;
485	}
486	int newindex123 = newindexoffset + pci1psize234 + pci2psize34;
487	for (int ci3=0; ci3<tsize3; ci3++) {
488	double *tmp_int_buffer = &int_buffer[newindex123];
489	double *tmp_shiftbuffer = &shiftbuffer[redundant_index];
490	for (int ci4=0; ci4<tsize4; ci4++) {
491	tmp_int_buffer[ci4] = tmp_shiftbuffer[ci4];
492	}
493	redundant_index += tsize4;
494	newindex123 += psize4;
495	}
496	}
497	}
498	}
499	}
500	else if (nc3 == 1 && nc4 == 1) {
501	// this special case collapses the 3rd and 4th indices together
502	int size34 = size3 * size4;
503	int size234 = size2 * size34;
504	double* redund_ints = shiftbuffer;
505	redundant_index = 0;
506	int newindex1 = ogc1size234 + ogc2size34 + ogc3*size4 + ogc4;
507	for (int ci1=0; ci1<tsize1; ci1++) {
508	int newindex12 = newindex1;
509	for (int ci2=0; ci2<tsize2; ci2++) {
510	double *tmp_int_buffer = &int_buffer[newindex12];
511	double *tmp_redund_ints = &redund_ints[redundant_index];
512	for (int ci34=0; ci34<size34; ci34++)
513	tmp_int_buffer[ci34] = tmp_redund_ints[ci34];
514	redundant_index += size34;
515	newindex12 += size34;
516	}
517	newindex1 += size234;
518	}
519	}
520	else {
521	int size34 = size3 * size4;
522	int size234 = size2 * size34;
523	double* redund_ints = shiftbuffer;
524	redundant_index = 0;
525	int newindex1 = ogc1size234 + ogc2size34 + ogc3*size4 + ogc4;
526	for (int ci1=0; ci1<tsize1; ci1++) {
527	int newindex12 = newindex1;
528	for (int ci2=0; ci2<tsize2; ci2++) {
529	int newindex123 = newindex12;
530	for (int ci3=0; ci3<tsize3; ci3++) {
531	double *tmp_int_buffer = &int_buffer[newindex123];
532	double *tmp_redund_ints = &redund_ints[redundant_index];
533	for (int ci4=0; ci4<tsize4; ci4++) {
534	tmp_int_buffer[ci4] = tmp_redund_ints[ci4];
535	}
536	redundant_index += tsize4;
537	newindex123 += size4;
538	}
539	newindex12 += size34;
540	}
541	newindex1 += size234;
542	}
543	}
544
545	/* End loop over generalized contractions. */
546	ogc4 += tsize4;
547	}
548	ogc3 += tsize3;
549	}
550	ogc2 += tsize2;
551	}
552	ogc1 += tsize1;
553	}
554
555	if ( !int_unit2
556	&& !int_unit4
557	&& int_integral_storage) {
558	#ifdef EREP_TIMING
559	tim_change("maybe store");
560	#endif
561	int_store_integral(sh1,sh2,sh3,sh4,p12,p34,p13p24,size);
562	}
563
564	/* We branch here if an integral was precomputed and the int_buffer
565	* has been already filled. */
566	post_computation:
567
568	#ifdef EREP_TIMING
569	tim_change("post");
570	#endif
571
572	/* Unpermute all of the permuted quantities. */
573	if ((!permute_)&&(p12\|\|p34\|\|p13p24)) {
574	if (p13p24) {
575	iswtch(&sh1,&sh3);iswtch(psh1,psh3);iswtch(&osh1,&osh3);
576	iswtch(&sh2,&sh4);iswtch(psh2,psh4);iswtch(&osh2,&osh4);
577	iswtch(&int_unit2,&int_unit4);
578	iswtch(&am1,&am3);
579	iswtch(&am2,&am4);
580	iswtch(&am12,&am34);
581	sswtch(&int_shell1,&int_shell3);
582	swtch(pbs1,pbs3);
583	sswtch(&int_shell2,&int_shell4);
584	swtch(pbs2,pbs4);
585	iswtch(&int_expweight1,&int_expweight3);
586	iswtch(&int_expweight2,&int_expweight4);
587	}
588	if (p34) {
589	iswtch(&sh3,&sh4);iswtch(psh3,psh4);iswtch(&osh3,&osh4);
590	iswtch(&am3,&am4);
591	sswtch(&int_shell3,&int_shell4);
592	swtch(pbs3,pbs4);
593	iswtch(&int_expweight3,&int_expweight4);
594	}
595	if (p12) {
596	iswtch(&sh1,&sh2);iswtch(psh1,psh2);iswtch(&osh1,&osh2);
597	iswtch(&am1,&am2);
598	sswtch(&int_shell1,&int_shell2);
599	swtch(pbs1,pbs2);
600	iswtch(&int_expweight1,&int_expweight2);
601	}
602	}
603
604	pbs1_ = 0;
605	pbs2_ = 0;
606	pbs3_ = 0;
607	pbs4_ = 0;
608
609	/* Transform to pure am (if requested in the centers structure). */
610	if (!(flags&INT_NOPURE)) {
611	transform_2e(integral_, int_buffer, int_buffer,
612	&bs1_->shell(sh1),
613	int_unit2?int_unit_shell:&bs2_->shell(sh2),
614	&bs3_->shell(sh3),
615	int_unit4?int_unit_shell:&bs4_->shell(sh4));
616	}
617
618	/* Remove the redundant integrals, unless redundant_ is set. */
619	if (!redundant_) {
620	int redundant_offset = 0;
621	int nonredundant_offset = 0;
622	if ((osh1 == osh4)&&(osh2 == osh3)&&(osh1 != osh2)) {
623	ExEnv::errn() << scprintf("nonredundant integrals cannot be generated\n");
624	fail();
625	}
626	e12 = (int_unit2?0:(osh1 == osh2));
627	e13e24 = ((osh1 == osh3)
628	&& ((int_unit2 && int_unit4)
629	\|\| ((int_unit2\|\|int_unit4)?0:(osh2 == osh4))));
630	e34 = (int_unit4?0:(osh3 == osh4));
631	if (e12\|\|e34\|\|e13e24) {
632	nonredundant_erep(int_buffer,e12,e34,e13e24,
633	int_shell1->nfunction(),
634	int_shell2->nfunction(),
635	int_shell3->nfunction(),
636	int_shell4->nfunction(),
637	&redundant_offset,
638	&nonredundant_offset);
639	}
640	}
641
642	#ifdef EREP_TIMING
643	tim_exit("post");
644	tim_exit(section);
645	#endif
646	}
647
648	/* This computes the two electron derivatives for all unique
649	* centers in the passed shell quartet. One center in
650	* the set of unique centers is not included. This can
651	* be computed as minus the sum of the other derivatives.
652	* The list of centers for which integrals were computed can
653	* be determined from the contents of der_centers.
654	* The results are put into the global integral buffer in the
655	* format:
656	* +------------------+
657	* \| dercenter1 +---+ \|
658	* \| \| x \| \|
659	* \| +---+ \|
660	* \| \| y \| \|
661	* \| +---+ \|
662	* \| \| z \| \|
663	* \| +---+ \|
664	* +------------------+
665	* \| dercenter2 +---+ \|
666	* \| \| x \| \|
667	* \| +---+ \|
668	* \| \| y \| \|
669	* \| +---+ \|
670	* \| \| z \| \|
671	* \| +---+ \|
672	* +------------------+
673	* \| dercenter3 +---+ \|
674	* \| \| x \| \|
675	* \| +---+ \|
676	* \| \| y \| \|
677	* \| +---+ \|
678	* \| \| z \| \|
679	* \| +---+ \|
680	* +------------------+
681	*/
682
683	void
684	Int2eV3::erep_all1der(int &psh1, int &psh2, int &psh3, int &psh4,
685	der_centersv3_t *der_centers)
686	{
687	double *current_buffer;
688	int nints;
689	double *user_int_buffer;
690	int omit;
691	GaussianBasisSet *cs[4];
692	int sh[4];
693	int n_unique;
694	int i,j;
695	GaussianShell shell1,shell2,shell3,shell4;
696	GaussianBasisSet ucs[4]; / The centers struct for the unique centers. */
697	int unum[4]; /* The number of times that this unique center occurs. */
698	int uam[4]; /* The total angular momentum on each unique center. */
699	int am[4];
700	int osh[4];
701	int cen[4];
702	int ucen[4];
703	int ncart;
704	double *current_pure_buffer;
705
706	cs[0] = bs1_.pointer();
707	cs[1] = bs2_.pointer();
708	cs[2] = bs3_.pointer();
709	cs[3] = bs4_.pointer();
710
711	sh[0] = psh1;
712	sh[1] = psh2;
713	sh[2] = psh3;
714	sh[3] = psh4;
715
716	/* Set up pointers to the current shells. */
717	shell1 = &bs1_->shell(psh1);
718	shell2 = &bs2_->shell(psh2);
719	shell3 = &bs3_->shell(psh3);
720	shell4 = &bs4_->shell(psh4);
721
722	/* Number of cartesian and pure integrals. */
723	ncart = shell1->ncartesian()*shell2->ncartesian()
724	shell3->ncartesian()shell4->ncartesian();
725	nints = shell1->nfunction()*shell2->nfunction()
726	shell3->nfunction()shell4->nfunction();
727
728	am[0] = shell1->max_am();
729	am[1] = shell2->max_am();
730	am[2] = shell3->max_am();
731	am[3] = shell4->max_am();
732
733	/* Compute the offset shell numbers. */
734	osh[0] = psh1 + bs1_shell_offset_;
735	osh[1] = psh2 + bs2_shell_offset_;
736	osh[2] = psh3 + bs3_shell_offset_;
737	osh[3] = psh4 + bs4_shell_offset_;
738
739	for (i=0; i<4; i++) cen[i] = cs[i]->shell_to_center(sh[i]);
740
741	/* This macro returns true if two shell centers are the same. */
742	#define SC(cs1,shc1,cs2,shc2) (((cs1)==(cs2))&&((shc1)==(shc2)))
743
744	/* Build the list of unique centers structures and shells. */
745	n_unique = 0;
746	for (i=0; i<4; i++) {
747	int unique = 1;
748	for (j=0; j<n_unique; j++) {
749	if (SC(ucs[j],ucen[j],cs[i],cen[i])) {
750	unique = 0;
751	uam[j] += am[i];
752	unum[j]++;
753	break;
754	}
755	}
756	if (unique) {
757	ucs[n_unique] = cs[i];
758	ucen[n_unique] = cen[i];
759	uam[n_unique] = am[i];
760	unum[n_unique] = 1;
761	n_unique++;
762	}
763	}
764
765	/* Choose which unique center is to be omitted from the calculation. */
766	if (n_unique == 1) {
767	omit = 0;
768	}
769	else if (n_unique == 2) {
770	if (unum[0]==3) omit = 0;
771	else if (unum[1]==3) omit = 1;
772	else if (uam[1]>uam[0]) omit = 1;
773	else omit = 0;
774	}
775	else if (n_unique == 3) {
776	if (unum[0]==2) omit = 0;
777	else if (unum[1]==2) omit = 1;
778	else omit = 2;
779	}
780	else {
781	int max = 0;
782	omit = 0;
783	for (i=0; i<4; i++) {
784	if (uam[i]>max) {
785	max = uam[i];
786	omit = i;
787	}
788	}
789	}
790
791	/* Save the location of the int_buffer. */
792	user_int_buffer = int_buffer;
793	int_buffer = int_derint_buffer;
794
795	/* Zero out the result integrals. */
796	for (i=0; i<3(n_unique-1)ncart; i++) user_int_buffer[i] = 0.0;
797
798	/* Loop thru the unique centers, computing the integrals and
799	* skip the derivative on the unique center specified by omit. */
800	der_centers->n = 0;
801	current_buffer = user_int_buffer;
802	for (i=0; i<n_unique; i++) {
803	if (i==omit) continue;
804
805	der_centers->cs[der_centers->n] = ucs[i];
806	der_centers->num[der_centers->n] = ucen[i];
807	der_centers->n++;
808
809	for (j=0; j<4; j++) {
810	if (SC(ucs[i],ucen[i],cs[j],cen[j])) {
811	int old_perm = permute();
812	set_permute(0);
813	compute_erep_1der(0,current_buffer,
814	&psh1,&psh2,&psh3,&psh4,j);
815	set_permute(old_perm);
816	}
817	}
818
819	current_buffer = &current_buffer[3*ncart];
820	}
821
822	/* Put the information about the omitted center into der_centers. */
823	der_centers->ocs = ucs[omit];
824	der_centers->onum = ucen[omit];
825
826	/* Transform to pure am. */
827	current_buffer = user_int_buffer;
828	current_pure_buffer = user_int_buffer;
829	for (i=0; i<3*der_centers->n; i++) {
830	transform_2e(integral_, current_buffer, current_pure_buffer,
831	shell1, shell2, shell3, shell4);
832	current_buffer = &current_buffer[ncart];
833	current_pure_buffer = &current_pure_buffer[nints];
834	}
835
836	/* Eliminate redundant integrals, unless flags specifies otherwise. */
837	current_buffer = user_int_buffer;
838	if (!redundant_) {
839	int redundant_offset = 0;
840	int nonredundant_offset = 0;
841	int e12,e13e24,e34;
842	int i;
843
844	if ((osh[0] == osh[3])&&(osh[1] == osh[2])&&(osh[0] != osh[1])) {
845	ExEnv::errn() << scprintf("nonredundant integrals cannot be generated (1der)\n");
846	fail();
847	}
848
849	/* Shell equivalence information. */
850	e12 = (osh[0] == osh[1]);
851	e13e24 = ((osh[0] == osh[2]) && (osh[1] == osh[3]));
852	e34 = (osh[2] == osh[3]);
853	if (e12\|\|e13e24\|\|e34) {
854	/* Repack x, y, and z integrals. */
855	for (i=0; i<3*der_centers->n; i++) {
856	nonredundant_erep(current_buffer,e12,e34,e13e24,
857	shell1->nfunction(),
858	shell2->nfunction(),
859	shell3->nfunction(),
860	shell4->nfunction(),
861	&redundant_offset,
862	&nonredundant_offset);
863	}
864	}
865	}
866
867	/* Return the integral buffers to their normal state. */
868	int_derint_buffer = int_buffer;
869	int_buffer = user_int_buffer;
870	}
871
872	/* This will call compute_erep
873	* to compute the derivatives in terms of order == 0 integrals.
874	* flags are the flags to the int_comp_erep routine
875	* psh1-4 are pointers to the shell numbers
876	* dercenter is 0, 1, 2, or 3 -- the center within the integral
877	* which we are taking the derivative with respect to.
878	* The results are accumulated in buffer, which cannot be the same
879	* as the current int_buffer.
880	*/
881	void
882	Int2eV3::compute_erep_1der(int flags, double *buffer,
883	int psh1, int psh2, int psh3, int psh4,
884	int dercenter)
885	{
886	int ii;
887	int index;
888	int size1,size2,size3,size4,size1234;
889	int sizem234,sizem34,sizem2,sizem3,sizem4;
890	int sizep234,sizep34,sizep2,sizep3,sizep4;
891	GaussianShell shell1,shell2,shell3,shell4;
892
893	#ifdef DER_TIMING
894	tim_enter("erep_1der");
895	#endif
896
897	/* Set up pointers to the current shells. */
898	shell1 = &bs1_->shell(*psh1);
899	shell2 = &bs2_->shell(*psh2);
900	shell3 = &bs3_->shell(*psh3);
901	shell4 = &bs4_->shell(*psh4);
902
903	if ((dercenter<0) \|\| (dercenter > 3)) {
904	ExEnv::errn() << scprintf("illegal derivative center -- must be 0, 1, 2, or 3\n");
905	fail();
906	}
907
908	/* Offsets for the intermediates with original angular momentum. */
909	for (ii=size1=0; ii<shell1->ncontraction(); ii++)
910	size1 += INT_NCART(shell1->am(ii));
911	for (ii=size2=0; ii<shell2->ncontraction(); ii++)
912	size2 += INT_NCART(shell2->am(ii));
913	for (ii=size3=0; ii<shell3->ncontraction(); ii++)
914	size3 += INT_NCART(shell3->am(ii));
915	for (ii=size4=0; ii<shell4->ncontraction(); ii++)
916	size4 += INT_NCART(shell4->am(ii));
917
918	size1234 = size1size2size3*size4;
919
920	#define DCTEST(n) ((dercenter==n)?1:0)
921	/* Offsets for the intermediates with angular momentum decremented. */
922	for (ii=sizem2=0; ii<shell2->ncontraction(); ii++)
923	sizem2 += INT_NCART(shell2->am(ii)-DCTEST(1));
924	for (ii=sizem3=0; ii<shell3->ncontraction(); ii++)
925	sizem3 += INT_NCART(shell3->am(ii)-DCTEST(2));
926	for (ii=sizem4=0; ii<shell4->ncontraction(); ii++)
927	sizem4 += INT_NCART(shell4->am(ii)-DCTEST(3));
928	sizem34 = sizem4 * sizem3;
929	sizem234 = sizem34 * sizem2;
930
931	/* Offsets for the intermediates with angular momentum incremented. */
932	for (ii=sizep2=0; ii<shell2->ncontraction(); ii++)
933	sizep2 += INT_NCART(shell2->am(ii)+DCTEST(1));
934	for (ii=sizep3=0; ii<shell3->ncontraction(); ii++)
935	sizep3 += INT_NCART(shell3->am(ii)+DCTEST(2));
936	for (ii=sizep4=0; ii<shell4->ncontraction(); ii++)
937	sizep4 += INT_NCART(shell4->am(ii)+DCTEST(3));
938	sizep34 = sizep4 * sizep3;
939	sizep234 = sizep34 * sizep2;
940
941	#ifdef DER_TIMING
942	tim_enter("- erep");
943	#endif
944
945	int old_perm = permute();
946	set_permute(0);
947	int old_red = redundant();
948	set_redundant(1);
949	compute_erep(flags\|INT_NOPURE,
950	psh1,psh2,psh3,psh4,
951	-DCTEST(0),
952	-DCTEST(1),
953	-DCTEST(2),
954	-DCTEST(3));
955	set_permute(old_perm);
956	set_redundant(old_red);
957
958	/* Trouble if cpp is nonANSI. */
959	#define DERLOOP(index,indexp1) \
960	oc##indexp1 = 0;\
961	for ( c##indexp1 =0; c##indexp1 <shell##indexp1->ncontraction(); c##indexp1 ++) {\
962	am[index] = shell##indexp1->am(c##indexp1);\
963	FOR_CART(i[index],j[index],k[index],am[index])
964
965	#define END_DERLOOP(index,indexp1,sign) \
966	END_FOR_CART\
967	oc##indexp1 += INT_NCART(am[index] sign DCTEST(index));\
968	}
969
970	#define ALLDERLOOPS\
971	DERLOOP(0,1)\
972	DERLOOP(1,2)\
973	DERLOOP(2,3)\
974	DERLOOP(3,4)
975
976	#define END_ALLDERLOOPS(sign)\
977	END_DERLOOP(3,4,sign)\
978	END_DERLOOP(2,3,sign)\
979	END_DERLOOP(1,2,sign)\
980	END_DERLOOP(0,1,sign)
981
982	/* Place the contributions into the user integral buffer. */
983	index = 0;
984
985	if (dercenter==0) {
986	int ogc1,ogc1m,gc1,i1,k1,f234,size234;
987	size234=size2size3size4;
988
989	#ifdef DER_TIMING
990	tim_change("- 0");
991	#endif
992	/* The center 0 d/dx integrals */
993	ogc1 = 0;
994	ogc1m = 0;
995	for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
996	int am1 = shell1->am(gc1);
997	// only integrals with x^n n>0 on center 0 contribute
998	// so skip over n==0
999	index += (am1+1)*size234;
1000	int c1 = am1+1;
1001	for (i1=1; i1<=am1; i1++) {
1002	double factor = -i1;
1003	for (k1=0; k1<=am1-i1; k1++) {
1004	int c1xm1 = c1-am1-1;//=INT_CARTINDEX(am1-1,i1-1,j1)
1005	double *tmp_buffer=&buffer[index];
1006	double tmp_int_buffer=&int_buffer[(ogc1m+c1xm1)size234];
1007	for (f234=0; f234<size234; f234++) {
1008	tmp_buffer[f234] += factor * tmp_int_buffer[f234];
1009	}
1010	index+=size234;
1011	c1++;
1012	}
1013	}
1014	ogc1 += c1;
1015	ogc1m += INT_NCART(am1-1);
1016	}
1017
1018	/* The center 0 d/dy integrals */
1019	ogc1 = 0;
1020	ogc1m = 0;
1021	for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1022	int am1 = shell1->am(gc1);
1023	// only integrals with y^n n>0 on center 0 contribute
1024	// so skip over n==0 (by making i1+k1<am1)
1025	int c1 = 0;
1026	for (i1=0; i1<=am1; i1++) {
1027	for (k1=0; k1<=am1-i1-1; k1++) {
1028	double factor = -(am1-i1-k1);
1029	int c1ym1 = c1-i1;//=INT_CARTINDEX(am1-1,i1,j1-1)
1030	double *tmp_buffer=&buffer[index];
1031	double tmp_int_buffer=&int_buffer[(ogc1m+c1ym1)size234];
1032	for (f234=0; f234<size234; f234++) {
1033	tmp_buffer[f234] += factor * tmp_int_buffer[f234];
1034	}
1035	index+=size234;
1036	c1++;
1037	}
1038	// account for the y^n n==0 case by an extra increment
1039	c1++;
1040	index+=size234;
1041	}
1042	ogc1 += c1;
1043	ogc1m += INT_NCART(am1-1);
1044	}
1045
1046	/* The center 0 d/dz integrals */
1047	ogc1 = 0;
1048	ogc1m = 0;
1049	for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1050	int am1 = shell1->am(gc1);
1051	int c1 = 0;
1052	for (i1=0; i1<=am1; i1++) {
1053	// only integrals with z^n n>0 on center 0 contribute
1054	// so skip over n==0
1055	c1++;
1056	index+=size234;
1057	for (k1=1; k1<=am1-i1; k1++) {
1058	double factor = -k1;
1059	int c1zm1 = c1-i1-1;//=INT_CARTINDEX(am1-1,i1,j1)
1060	double *tmp_buffer=&buffer[index];
1061	double tmp_int_buffer=&int_buffer[(ogc1m+c1zm1)size234];
1062	for (f234=0; f234<size234; f234++) {
1063	tmp_buffer[f234] += factor * tmp_int_buffer[f234];
1064	}
1065	index+=size234;
1066	c1++;
1067	}
1068	}
1069	ogc1 += c1;
1070	ogc1m += INT_NCART(am1-1);
1071	}
1072	}
1073	else if (dercenter == 1) {
1074	int ogc2,ogc2m,gc2,i2,k2,f1,f34,size34,size234;
1075	size34 = size3*size4;
1076	size234 = size2size3size4;
1077
1078	#ifdef DER_TIMING
1079	tim_change("- 1");
1080	#endif
1081	/* The center 1 d/dx integrals */
1082	ogc2 = 0;
1083	ogc2m = 0;
1084	for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1085	int am2 = shell2->am(gc2);
1086	// only integrals with x^n n>0 on center 1 contribute
1087	// so skip over n==0
1088	int c2 = am2+1;
1089	for (i2=1; i2<=am2; i2++) {
1090	double factor = -i2;
1091	for (k2=0; k2<=am2-i2; k2++) {
1092	int c2xm1 = c2-am2-1;//=INT_CARTINDEX(am2-1,i2-1,j2)
1093	int buffer_index = (ogc2+c2)*size34;
1094	int int_buffer_index = (ogc2m+c2xm1)*size34;
1095	for (f1=0; f1<size1; f1++) {
1096	double *tmp_buffer=&buffer[buffer_index];
1097	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1098	for (f34=0; f34<size34; f34++) {
1099	tmp_buffer[f34] += factor * tmp_int_buffer[f34];
1100	}
1101	buffer_index += size234;
1102	int_buffer_index += sizem234;
1103	}
1104	c2++;
1105	}
1106	}
1107	ogc2 += c2;
1108	ogc2m += INT_NCART(am2-1);
1109	}
1110	index += size1234;
1111
1112	/* The center 1 d/dy integrals */
1113	ogc2 = 0;
1114	ogc2m = 0;
1115	for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1116	int am2 = shell2->am(gc2);
1117	// only integrals with y^n n>0 on center 1 contribute
1118	// so skip over n==0
1119	int c2 = 0;
1120	for (i2=0; i2<=am2; i2++) {
1121	for (k2=0; k2<=am2-i2-1; k2++) {
1122	double factor = -(am2-k2-i2);
1123	int c2ym1 = c2-i2;//=INT_CARTINDEX(am2-1,i2,j2-1)
1124	int buffer_index = size1234 + (ogc2+c2)*size34;
1125	int int_buffer_index = (ogc2m+c2ym1)*size34;
1126	for (f1=0; f1<size1; f1++) {
1127	double *tmp_buffer=&buffer[buffer_index];
1128	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1129	for (f34=0; f34<size34; f34++) {
1130	tmp_buffer[f34] += factor * tmp_int_buffer[f34];
1131	}
1132	buffer_index += size234;
1133	int_buffer_index += sizem234;
1134	}
1135	c2++;
1136	}
1137	// account for the y^n n==0 case by an extra increment
1138	c2++;
1139	}
1140	ogc2 += c2;
1141	ogc2m += INT_NCART(am2-1);
1142	}
1143	index += size1234;
1144
1145	/* The center 1 d/dz integrals */
1146	ogc2 = 0;
1147	ogc2m = 0;
1148	for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1149	int am2 = shell2->am(gc2);
1150	// only integrals with z^n n>0 on center 1 contribute
1151	// so skip over n==0
1152	int c2 = 0;
1153	for (i2=0; i2<=am2; i2++) {
1154	// account for the z^n n==0 case by an extra increment
1155	c2++;
1156	for (k2=1; k2<=am2-i2; k2++) {
1157	double factor = -k2;
1158	int c2zm1 = c2-i2-1;//=INT_CARTINDEX(am2-1,i2,j2-1)
1159	int buffer_index = size1234+size1234+(ogc2+c2)*size34;
1160	int int_buffer_index = (ogc2m+c2zm1)*size34;
1161	for (f1=0; f1<size1; f1++) {
1162	double *tmp_buffer=&buffer[buffer_index];
1163	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1164	for (f34=0; f34<size34; f34++) {
1165	tmp_buffer[f34] += factor * tmp_int_buffer[f34];
1166	}
1167	buffer_index += size234;
1168	int_buffer_index += sizem234;
1169	}
1170	c2++;
1171	}
1172	}
1173	ogc2 += c2;
1174	ogc2m += INT_NCART(am2-1);
1175	}
1176	index += size1234;
1177	}
1178	else if (dercenter == 2) {
1179	int ogc3,ogc3m,gc3,i3,k3,f12,f4,size12,size34;
1180	size12 = size1*size2;
1181	size34 = size3*size4;
1182
1183	#ifdef DER_TIMING
1184	tim_change("- 2");
1185	#endif
1186	/* The center 2 d/dx integrals */
1187	ogc3 = 0;
1188	ogc3m = 0;
1189	for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1190	int am3 = shell3->am(gc3);
1191	// only integrals with x^n n>0 on center 2 contribute
1192	// so skip over n==0
1193	int c3 = am3+1;
1194	for (i3=1; i3<=am3; i3++) {
1195	double factor = -i3;
1196	for (k3=0; k3<=am3-i3; k3++) {
1197	int c3xm1 = c3-am3-1;//=INT_CARTINDEX(am3-1,i3-1,j3)
1198	int buffer_index = (ogc3+c3)*size4;
1199	int int_buffer_index = (ogc3m+c3xm1)*size4;
1200	for (f12=0; f12<size12; f12++) {
1201	double *tmp_buffer=&buffer[buffer_index];
1202	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1203	for (f4=0; f4<size4; f4++) {
1204	tmp_buffer[f4] += factor * tmp_int_buffer[f4];
1205	}
1206	buffer_index += size34;
1207	int_buffer_index += sizem34;
1208	}
1209	c3++;
1210	}
1211	}
1212	ogc3 += c3;
1213	ogc3m += INT_NCART(am3-1);
1214	}
1215	index += size1234;
1216
1217	/* The center 2 d/dy integrals */
1218	ogc3 = 0;
1219	ogc3m = 0;
1220	for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1221	int am3 = shell3->am(gc3);
1222	// only integrals with y^n n>0 on center 2 contribute
1223	// so skip over n==0
1224	int c3 = 0;
1225	for (i3=0; i3<=am3; i3++) {
1226	for (k3=0; k3<=am3-i3-1; k3++) {
1227	double factor = -(am3-k3-i3);
1228	int c3ym1 = c3-i3;//=INT_CARTINDEX(am3-1,i3,j3-1)
1229	int buffer_index = size1234 + (ogc3+c3)*size4;
1230	int int_buffer_index = (ogc3m+c3ym1)*size4;
1231	for (f12=0; f12<size12; f12++) {
1232	double *tmp_buffer=&buffer[buffer_index];
1233	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1234	for (f4=0; f4<size4; f4++) {
1235	tmp_buffer[f4] += factor * tmp_int_buffer[f4];
1236	}
1237	buffer_index += size34;
1238	int_buffer_index += sizem34;
1239	}
1240	c3++;
1241	}
1242	// account for the y^n n==0 case by an extra increment
1243	c3++;
1244	}
1245	ogc3 += c3;
1246	ogc3m += INT_NCART(am3-1);
1247	}
1248	index += size1234;
1249
1250	/* The center 2 d/dz integrals */
1251	ogc3 = 0;
1252	ogc3m = 0;
1253	for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1254	int am3 = shell3->am(gc3);
1255	// only integrals with z^n n>0 on center 2 contribute
1256	// so skip over n==0
1257	int c3 = 0;
1258	for (i3=0; i3<=am3; i3++) {
1259	// account for the z^n n==0 case by an extra increment
1260	c3++;
1261	for (k3=1; k3<=am3-i3; k3++) {
1262	double factor = -k3;
1263	int c3zm1 = c3-i3-1;//=INT_CARTINDEX(am3-1,i3,j3)
1264	int buffer_index = size1234+size1234+(ogc3+c3)*size4;
1265	int int_buffer_index = (ogc3m+c3zm1)*size4;
1266	for (f12=0; f12<size12; f12++) {
1267	double *tmp_buffer=&buffer[buffer_index];
1268	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1269	for (f4=0; f4<size4; f4++) {
1270	tmp_buffer[f4] += factor * tmp_int_buffer[f4];
1271	}
1272	buffer_index += size34;
1273	int_buffer_index += sizem34;
1274	}
1275	c3++;
1276	}
1277	}
1278	ogc3 += c3;
1279	ogc3m += INT_NCART(am3-1);
1280	}
1281	index += size1234;
1282	}
1283	else if (dercenter == 3) {
1284	int ogc4,ogc4m,gc4,i4,k4,f123,size123;
1285	size123 = size1size2size3;
1286
1287	#ifdef DER_TIMING
1288	tim_change("- 3");
1289	#endif
1290	/* The center 3 d/dx integrals */
1291	ogc4 = 0;
1292	ogc4m = 0;
1293	for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1294	int am4 = shell4->am(gc4);
1295	// only integrals with x^n n>0 on center 3 contribute
1296	// so skip over n==0
1297	int c4 = am4+1;
1298	for (i4=1; i4<=am4; i4++) {
1299	double factor = -i4;
1300	for (k4=0; k4<=am4-i4; k4++) {
1301	int c4xm1 = c4-am4-1;//=INT_CARTINDEX(am4-1,i4-1,j4)
1302	int buffer_index = ogc4+c4;
1303	int int_buffer_index = ogc4m+c4xm1;
1304	for (f123=0; f123<size123; f123++) {
1305	buffer[buffer_index] += factor * int_buffer[int_buffer_index];
1306	buffer_index += size4;
1307	int_buffer_index += sizem4;
1308	}
1309	c4++;
1310	}
1311	}
1312	ogc4 += c4;
1313	ogc4m += INT_NCART(am4-1);
1314	}
1315	index += size1234;
1316
1317	/* The center 3 d/dy integrals */
1318	ogc4 = 0;
1319	ogc4m = 0;
1320	for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1321	int am4 = shell4->am(gc4);
1322	// only integrals with y^n n>0 on center 3 contribute
1323	// so skip over n==0
1324	int c4 = 0;
1325	for (i4=0; i4<=am4; i4++) {
1326	for (k4=0; k4<=am4-i4-1; k4++) {
1327	double factor = -(am4-k4-i4);
1328	int c4ym1 = c4-i4;//=INT_CARTINDEX(am4-1,i4,j4-1)
1329	int buffer_index = size1234 + ogc4+c4;
1330	int int_buffer_index = ogc4m+c4ym1;
1331	for (f123=0; f123<size123; f123++) {
1332	buffer[buffer_index] += factor * int_buffer[int_buffer_index];
1333	buffer_index += size4;
1334	int_buffer_index += sizem4;
1335	}
1336	c4++;
1337	}
1338	// account for the y^n n==0 case by an extra increment
1339	c4++;
1340	}
1341	ogc4 += c4;
1342	ogc4m += INT_NCART(am4-1);
1343	}
1344	index += size1234;
1345
1346	/* The center 3 d/dz integrals */
1347	ogc4 = 0;
1348	ogc4m = 0;
1349	for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1350	int am4 = shell4->am(gc4);
1351	// only integrals with z^n n>0 on center 3 contribute
1352	// so skip over n==0
1353	int c4 = 0;
1354	for (i4=0; i4<=am4; i4++) {
1355	// account for the z^n n==0 case by an extra increment
1356	c4++;
1357	for (k4=1; k4<=am4-i4; k4++) {
1358	double factor = -k4;
1359	int c4zm1 = c4-i4-1;//=INT_CARTINDEX(am4-1,i4,j4-1)
1360	int buffer_index = size1234+size1234+ogc4+c4;
1361	int int_buffer_index = ogc4m+c4zm1;
1362	for (f123=0; f123<size123; f123++) {
1363	buffer[buffer_index] += factor * int_buffer[int_buffer_index];
1364	buffer_index += size4;
1365	int_buffer_index += sizem4;
1366	}
1367	c4++;
1368	}
1369	}
1370	ogc4 += c4;
1371	ogc4m += INT_NCART(am4-1);
1372	}
1373	index += size1234;
1374	}
1375
1376	#ifdef DER_TIMING
1377	tim_change("+ erep");
1378	#endif
1379
1380	/* Compute the next contribution to the integrals. */
1381	/* Tell the build routine that we need an exponent weighted contraction
1382	* with the exponents taken from the dercenter and adjust the
1383	* angular momentum of dercenter to the needed value. */
1384	if (dercenter==0) int_expweight1 = 1;
1385	else if (dercenter==1) int_expweight2 = 1;
1386	else if (dercenter==2) int_expweight3 = 1;
1387	else if (dercenter==3) int_expweight4 = 1;
1388	old_perm = permute();
1389	set_permute(0);
1390	old_red = redundant();
1391	set_redundant(1);
1392	compute_erep(flags\|INT_NOPURE,
1393	psh1,psh2,psh3,psh4,
1394	DCTEST(0),
1395	DCTEST(1),
1396	DCTEST(2),
1397	DCTEST(3));
1398	set_permute(old_perm);
1399	set_redundant(old_red);
1400	if (dercenter==0) int_expweight1 = 0;
1401	else if (dercenter==1) int_expweight2 = 0;
1402	else if (dercenter==2) int_expweight3 = 0;
1403	else if (dercenter==3) int_expweight4 = 0;
1404
1405	/* Place the contributions into the user integral buffer. */
1406	index = 0;
1407	if (dercenter==0) {
1408	int ogc1,ogc1p,gc1,i1,k1,f234,size234;
1409	size234=size2size3size4;
1410
1411	#ifdef DER_TIMING
1412	tim_change("+ 0");
1413	#endif
1414	/* The center 0 d/dx integrals */
1415	ogc1 = 0;
1416	ogc1p = 0;
1417	for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1418	int am1 = shell1->am(gc1);
1419	int c1 = 0;
1420	for (i1=0; i1<=am1; i1++) {
1421	for (k1=0; k1<=am1-i1; k1++) {
1422	int c1xp1 = c1+am1+2;//=INT_CARTINDEX(am1+1,i1+1,j1)
1423	double *tmp_buffer=&buffer[index];
1424	double tmp_int_buffer=&int_buffer[(ogc1p+c1xp1)size234];
1425	for (f234=0; f234<size234; f234++) {
1426	tmp_buffer[f234] += tmp_int_buffer[f234];
1427	}
1428	index+=size234;
1429	c1++;
1430	}
1431	}
1432	ogc1 += c1;
1433	ogc1p += INT_NCART(am1+1);
1434	}
1435
1436	/* The center 0 d/dy integrals */
1437	ogc1 = 0;
1438	ogc1p = 0;
1439	for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1440	int am1 = shell1->am(gc1);
1441	int c1 = 0;
1442	for (i1=0; i1<=am1; i1++) {
1443	for (k1=0; k1<=am1-i1; k1++) {
1444	int c1yp1 = c1+i1;//=INT_CARTINDEX(am1+1,i1,j1+1)
1445	double *tmp_buffer=&buffer[index];
1446	double tmp_int_buffer=&int_buffer[(ogc1p+c1yp1)size234];
1447	for (f234=0; f234<size234; f234++) {
1448	tmp_buffer[f234] += tmp_int_buffer[f234];
1449	}
1450	index+=size234;
1451	c1++;
1452	}
1453	}
1454	ogc1 += c1;
1455	ogc1p += INT_NCART(am1+1);
1456	}
1457
1458	/* The center 0 d/dz integrals */
1459	ogc1 = 0;
1460	ogc1p = 0;
1461	for (gc1=0; gc1<shell1->ncontraction(); gc1++) {
1462	int am1 = shell1->am(gc1);
1463	int c1 = 0;
1464	for (i1=0; i1<=am1; i1++) {
1465	for (k1=0; k1<=am1-i1; k1++) {
1466	int c1zp1 = c1+i1+1;//=INT_CARTINDEX(am1+1,i1,j1)
1467	double *tmp_buffer=&buffer[index];
1468	double tmp_int_buffer=&int_buffer[(ogc1p+c1zp1)size234];
1469	for (f234=0; f234<size234; f234++) {
1470	tmp_buffer[f234] += tmp_int_buffer[f234];
1471	}
1472	index+=size234;
1473	c1++;
1474	}
1475	}
1476	ogc1 += c1;
1477	ogc1p += INT_NCART(am1+1);
1478	}
1479	}
1480	else if (dercenter == 1) {
1481	int ogc2,ogc2p,gc2,i2,k2,f1,f34,size34,size234;
1482	size34 = size3*size4;
1483	size234 = size2size3size4;
1484
1485	#ifdef DER_TIMING
1486	tim_change("+ 1");
1487	#endif
1488	/* The center 1 d/dx integrals */
1489	ogc2 = 0;
1490	ogc2p = 0;
1491	for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1492	int am2 = shell2->am(gc2);
1493	int c2=0;
1494	for (i2=0; i2<=am2; i2++) {
1495	for (k2=0; k2<=am2-i2; k2++) {
1496	int c2xp1 = c2+am2+2;//=INT_CARTINDEX(am2+1,i2+1,j2)
1497	int buffer_index = (ogc2+c2)*size34;
1498	int int_buffer_index = (ogc2p+c2xp1)*size34;
1499	for (f1=0; f1<size1; f1++) {
1500	double *tmp_buffer=&buffer[buffer_index];
1501	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1502	for (f34=0; f34<size34; f34++) {
1503	tmp_buffer[f34] += tmp_int_buffer[f34];
1504	}
1505	buffer_index += size234;
1506	int_buffer_index += sizep234;
1507	}
1508	c2++;
1509	}
1510	}
1511	ogc2 += c2;
1512	ogc2p += INT_NCART(am2+1);
1513	}
1514	index += size1234;
1515
1516	/* The center 1 d/dy integrals */
1517	ogc2 = 0;
1518	ogc2p = 0;
1519	for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1520	int am2 = shell2->am(gc2);
1521	int c2 = 0;
1522	for (i2=0; i2<=am2; i2++) {
1523	for (k2=0; k2<=am2-i2; k2++) {
1524	int c2yp1 = c2+i2;//=INT_CARTINDEX(am2+1,i2,j2+1)
1525	int buffer_index = size1234 + (ogc2+c2)*size34;
1526	int int_buffer_index = (ogc2p+c2yp1)*size34;
1527	for (f1=0; f1<size1; f1++) {
1528	double *tmp_buffer=&buffer[buffer_index];
1529	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1530	for (f34=0; f34<size34; f34++) {
1531	tmp_buffer[f34] += tmp_int_buffer[f34];
1532	}
1533	buffer_index += size234;
1534	int_buffer_index += sizep234;
1535	}
1536	c2++;
1537	}
1538	}
1539	ogc2 += c2;
1540	ogc2p += INT_NCART(am2+1);
1541	}
1542	index += size1234;
1543
1544	/* The center 1 d/dz integrals */
1545	ogc2 = 0;
1546	ogc2p = 0;
1547	for (gc2=0; gc2<shell2->ncontraction(); gc2++) {
1548	int am2 = shell2->am(gc2);
1549	int c2 = 0;
1550	for (i2=0; i2<=am2; i2++) {
1551	for (k2=0; k2<=am2-i2; k2++) {
1552	int c2zp1 = c2+i2+1;//=INT_CARTINDEX(am2+1,i2,j2+1)
1553	int buffer_index = size1234+size1234+(ogc2+c2)*size34;
1554	int int_buffer_index = (ogc2p+c2zp1)*size34;
1555	for (f1=0; f1<size1; f1++) {
1556	double *tmp_buffer=&buffer[buffer_index];
1557	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1558	for (f34=0; f34<size34; f34++) {
1559	tmp_buffer[f34] += tmp_int_buffer[f34];
1560	}
1561	buffer_index += size234;
1562	int_buffer_index += sizep234;
1563	}
1564	c2++;
1565	}
1566	}
1567	ogc2 += c2;
1568	ogc2p += INT_NCART(am2+1);
1569	}
1570	index += size1234;
1571	}
1572	else if (dercenter == 2) {
1573	int ogc3,ogc3p,gc3,i3,k3,f12,f4,size12,size34;
1574	size12 = size1*size2;
1575	size34 = size3*size4;
1576
1577	#ifdef DER_TIMING
1578	tim_change("+ 2");
1579	#endif
1580	/* The center 2 d/dx integrals */
1581	ogc3 = 0;
1582	ogc3p = 0;
1583	for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1584	int am3 = shell3->am(gc3);
1585	int c3 = 0;
1586	for (i3=0; i3<=am3; i3++) {
1587	for (k3=0; k3<=am3-i3; k3++) {
1588	int c3xp1 = c3+am3+2;//=INT_CARTINDEX(am3+1,i3+1,j3)
1589	int buffer_index = (ogc3+c3)*size4;
1590	int int_buffer_index = (ogc3p+c3xp1)*size4;
1591	for (f12=0; f12<size12; f12++) {
1592	double *tmp_buffer=&buffer[buffer_index];
1593	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1594	for (f4=0; f4<size4; f4++) {
1595	tmp_buffer[f4] += tmp_int_buffer[f4];
1596	}
1597	buffer_index += size34;
1598	int_buffer_index += sizep34;
1599	}
1600	c3++;
1601	}
1602	}
1603	ogc3 += c3;
1604	ogc3p += INT_NCART(am3+1);
1605	}
1606	index += size1234;
1607
1608	/* The center 2 d/dy integrals */
1609	ogc3 = 0;
1610	ogc3p = 0;
1611	for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1612	int am3 = shell3->am(gc3);
1613	int c3 = 0;
1614	for (i3=0; i3<=am3; i3++) {
1615	for (k3=0; k3<=am3-i3; k3++) {
1616	int c3yp1 = c3+i3;//=INT_CARTINDEX(am3+1,i3,j3+1)
1617	int buffer_index = size1234 + (ogc3+c3)*size4;
1618	int int_buffer_index = (ogc3p+c3yp1)*size4;
1619	for (f12=0; f12<size12; f12++) {
1620	double *tmp_buffer=&buffer[buffer_index];
1621	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1622	for (f4=0; f4<size4; f4++) {
1623	tmp_buffer[f4] += tmp_int_buffer[f4];
1624	}
1625	buffer_index += size34;
1626	int_buffer_index += sizep34;
1627	}
1628	c3++;
1629	}
1630	}
1631	ogc3 += c3;
1632	ogc3p += INT_NCART(am3+1);
1633	}
1634	index += size1234;
1635
1636	/* The center 2 d/dz integrals */
1637	ogc3 = 0;
1638	ogc3p = 0;
1639	for (gc3=0; gc3<shell3->ncontraction(); gc3++) {
1640	int am3 = shell3->am(gc3);
1641	int c3 = 0;
1642	for (i3=0; i3<=am3; i3++) {
1643	for (k3=0; k3<=am3-i3; k3++) {
1644	int c3zp1 = c3+i3+1;//=INT_CARTINDEX(am3+1,i3,j3)
1645	int buffer_index = size1234+size1234+(ogc3+c3)*size4;
1646	int int_buffer_index = (ogc3p+c3zp1)*size4;
1647	for (f12=0; f12<size12; f12++) {
1648	double *tmp_buffer=&buffer[buffer_index];
1649	double *tmp_int_buffer=&int_buffer[int_buffer_index];
1650	for (f4=0; f4<size4; f4++) {
1651	tmp_buffer[f4] += tmp_int_buffer[f4];
1652	}
1653	buffer_index += size34;
1654	int_buffer_index += sizep34;
1655	}
1656	c3++;
1657	}
1658	}
1659	ogc3 += c3;
1660	ogc3p += INT_NCART(am3+1);
1661	}
1662	index += size1234;
1663	}
1664	else if (dercenter == 3) {
1665	int ogc4,ogc4p,gc4,i4,k4,f123,size123;
1666	size123 = size1size2size3;
1667
1668	#ifdef DER_TIMING
1669	tim_change("+ 3");
1670	#endif
1671	/* The center 3 d/dx integrals */
1672	ogc4 = 0;
1673	ogc4p = 0;
1674	for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1675	int am4 = shell4->am(gc4);
1676	int c4 = 0;
1677	for (i4=0; i4<=am4; i4++) {
1678	for (k4=0; k4<=am4-i4; k4++) {
1679	int c4xp1 = c4+am4+2;//=INT_CARTINDEX(am4+1,i4+1,j4)
1680	int buffer_index = ogc4+c4;
1681	int int_buffer_index = ogc4p+c4xp1;
1682	for (f123=0; f123<size123; f123++) {
1683	buffer[buffer_index] += int_buffer[int_buffer_index];
1684	buffer_index += size4;
1685	int_buffer_index += sizep4;
1686	}
1687	c4++;
1688	}
1689	}
1690	ogc4 += c4;
1691	ogc4p += INT_NCART(am4+1);
1692	}
1693	index += size1234;
1694
1695	/* The center 3 d/dy integrals */
1696	ogc4 = 0;
1697	ogc4p = 0;
1698	for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1699	int am4 = shell4->am(gc4);
1700	int c4 = 0;
1701	for (i4=0; i4<=am4; i4++) {
1702	for (k4=0; k4<=am4-i4; k4++) {
1703	int c4yp1 = c4+i4;//=INT_CARTINDEX(am4+1,i4,j4+1)
1704	int buffer_index = size1234 + ogc4+c4;
1705	int int_buffer_index = ogc4p+c4yp1;
1706	for (f123=0; f123<size123; f123++) {
1707	buffer[buffer_index] += int_buffer[int_buffer_index];
1708	buffer_index += size4;
1709	int_buffer_index += sizep4;
1710	}
1711	c4++;
1712	}
1713	}
1714	ogc4 += c4;
1715	ogc4p += INT_NCART(am4+1);
1716	}
1717	index += size1234;
1718
1719	/* The center 3 d/dz integrals */
1720	ogc4 = 0;
1721	ogc4p = 0;
1722	for (gc4=0; gc4<shell4->ncontraction(); gc4++) {
1723	int am4 = shell4->am(gc4);
1724	int c4 = 0;
1725	for (i4=0; i4<=am4; i4++) {
1726	for (k4=0; k4<=am4-i4; k4++) {
1727	int c4zp1 = c4+i4+1;//=INT_CARTINDEX(am4+1,i4,j4)
1728	int buffer_index = size1234+size1234+ogc4+c4;
1729	int int_buffer_index = ogc4p+c4zp1;
1730	for (f123=0; f123<size123; f123++) {
1731	buffer[buffer_index] += int_buffer[int_buffer_index];
1732	buffer_index += size4;
1733	int_buffer_index += sizep4;
1734	}
1735	c4++;
1736	}
1737	}
1738	ogc4 += c4;
1739	ogc4p += INT_NCART(am4+1);
1740	}
1741	index += size1234;
1742	}
1743	#ifdef DER_TIMING
1744	tim_exit(0);
1745	tim_exit(0);
1746	#endif
1747	}
1748
1749	void
1750	Int2eV3::nonredundant_erep(double *buffer, int e12, int e34, int e13e24,
1751	int n1, int n2, int n3, int n4,
1752	int red_off, int nonred_off)
1753	{
1754	int nonredundant_index;
1755	int i,j,k,l;
1756
1757	double redundant_ptr = &buffer[red_off];
1758	double nonredundant_ptr = &buffer[nonred_off];
1759
1760	nonredundant_index = 0;
1761	int n34 = n3*n4;
1762	for (i=0; i<n1; i++) {
1763	int jmax = INT_MAX2(e12,i,n2);
1764	for (j=0; j<=jmax; j++) {
1765	int kmax = INT_MAX3(e13e24,i,n3);
1766	for (k=0; k<=kmax; k++) {
1767	int lmax = INT_MAX4(e13e24,e34,i,j,k,n4);
1768	for (l=0; l<=lmax; l++) {
1769	nonredundant_ptr[l] = redundant_ptr[l];
1770	}
1771	redundant_ptr += n4;
1772	nonredundant_index += lmax+1;
1773	nonredundant_ptr += lmax+1;
1774	}
1775	redundant_ptr += (n3-(kmax+1))*n4;
1776	}
1777	redundant_ptr += (n2-(jmax+1))*n34;
1778	}
1779	red_off += n1n2*n34;
1780	*nonred_off += nonredundant_index;
1781	}
1782
1783	/* This is an alternate interface to int_erep_all1der. It takes
1784	* as arguments the flags, an integer vector of shell numbers
1785	* and an integer vector which will be filled in with size
1786	* information, if it is non-NULL, and the dercenters pointer. */
1787	void
1788	Int2eV3::erep_all1der(int shells, int sizes,
1789	der_centersv3_t *dercenters)
1790	{
1791	erep_all1der(shells[0],shells[1],shells[2],shells[3],
1792	dercenters);
1793	if (sizes) {
1794	sizes[0] = bs1_->shell(shells[0]).nfunction();
1795	sizes[1] = bs2_->shell(shells[1]).nfunction();
1796	sizes[2] = bs3_->shell(shells[2]).nfunction();
1797	sizes[3] = bs4_->shell(shells[3]).nfunction();
1798	}
1799	}
1800
1801	void
1802	Int2eV3::int_erep_bound1der(int flags, int bsh1, int bsh2, int *size)
1803	{
1804	double *current_buffer;
1805	int nints;
1806	double *user_int_buffer;
1807	int i;
1808	GaussianShell shell1,shell2,shell3,shell4;
1809	int osh[4];
1810	int sh1 = bsh1;
1811	int sh2 = bsh2;
1812	int sh3 = bsh1;
1813	int sh4 = bsh2;
1814	int *psh1 = &sh1;
1815	int *psh2 = &sh2;
1816	int *psh3 = &sh3;
1817	int *psh4 = &sh4;
1818	int ncart;
1819	double *current_pure_buffer;
1820
1821	/* Set up pointers to the current shells. */
1822	shell1 = &bs1_->shell(*psh1);
1823	shell2 = &bs2_->shell(*psh2);
1824	shell3 = &bs3_->shell(*psh3);
1825	shell4 = &bs4_->shell(*psh4);
1826
1827	/* Number of cartesian and pure integrals. */
1828	ncart = shell1->ncartesian()*shell2->ncartesian()
1829	shell3->ncartesian()shell4->ncartesian();
1830	nints = shell1->nfunction() * shell2->nfunction()
1831	* shell3->nfunction() * shell4->nfunction();
1832
1833	/* Compute the offset shell numbers. */
1834	osh[0] = *psh1 + bs1_shell_offset_;
1835	osh[1] = *psh2 + bs2_shell_offset_;
1836	osh[2] = *psh3 + bs3_shell_offset_;
1837	osh[3] = *psh4 + bs4_shell_offset_;
1838
1839	/* Save the location of the int_buffer. */
1840	user_int_buffer = int_buffer;
1841	int_buffer = int_derint_buffer;
1842
1843	/* Zero out the result integrals. */
1844	for (i=0; i<3*ncart; i++) user_int_buffer[i] = 0.0;
1845
1846	/* Set the size so it is available to the caller. */
1847	size = nints 3;
1848
1849	current_buffer = user_int_buffer;
1850	int old_perm = permute();
1851	compute_erep_bound1der(flags\|INT_NOPURE,current_buffer,
1852	psh1,psh2,psh3,psh4);
1853	set_permute(old_perm);
1854
1855	/* Transform to pure am. */
1856	current_buffer = user_int_buffer;
1857	current_pure_buffer = user_int_buffer;
1858	for (i=0; i<3; i++) {
1859	transform_2e(integral_, current_buffer, current_pure_buffer,
1860	&bs1_->shell(sh1),
1861	&bs2_->shell(sh2),
1862	&bs3_->shell(sh3),
1863	&bs4_->shell(sh4));
1864	current_buffer = &current_buffer[ncart];
1865	current_pure_buffer = &current_pure_buffer[nints];
1866	}
1867
1868	/* Eliminate redundant integrals, unless flags specifies otherwise. */
1869	current_buffer = user_int_buffer;
1870	if (!redundant_) {
1871	int redundant_offset = 0;
1872	int nonredundant_offset = 0;
1873	int e12,e13e24,e34;
1874	int i;
1875
1876	if ((osh[0] == osh[3])&&(osh[1] == osh[2])&&(osh[0] != osh[1])) {
1877	ExEnv::errn() << scprintf("nonredundant integrals cannot be generated (1der)\n");
1878	fail();
1879	}
1880
1881	/* Shell equivalence information. */
1882	e12 = (osh[0] == osh[1]);
1883	e13e24 = ((osh[0] == osh[2]) && (osh[1] == osh[3]));
1884	e34 = (osh[2] == osh[3]);
1885	/* Repack x, y, and z integrals. */
1886	for (i=0; i<3; i++) {
1887	nonredundant_erep(current_buffer,e12,e34,e13e24,
1888	shell1->nfunction(),
1889	shell2->nfunction(),
1890	shell3->nfunction(),
1891	shell4->nfunction(),
1892	&redundant_offset,
1893	&nonredundant_offset);
1894	}
1895	}
1896
1897	/* Return the integral buffers to their normal state. */
1898	int_derint_buffer = int_buffer;
1899	int_buffer = user_int_buffer;
1900	}
1901
1902	/* This routine computes a quantity needed to compute the
1903	* derivative integral bounds.
1904	* It fills int_buffer with (sh1+i sh2\|sh1+i sh2).
1905	*/
1906	void
1907	Int2eV3::compute_erep_bound1der(int flags, double *buffer,
1908	int psh1, int psh2, int psh3, int psh4)
1909	{
1910	int oc1,oc2,oc3,oc4;
1911	int ii;
1912	int c1,c2,c3,c4;
1913	int i[4],j[4],k[4],am[4];
1914	int index;
1915	int sizem234,sizem34,sizem2,sizem3,sizem4;
1916	int sizep234,sizep34,sizep2,sizep3,sizep4;
1917	int sizepm234,sizepm34,sizepm2,sizepm3,sizepm4;
1918	GaussianShell shell1,shell2,shell3,shell4;
1919
1920	/* Set up pointers to the current shells. */
1921	shell1 = &bs1_->shell(*psh1);
1922	shell2 = &bs2_->shell(*psh2);
1923	shell3 = &bs3_->shell(*psh3);
1924	shell4 = &bs4_->shell(*psh4);
1925
1926	#define DCT1(n) ((((n)==0)\|\|((n)==2))?1:0)
1927	#define DCT2(n) ((((n)==0)\|\|((n)==2))?((((n)==0))?1:-1):0)
1928	/* Offsets for the intermediates with angular momentum decremented. */
1929	for (ii=sizem2=0; ii<shell2->ncontraction(); ii++)
1930	sizem2 += INT_NCART(shell2->am(ii)-DCT1(1));
1931	for (ii=sizem3=0; ii<shell3->ncontraction(); ii++)
1932	sizem3 += INT_NCART(shell3->am(ii)-DCT1(2));
1933	for (ii=sizem4=0; ii<shell4->ncontraction(); ii++)
1934	sizem4 += INT_NCART(shell4->am(ii)-DCT1(3));
1935	sizem34 = sizem4 * sizem3;
1936	sizem234 = sizem34 * sizem2;
1937
1938	/* Offsets for the intermediates with angular momentum incremented. */
1939	for (ii=sizep2=0; ii<shell2->ncontraction(); ii++)
1940	sizep2 += INT_NCART(shell2->am(ii)+DCT1(1));
1941	for (ii=sizep3=0; ii<shell3->ncontraction(); ii++)
1942	sizep3 += INT_NCART(shell3->am(ii)+DCT1(2));
1943	for (ii=sizep4=0; ii<shell4->ncontraction(); ii++)
1944	sizep4 += INT_NCART(shell4->am(ii)+DCT1(3));
1945	sizep34 = sizep4 * sizep3;
1946	sizep234 = sizep34 * sizep2;
1947
1948	/* Offsets for the intermediates with angular momentum incremented and
1949	* decremented. */
1950	for (ii=sizepm2=0; ii<shell2->ncontraction(); ii++)
1951	sizepm2 += INT_NCART(shell2->am(ii)+DCT2(1));
1952	for (ii=sizepm3=0; ii<shell3->ncontraction(); ii++)
1953	sizepm3 += INT_NCART(shell3->am(ii)+DCT2(2));
1954	for (ii=sizepm4=0; ii<shell4->ncontraction(); ii++)
1955	sizepm4 += INT_NCART(shell4->am(ii)+DCT2(3));
1956	sizepm34 = sizepm4 * sizepm3;
1957	sizepm234 = sizepm34 * sizepm2;
1958
1959	/* END_DERLOOP must be redefined here because it previously depended
1960	* on the DCTEST macro */
1961	#undef END_DERLOOP
1962	#define END_DERLOOP(index,indexp1,sign) \
1963	END_FOR_CART\
1964	oc##indexp1 += INT_NCART(am[index] sign DCT1(index));\
1965	}
1966
1967	#undef END_ALLDERLOOPS
1968	#define END_ALLDERLOOPS(sign)\
1969	END_DERLOOP(3,4,sign)\
1970	END_DERLOOP(2,3,sign)\
1971	END_DERLOOP(1,2,sign)\
1972	END_DERLOOP(0,1,sign)
1973
1974	int old_perm = permute();
1975	set_permute(0);
1976	int old_red = redundant();
1977	set_redundant(1);
1978	compute_erep(flags,psh1,psh2,psh3,psh4,
1979	-DCT1(0),
1980	-DCT1(1),
1981	-DCT1(2),
1982	-DCT1(3));
1983	set_permute(old_perm);
1984	set_redundant(old_red);
1985
1986	/* Place the contributions into the user integral buffer. */
1987	index = 0;
1988	/* The d/dx integrals */
1989	ALLDERLOOPS
1990	if (i[0]>0 && i[2]>0) {
1991	buffer[index] += i[0] * i[2] * int_buffer[
1992	(oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0]-DCT1(0),j[0])) * sizem234
1993	+(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1]-DCT1(1),j[1])) * sizem34
1994	+(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2]-DCT1(2),j[2])) * sizem4
1995	+(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3]-DCT1(3),j[3]))
1996	];
1997	}
1998	index++;
1999	END_ALLDERLOOPS(-)
2000
2001	/* The d/dy integrals */
2002	ALLDERLOOPS
2003	if (j[0]>0 && j[2]>0) {
2004	buffer[index] += j[0] * j[2] * int_buffer[
2005	(oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0],j[0]-DCT1(0))) * sizem234
2006	+(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1],j[1]-DCT1(1))) * sizem34
2007	+(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2],j[2]-DCT1(2))) * sizem4
2008	+(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3],j[3]-DCT1(3)))
2009	];
2010	}
2011	index++;
2012	END_ALLDERLOOPS(-)
2013
2014	/* The d/dz integrals */
2015	ALLDERLOOPS
2016	if (k[0]>0 && k[2]>0) {
2017	buffer[index] += k[0] * k[2] * int_buffer[
2018	(oc1 + INT_CARTINDEX(am[0]-DCT1(0),i[0],j[0])) * sizem234
2019	+(oc2 + INT_CARTINDEX(am[1]-DCT1(1),i[1],j[1])) * sizem34
2020	+(oc3 + INT_CARTINDEX(am[2]-DCT1(2),i[2],j[2])) * sizem4
2021	+(oc4 + INT_CARTINDEX(am[3]-DCT1(3),i[3],j[3]))
2022	];
2023	}
2024	index++;
2025	END_ALLDERLOOPS(-)
2026
2027	/* Compute the next contribution to the integrals. */
2028	/* Tell the build routine that we need an exponent weighted contraction
2029	* with the exponents taken from the dercenter and adjust the
2030	* angular momentum of dercenter to the needed value. */
2031	int_expweight1 = 1;
2032	int_expweight3 = 1;
2033	old_perm = permute();
2034	set_permute(0);
2035	old_red = redundant();
2036	set_redundant(1);
2037	compute_erep(flags,psh1,psh2,psh3,psh4,
2038	DCT1(0),
2039	DCT1(1),
2040	DCT1(2),
2041	DCT1(3));
2042	set_permute(old_perm);
2043	set_redundant(old_red);
2044	int_expweight1 = 0;
2045	int_expweight3 = 0;
2046
2047	/* Place the contributions into the user integral buffer. */
2048	index = 0;
2049	/* The d/dx integrals */
2050	ALLDERLOOPS
2051	buffer[index] += int_buffer[
2052	(oc1+INT_CARTINDEX(am[0]+DCT1(0),i[0]+DCT1(0),j[0]))*sizep234
2053	+(oc2+INT_CARTINDEX(am[1]+DCT1(1),i[1]+DCT1(1),j[1]))*sizep34
2054	+(oc3+INT_CARTINDEX(am[2]+DCT1(2),i[2]+DCT1(2),j[2]))*sizep4
2055	+(oc4+INT_CARTINDEX(am[3]+DCT1(3),i[3]+DCT1(3),j[3]))
2056	];
2057	index++;
2058	END_ALLDERLOOPS(+)
2059
2060	/* The d/dy integrals */
2061	ALLDERLOOPS
2062	buffer[index] += int_buffer[
2063	(oc1+INT_CARTINDEX(am[0]+DCT1(0),i[0],j[0]+DCT1(0)))*sizep234
2064	+(oc2+INT_CARTINDEX(am[1]+DCT1(1),i[1],j[1]+DCT1(1)))*sizep34
2065	+(oc3+INT_CARTINDEX(am[2]+DCT1(2),i[2],j[2]+DCT1(2)))*sizep4
2066	+(oc4+INT_CARTINDEX(am[3]+DCT1(3),i[3],j[3]+DCT1(3)))
2067	];
2068	index++;
2069	END_ALLDERLOOPS(+)
2070
2071	/* The d/dz integrals */
2072	ALLDERLOOPS
2073	buffer[index] += int_buffer[
2074	(oc1 + INT_CARTINDEX(am[0]+DCT1(0),i[0],j[0])) * sizep234
2075	+(oc2 + INT_CARTINDEX(am[1]+DCT1(1),i[1],j[1])) * sizep34
2076	+(oc3 + INT_CARTINDEX(am[2]+DCT1(2),i[2],j[2])) * sizep4
2077	+(oc4 + INT_CARTINDEX(am[3]+DCT1(3),i[3],j[3]))
2078	];
2079	index++;
2080	END_ALLDERLOOPS(+)
2081
2082	/* END_DERLOOP must be redefined here because it previously depended
2083	* on the DCT1 macro */
2084	#undef END_DERLOOP
2085	#define END_DERLOOP(index,indexp1,sign) \
2086	END_FOR_CART\
2087	oc##indexp1 += INT_NCART(am[index] sign DCT2(index));\
2088	}
2089
2090	#undef END_ALLDERLOOPS
2091	#define END_ALLDERLOOPS(sign)\
2092	END_DERLOOP(3,4,sign)\
2093	END_DERLOOP(2,3,sign)\
2094	END_DERLOOP(1,2,sign)\
2095	END_DERLOOP(0,1,sign)
2096
2097	/* Compute the next contribution to the integrals. */
2098	/* Tell the build routine that we need an exponent weighted contraction
2099	* with the exponents taken from the dercenter and adjust the
2100	* angular momentum of dercenter to the needed value. */
2101	int_expweight1 = 1;
2102	old_perm = permute();
2103	set_permute(0);
2104	old_red = redundant();
2105	set_redundant(1);
2106	compute_erep(flags,psh1,psh2,psh3,psh4,
2107	DCT2(0),
2108	DCT2(1),
2109	DCT2(2),
2110	DCT2(3));
2111	set_permute(old_perm);
2112	set_redundant(old_red);
2113	int_expweight1 = 0;
2114
2115	/* Place the contributions into the user integral buffer. */
2116	index = 0;
2117	/* The d/dx integrals */
2118	ALLDERLOOPS
2119	if (i[2] > 0) {
2120	buffer[index] -= 2.0 * int_buffer[
2121	(oc1+INT_CARTINDEX(am[0]+DCT2(0),i[0]+DCT2(0),j[0]))*sizepm234
2122	+(oc2+INT_CARTINDEX(am[1]+DCT2(1),i[1]+DCT2(1),j[1]))*sizepm34
2123	+(oc3+INT_CARTINDEX(am[2]+DCT2(2),i[2]+DCT2(2),j[2]))*sizepm4
2124	+(oc4+INT_CARTINDEX(am[3]+DCT2(3),i[3]+DCT2(3),j[3]))
2125	];
2126	}
2127	index++;
2128	END_ALLDERLOOPS(+)
2129
2130	/* The d/dy integrals */
2131	ALLDERLOOPS
2132	if (j[2] > 0) {
2133	buffer[index] -= 2.0 * int_buffer[
2134	(oc1+INT_CARTINDEX(am[0]+DCT2(0),i[0],j[0]+DCT2(0)))*sizepm234
2135	+(oc2+INT_CARTINDEX(am[1]+DCT2(1),i[1],j[1]+DCT2(1)))*sizepm34
2136	+(oc3+INT_CARTINDEX(am[2]+DCT2(2),i[2],j[2]+DCT2(2)))*sizepm4
2137	+(oc4+INT_CARTINDEX(am[3]+DCT2(3),i[3],j[3]+DCT2(3)))
2138	];
2139	}
2140	index++;
2141	END_ALLDERLOOPS(+)
2142
2143	/* The d/dz integrals */
2144	ALLDERLOOPS
2145	if (k[2] > 0) {
2146	buffer[index] -= 2.0 * int_buffer[
2147	(oc1 + INT_CARTINDEX(am[0]+DCT2(0),i[0],j[0])) * sizepm234
2148	+(oc2 + INT_CARTINDEX(am[1]+DCT2(1),i[1],j[1])) * sizepm34
2149	+(oc3 + INT_CARTINDEX(am[2]+DCT2(2),i[2],j[2])) * sizepm4
2150	+(oc4 + INT_CARTINDEX(am[3]+DCT2(3),i[3],j[3]))
2151	];
2152	}
2153	index++;
2154	END_ALLDERLOOPS(+)
2155	}
2156
2157	/////////////////////////////////////////////////////////////////////////////
2158
2159	// Local Variables:
2160	// mode: c++
2161	// c-file-style: "CLJ-CONDENSED"
2162	// End:

Note: See TracBrowser for help on using the repository browser.

Download in other formats: