Context Navigation

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

Visit:

Candidate_v1.7.0 stable

Last change on this file since 1ca493a was 860145, checked in by Frederik Heber <heber@…>, 9 years ago
Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
Property mode set to `100644`
File size: 65.4 KB

Rev	Line
[0b990d]	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: