Context Navigation

source: ThirdParty/mpqc_open/src/lib/chemistry/qc/intv3/comp1e.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: 85.8 KB

Rev	Line
[0b990d]	1	//
	2	// comp1e.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 <stdlib.h>
	29	#include <math.h>
	30
	31	#include <util/misc/formio.h>
	32	#include <util/misc/math.h>
	33	#include <chemistry/qc/intv3/macros.h>
	34	#include <chemistry/qc/intv3/fjt.h>
	35	#include <chemistry/qc/intv3/utils.h>
	36	#include <chemistry/qc/intv3/int1e.h>
	37	#include <chemistry/qc/intv3/tformv3.h>
	38
	39	using namespace std;
	40	using namespace sc;
	41
	42	#define IN(i,j) ((i)==(j)?1:0)
	43	#define SELECT(x1,x2,x3,s) (((s)==0)?x1:(((s)==1)?(x2):(x3)))
	44
	45	#define DEBUG_NUC_SHELL_DER 0
	46	#define DEBUG_NUC_PRIM 0
	47	#define DEBUG_NUC_SHELL 0
	48
	49	#define DEBUG_EFIELD_PRIM 0
	50
	51	/* ------------ Initialization of 1e routines. ------------------- */
	52	/* This routine returns a buffer large enough to hold a shell doublet
	53	* of integrals (if order == 0) or derivative integrals (if order == 1).
	54	*/
	55	void
	56	Int1eV3::int_initialize_1e(int flags, int order)
	57	{
	58	int jmax1,jmax2,jmax;
	59	int scratchsize=0,nshell2;
	60
	61	/* The efield routines look like derivatives so bump up order if
	62	* it is zero to allow efield integrals to be computed.
	63	*/
	64	if (order == 0) order = 1;
	65
	66	jmax1 = bs1_->max_angular_momentum();
	67	jmax2 = bs2_->max_angular_momentum();
	68	jmax = jmax1 + jmax2;
	69
	70	fjt_ = new FJT(jmax + 2*order);
	71
	72	nshell2 = bs1_->max_ncartesian_in_shell()*bs2_->max_ncartesian_in_shell();
	73
	74	if (order == 0) {
	75	init_order = 0;
	76	scratchsize = nshell2;
	77	}
	78	else if (order == 1) {
	79	init_order = 1;
	80	scratchsize = nshell2*3;
	81	}
	82	else {
	83	ExEnv::errn() << scprintf("int_initialize_1e: invalid order: %d\n",order);
	84	exit(1);
	85	}
	86
	87	buff = new double[scratchsize];
	88	cartesianbuffer = new double[scratchsize];
	89	cartesianbuffer_scratch = new double[scratchsize];
	90
	91	inter.set_dim(jmax1+order+1,jmax2+order+1,jmax+2*order+1);
	92	efield_inter.set_dim(jmax1+order+1,jmax2+order+1,jmax+2*order+1);
	93	int i,j,m;
	94	for (i=0; i<=jmax1+order; i++) {
	95	int sizei = INT_NCART_NN(i);
	96	for (j=0; j<=jmax2+order; j++) {
	97	int sizej = INT_NCART_NN(j);
	98	for (m=0; m<=jmax+2*order-i-j; m++) {
	99	inter(i,j,m) = new double[sizei*sizej];
	100	efield_inter(i,j,m) = new double[sizeisizej3];
	101	}
	102	for (; m<=jmax+2*order; m++) {
	103	inter(i,j,m) = 0;
	104	efield_inter(i,j,m) = 0;
	105	}
	106	}
	107	}
	108
	109	}
	110
	111	void
	112	Int1eV3::int_done_1e()
	113	{
	114	init_order = -1;
	115	delete[] buff;
	116	delete[] cartesianbuffer;
	117	delete[] cartesianbuffer_scratch;
	118	buff = 0;
	119	cartesianbuffer = 0;
	120	inter.delete_data();
	121	efield_inter.delete_data();
	122	}
	123
	124
	125	/* --------------------------------------------------------------- */
	126	/* ------------- Routines for the overlap matrix ----------------- */
	127	/* --------------------------------------------------------------- */
	128
	129	/* This computes the overlap integrals between functions in two shells.
	130	* The result is placed in the buffer.
	131	*/
	132	void
	133	Int1eV3::overlap(int ish, int jsh)
	134	{
	135	int c1,i1,j1,k1,c2,i2,j2,k2;
	136	int gc1,gc2;
	137	int index,index1,index2;
	138
	139	c1 = bs1_->shell_to_center(ish);
	140	c2 = bs2_->shell_to_center(jsh);
	141	for (int xyz=0; xyz<3; xyz++) {
	142	A[xyz] = bs1_->r(c1,xyz);
	143	B[xyz] = bs2_->r(c2,xyz);
	144	}
	145	gshell1 = &bs1_->shell(ish);
	146	gshell2 = &bs2_->shell(jsh);
	147	index = 0;
	148	FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
	149	FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
	150	cartesianbuffer[index] = comp_shell_overlap(gc1,i1,j1,k1,gc2,i2,j2,k2);
	151	index++;
	152	END_FOR_GCCART_GS(index2)
	153	END_FOR_GCCART_GS(index1)
	154
	155	transform_1e(integral_, cartesianbuffer, buff, gshell1, gshell2);
	156	}
	157
	158	/* This computes the overlap ints between functions in two shells.
	159	* The result is placed in the buffer.
	160	*/
	161	void
	162	Int1eV3::overlap_1der(int ish, int jsh,
	163	int idercs, int centernum)
	164	{
	165	int i;
	166	int c1,c2;
	167	int ni,nj;
	168
	169	if (!(init_order >= 0)) {
	170	ExEnv::errn() << scprintf("int_shell_overlap: one electron routines are not init'ed\n");
	171	exit(1);
	172	}
	173
	174	Ref<GaussianBasisSet> dercs;
	175	if (idercs == 0) dercs = bs1_;
	176	else dercs = bs2_;
	177
	178	c1 = bs1_->shell_to_center(ish);
	179	c2 = bs2_->shell_to_center(jsh);
	180	for (int xyz=0; xyz<3; xyz++) {
	181	A[xyz] = bs1_->r(c1,xyz);
	182	B[xyz] = bs2_->r(c2,xyz);
	183	}
	184	gshell1 = &bs1_->shell(ish);
	185	gshell2 = &bs2_->shell(jsh);
	186
	187	ni = gshell1->nfunction();
	188	nj = gshell2->nfunction();
	189
	190	#if 0
	191	ExEnv::outn() << scprintf("zeroing %d%d3 elements of buff\n",ni,nj);
	192	#endif
	193	for (i=0; i<ninj3; i++) {
	194	buff[i] = 0.0;
	195	}
	196
	197	int_accum_shell_overlap_1der(ish,jsh,dercs,centernum);
	198	}
	199
	200	/* This computes the overlap derivative integrals between functions
	201	* in two shells. The result is accumulated in the buffer which is ordered
	202	* atom 0 x, y, z, atom 1, ... .
	203	*/
	204	void
	205	Int1eV3::int_accum_shell_overlap_1der(int ish, int jsh,
	206	Ref<GaussianBasisSet> dercs, int centernum)
	207	{
	208	accum_shell_1der(buff,ish,jsh,dercs,centernum,&Int1eV3::comp_shell_overlap);
	209	}
	210
	211	/* Compute the overlap for the shell. The arguments are the
	212	* cartesian exponents for centers 1 and 2. The shell1 and shell2
	213	* globals are used. */
	214	double
	215	Int1eV3::comp_shell_overlap(int gc1, int i1, int j1, int k1,
	216	int gc2, int i2, int j2, int k2)
	217	{
	218	double exp1,exp2;
	219	int i,j,xyz;
	220	double result;
	221	double Pi;
	222	double oozeta;
	223	double AmB,AmB2;
	224	double tmp;
	225
	226	if ((i1<0)\|\|(j1<0)\|\|(k1<0)\|\|(i2<0)\|\|(j2<0)\|\|(k2<0)) return 0.0;
	227
	228	/* Loop over the primitives in the shells. */
	229	result = 0.0;
	230	for (i=0; i<gshell1->nprimitive(); i++) {
	231	for (j=0; j<gshell2->nprimitive(); j++) {
	232
	233	/* Compute the intermediates. */
	234	exp1 = gshell1->exponent(i);
	235	exp2 = gshell2->exponent(j);
	236	oozeta = 1.0/(exp1 + exp2);
	237	oo2zeta = 0.5*oozeta;
	238	AmB2 = 0.0;
	239	for (xyz=0; xyz<3; xyz++) {
	240	Pi = oozeta(exp1 A[xyz] + exp2 * B[xyz]);
	241	PmA[xyz] = Pi - A[xyz];
	242	PmB[xyz] = Pi - B[xyz];
	243	AmB = A[xyz] - B[xyz];
	244	AmB2 += AmB*AmB;
	245	}
	246	ss = pow(M_PI/(exp1+exp2),1.5)
	247	* exp(- oozeta * exp1 * exp2 * AmB2);
	248	tmp = gshell1->coefficient_unnorm(gc1,i)
	249	* gshell2->coefficient_unnorm(gc2,j)
	250	* comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	251	if (exponent_weighted == 0) tmp *= exp1;
	252	else if (exponent_weighted == 1) tmp *= exp2;
	253	result += tmp;
	254	}
	255	}
	256
	257	return result;
	258	}
	259
	260	/* Compute the overlap between two primitive functions. */
	261	#if 0
	262	double
	263	Int1eV3::int_prim_overlap(shell_t pshell1, shell_t pshell2,
	264	double pA, double pB,
	265	int prim1, int prim2,
	266	int i1, int j1, int k1,
	267	int i2, int j2, int k2)
	268	{
	269	int xyz;
	270	double Pi;
	271	double oozeta;
	272	double AmB,AmB2;
	273
	274	/* Compute the intermediates. */
	275	oozeta = 1.0/(gshell1->exponent(prim1) + gshell2->exponent(prim2));
	276	oo2zeta = 0.5*oozeta;
	277	AmB2 = 0.0;
	278	for (xyz=0; xyz<3; xyz++) {
	279	Pi = oozeta(gshell1->exponent(prim1) A[xyz]
	280	+ gshell2->exponent(prim2) * B[xyz]);
	281	PmA[xyz] = Pi - A[xyz];
	282	PmB[xyz] = Pi - B[xyz];
	283	AmB = A[xyz] - B[xyz];
	284	AmB2 += AmB*AmB;
	285	}
	286	ss = pow(M_PI/(gshell1->exponent(prim1)
	287	+gshell2->exponent(prim2)),1.5)
	288	* exp(- oozeta * gshell1->exponent(prim1)
	289	* gshell2->exponent(prim2) * AmB2);
	290	return comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	291	}
	292	#endif
	293
	294	double
	295	Int1eV3::comp_prim_overlap(int i1, int j1, int k1,
	296	int i2, int j2, int k2)
	297	{
	298	double result;
	299
	300	if (i1) {
	301	result = PmA[0] * comp_prim_overlap(i1-1,j1,k1,i2,j2,k2);
	302	if (i1>1) result += oo2zeta(i1-1) comp_prim_overlap(i1-2,j1,k1,i2,j2,k2);
	303	if (i2>0) result += oo2zetai2 comp_prim_overlap(i1-1,j1,k1,i2-1,j2,k2);
	304	return result;
	305	}
	306	if (j1) {
	307	result = PmA[1] * comp_prim_overlap(i1,j1-1,k1,i2,j2,k2);
	308	if (j1>1) result += oo2zeta(j1-1) comp_prim_overlap(i1,j1-2,k1,i2,j2,k2);
	309	if (j2>0) result += oo2zetaj2 comp_prim_overlap(i1,j1-1,k1,i2,j2-1,k2);
	310	return result;
	311	}
	312	if (k1) {
	313	result = PmA[2] * comp_prim_overlap(i1,j1,k1-1,i2,j2,k2);
	314	if (k1>1) result += oo2zeta(k1-1) comp_prim_overlap(i1,j1,k1-2,i2,j2,k2);
	315	if (k2>0) result += oo2zetak2 comp_prim_overlap(i1,j1,k1-1,i2,j2,k2-1);
	316	return result;
	317	}
	318
	319	if (i2) {
	320	result = PmB[0] * comp_prim_overlap(i1,j1,k1,i2-1,j2,k2);
	321	if (i2>1) result += oo2zeta(i2-1) comp_prim_overlap(i1,j1,k1,i2-2,j2,k2);
	322	if (i1>0) result += oo2zetai1 comp_prim_overlap(i1-1,j1,k1,i2-1,j2,k2);
	323	return result;
	324	}
	325	if (j2) {
	326	result = PmB[1] * comp_prim_overlap(i1,j1,k1,i2,j2-1,k2);
	327	if (j2>1) result += oo2zeta(j2-1) comp_prim_overlap(i1,j1,k1,i2,j2-2,k2);
	328	if (j1>0) result += oo2zetaj1 comp_prim_overlap(i1,j1-1,k1,i2,j2-1,k2);
	329	return result;
	330	}
	331	if (k2) {
	332	result = PmB[2] * comp_prim_overlap(i1,j1,k1,i2,j2,k2-1);
	333	if (k2>1) result += oo2zeta(k2-1) comp_prim_overlap(i1,j1,k1,i2,j2,k2-2);
	334	if (k1>0) result += oo2zetak1 comp_prim_overlap(i1,j1,k1-1,i2,j2,k2-1);
	335	return result;
	336	}
	337
	338	return ss;
	339	}
	340
	341	/* --------------------------------------------------------------- */
	342	/* ------------- Routines for the kinetic energy ----------------- */
	343	/* --------------------------------------------------------------- */
	344
	345	/* This computes the kinetic energy integrals between functions in two shells.
	346	* The result is placed in the buffer.
	347	*/
	348	void
	349	Int1eV3::kinetic(int ish, int jsh)
	350	{
	351	int c1,i1,j1,k1,c2,i2,j2,k2;
	352	int cart1,cart2;
	353	int index;
	354	int gc1,gc2;
	355
	356	c1 = bs1_->shell_to_center(ish);
	357	c2 = bs2_->shell_to_center(jsh);
	358	for (int xyz=0; xyz<3; xyz++) {
	359	A[xyz] = bs1_->r(c1,xyz);
	360	B[xyz] = bs2_->r(c2,xyz);
	361	}
	362	gshell1 = &bs1_->shell(ish);
	363	gshell2 = &bs2_->shell(jsh);
	364	index = 0;
	365	FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
	366	FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
	367	cartesianbuffer[index] = comp_shell_kinetic(gc1,i1,j1,k1,gc2,i2,j2,k2);
	368	index++;
	369	END_FOR_GCCART_GS(cart2)
	370	END_FOR_GCCART_GS(cart1)
	371
	372	transform_1e(integral_, cartesianbuffer, buff, gshell1, gshell2);
	373	}
	374
	375	void
	376	Int1eV3::int_accum_shell_kinetic(int ish, int jsh)
	377	{
	378	int c1,i1,j1,k1,c2,i2,j2,k2;
	379	int cart1,cart2;
	380	int index;
	381	int gc1,gc2;
	382
	383	c1 = bs1_->shell_to_center(ish);
	384	c2 = bs2_->shell_to_center(jsh);
	385	for (int xyz=0; xyz<3; xyz++) {
	386	A[xyz] = bs1_->r(c1,xyz);
	387	B[xyz] = bs2_->r(c2,xyz);
	388	}
	389	gshell1 = &bs1_->shell(ish);
	390	gshell2 = &bs2_->shell(jsh);
	391	index = 0;
	392
	393	FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
	394	FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
	395	cartesianbuffer[index] = comp_shell_kinetic(gc1,i1,j1,k1,gc2,i2,j2,k2);
	396	index++;
	397	END_FOR_GCCART_GS(cart2)
	398	END_FOR_GCCART_GS(cart1)
	399	accum_transform_1e(integral_, cartesianbuffer, buff, gshell1, gshell2);
	400	}
	401
	402	/* This computes the kinetic energy derivative integrals between functions
	403	* in two shells. The result is accumulated in the buffer which is ordered
	404	* atom 0 x, y, z, atom 1, ... .
	405	*/
	406	void
	407	Int1eV3::int_accum_shell_kinetic_1der(int ish, int jsh,
	408	Ref<GaussianBasisSet> dercs, int centernum)
	409	{
	410	accum_shell_1der(buff,ish,jsh,dercs,centernum,&Int1eV3::comp_shell_kinetic);
	411	}
	412
	413	/* This computes the basis function part of
	414	* generic derivative integrals between functions
	415	* in two shells. The result is accumulated in the buffer which is ordered
	416	* atom 0 x, y, z, atom 1, ... .
	417	* The function used to compute the nonderivative integrals is shell_function.
	418	*/
	419	void
	420	Int1eV3::accum_shell_1der(double *buff, int ish, int jsh,
	421	Ref<GaussianBasisSet> dercs, int centernum,
	422	double (Int1eV3::*shell_function)
	423	(int,int,int,int,int,int,int,int))
	424	{
	425	int i;
	426	int gc1,gc2;
	427	int c1,i1,j1,k1,c2,i2,j2,k2;
	428	int index1,index2;
	429	double tmp[3];
	430	double *ctmp = cartesianbuffer;
	431
	432	c1 = bs1_->shell_to_center(ish);
	433	c2 = bs2_->shell_to_center(jsh);
	434	for (int xyz=0; xyz<3; xyz++) {
	435	A[xyz] = bs1_->r(c1,xyz);
	436	B[xyz] = bs2_->r(c2,xyz);
	437	}
	438	gshell1 = &bs1_->shell(ish);
	439	gshell2 = &bs2_->shell(jsh);
	440	FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
	441	FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
	442	if ((bs1_==bs2_)&&(c1==c2)) {
	443	if ( three_center
	444	&& !((bs1_==third_centers)&&(c1==third_centernum))
	445	&& ((bs1_==dercs)&&(c1==centernum))) {
	446	for (i=0; i<3; i++) {
	447	/* Derivative wrt first shell. */
	448	exponent_weighted = 0;
	449	tmp[i] = 2.0 *
	450	(this->*shell_function)(gc1,i1+IN(i,0),j1+IN(i,1),k1+IN(i,2),gc2,i2,j2,k2);
	451	exponent_weighted = -1;
	452	if (SELECT(i1,j1,k1,i)) {
	453	tmp[i] -= SELECT(i1,j1,k1,i) *
	454	(this->*shell_function)(gc1,i1-IN(i,0),j1-IN(i,1),k1-IN(i,2),gc2,i2,j2,k2);
	455	}
	456	/* Derviative wrt second shell. */
	457	exponent_weighted = 1;
	458	tmp[i] += 2.0 *
	459	(this->*shell_function)(gc1,i1,j1,k1,gc2,i2+IN(i,0),j2+IN(i,1),k2+IN(i,2));
	460	exponent_weighted = -1;
	461	if (SELECT(i2,j2,k2,i)) {
	462	tmp[i] -= SELECT(i2,j2,k2,i) *
	463	(this->*shell_function)(gc1,i1,j1,k1,gc2,i2-IN(i,0),j2-IN(i,1),k2-IN(i,2));
	464	}
	465	}
	466	}
	467	else {
	468	/* If there are two centers and they are the same, then we
	469	* use translational invariance to get a net contrib of 0.0 */
	470	for (i=0; i<3; i++) tmp[i] = 0.0;
	471	}
	472	}
	473	else if ((bs1_==dercs)&&(c1==centernum)) {
	474	for (i=0; i<3; i++) {
	475	exponent_weighted = 0;
	476	tmp[i] = 2.0 *
	477	(this->*shell_function)(gc1,i1+IN(i,0),j1+IN(i,1),k1+IN(i,2),gc2,i2,j2,k2);
	478	exponent_weighted = -1;
	479	if (SELECT(i1,j1,k1,i)) {
	480	tmp[i] -= SELECT(i1,j1,k1,i) *
	481	(this->*shell_function)(gc1,i1-IN(i,0),j1-IN(i,1),k1-IN(i,2),gc2,i2,j2,k2);
	482	}
	483	}
	484	}
	485	else if ((bs2_==dercs)&&(c2==centernum)) {
	486	for (i=0; i<3; i++) {
	487	exponent_weighted = 1;
	488	tmp[i] = 2.0 *
	489	(this->*shell_function)(gc1,i1,j1,k1,gc2,i2+IN(i,0),j2+IN(i,1),k2+IN(i,2));
	490	exponent_weighted = -1;
	491	if (SELECT(i2,j2,k2,i)) {
	492	tmp[i] -= SELECT(i2,j2,k2,i) *
	493	(this->*shell_function)(gc1,i1,j1,k1,gc2,i2-IN(i,0),j2-IN(i,1),k2-IN(i,2));
	494	}
	495	}
	496	}
	497	else {
	498	for (i=0; i<3; i++) tmp[i] = 0.0;
	499	}
	500
	501	if (scale_shell_result) {
	502	for (i=0; i<3; i++) tmp[i] *= result_scale_factor;
	503	}
	504
	505	for (i=0; i<3; i++) ctmp[i] = tmp[i];
	506
	507	/* Increment the pointer to the xyz for the next atom. */
	508	ctmp += 3;
	509	END_FOR_GCCART_GS(index2)
	510	END_FOR_GCCART_GS(index1)
	511
	512	accum_transform_1e_xyz(integral_,
	513	cartesianbuffer, buff, gshell1, gshell2);
	514	}
	515
	516	/* This computes the basis function part of
	517	* generic derivative integrals between functions
	518	* in two shells. The result is accumulated in the buffer which is ordered
	519	* atom 0 x, y, z, atom 1, ... .
	520	* The function used to compute the nonderivative integrals is shell_function.
	521	*/
	522	void
	523	Int1eV3::accum_shell_block_1der(double *buff, int ish, int jsh,
	524	Ref<GaussianBasisSet> dercs, int centernum,
	525	void (Int1eV3::*shell_block_function)
	526	(int gc1, int a, int gc2, int b,
	527	int gcsize2, int gcoff1, int gcoff2,
	528	double coef, double *buffer))
	529	{
	530	int i;
	531	int gc1,gc2;
	532	int c1,i1,j1,k1,c2,i2,j2,k2;
	533	int index1,index2;
	534
	535	c1 = bs1_->shell_to_center(ish);
	536	c2 = bs2_->shell_to_center(jsh);
	537
	538	int docenter1=0, docenter2=0;
	539	int equiv12 = (bs1_==bs2_)&&(c1==c2);
	540	int der1 = (bs1_==dercs)&&(c1==centernum);
	541	int der2 = (bs2_==dercs)&&(c2==centernum);
	542	if (!equiv12) {
	543	docenter1 = der1;
	544	docenter2 = der2;
	545	}
	546	else if (three_center) {
	547	int equiv123 = (bs1_==third_centers)&&(c1==third_centernum);
	548	if (!equiv123) {
	549	docenter1 = der1;
	550	docenter2 = der2;
	551	}
	552	}
	553
	554	gshell1 = &bs1_->shell(ish);
	555	gshell2 = &bs2_->shell(jsh);
	556	int gcsize1 = gshell1->ncartesian();
	557	int gcsize2 = gshell2->ncartesian();
	558	memset(cartesianbuffer,0,sizeof(double)gcsize1gcsize2*3);
	559
	560	if (!docenter1 && !docenter2) return;
	561
	562	double coef;
	563	if (scale_shell_result) {
	564	coef = result_scale_factor;
	565	}
	566	else coef = 1.0;
	567
	568	for (int xyz=0; xyz<3; xyz++) {
	569	A[xyz] = bs1_->r(c1,xyz);
	570	B[xyz] = bs2_->r(c2,xyz);
	571	}
	572	int gcoff1 = 0;
	573	for (gc1=0; gc1<gshell1->ncontraction(); gc1++) {
	574	int a = gshell1->am(gc1);
	575	int sizea = INT_NCART_NN(a);
	576	int sizeap1 = INT_NCART_NN(a+1);
	577	int sizeam1 = INT_NCART(a-1);
	578	int gcoff2 = 0;
	579	for (gc2=0; gc2<gshell2->ncontraction(); gc2++) {
	580	int b = gshell2->am(gc2);
	581	int sizeb = INT_NCART_NN(b);
	582	int sizebp1 = INT_NCART_NN(b+1);
	583	int sizebm1 = INT_NCART(b-1);
	584	/* Derivative wrt first shell. */
	585	if (docenter1) {
	586	exponent_weighted = 0;
	587	memset(cartesianbuffer_scratch,0,sizeof(double)sizeap1sizeb);
	588	(this->*shell_block_function)(gc1, a+1, gc2, b,
	589	sizeb, 0, 0,
	590	coef, cartesianbuffer_scratch);
	591	index1=0;
	592	FOR_CART(i1,j1,k1,a) {
	593	index2=0;
	594	FOR_CART(i2,j2,k2,b) {
	595	double ctmp = &cartesianbuffer[((index1+gcoff1)gcsize2
	596	+index2+gcoff2)*3];
	597	for (i=0; i<3; i++) {
	598	int ind = INT_CARTINDEX(a+1,i1+IN(i,0),j1+IN(i,1));
	599	ctmp++ += 2.0cartesianbuffer_scratch[ind*sizeb+index2];
	600	}
	601	index2++;
	602	} END_FOR_CART;
	603	index1++;
	604	} END_FOR_CART;
	605
	606	if (a) {
	607	exponent_weighted = -1;
	608	memset(cartesianbuffer_scratch,0,sizeof(double)sizeam1sizeb);
	609	(this->*shell_block_function)(gc1, a-1, gc2, b,
	610	sizeb, 0, 0,
	611	coef, cartesianbuffer_scratch);
	612	index1=0;
	613	FOR_CART(i1,j1,k1,a) {
	614	index2=0;
	615	FOR_CART(i2,j2,k2,b) {
	616	double ctmp = &cartesianbuffer[((index1+gcoff1)gcsize2
	617	+index2+gcoff2)*3];
	618	for (i=0; i<3; i++) {
	619	int sel = SELECT(i1,j1,k1,i);
	620	if (sel) {
	621	int ind = INT_CARTINDEX(a-1,i1-IN(i,0),j1-IN(i,1));
	622	ctmp[i] -= sel * cartesianbuffer_scratch[ind*sizeb+index2];
	623	}
	624	}
	625	index2++;
	626	} END_FOR_CART;
	627	index1++;
	628	} END_FOR_CART;
	629	}
	630	}
	631	if (docenter2) {
	632	/* Derviative wrt second shell. */
	633	exponent_weighted = 1;
	634	memset(cartesianbuffer_scratch,0,sizeof(double)sizeasizebp1);
	635	(this->*shell_block_function)(gc1, a, gc2, b+1,
	636	sizebp1, 0, 0,
	637	coef, cartesianbuffer_scratch);
	638	index1=0;
	639	FOR_CART(i1,j1,k1,a) {
	640	index2=0;
	641	FOR_CART(i2,j2,k2,b) {
	642	double ctmp = &cartesianbuffer[((index1+gcoff1)gcsize2
	643	+index2+gcoff2)*3];
	644	for (i=0; i<3; i++) {
	645	int ind = INT_CARTINDEX(b+1,i2+IN(i,0),j2+IN(i,1));
	646	ctmp++ += 2.0cartesianbuffer_scratch[index1*sizebp1+ind];
	647	}
	648	index2++;
	649	} END_FOR_CART;
	650	index1++;
	651	} END_FOR_CART;
	652
	653	if (b) {
	654	exponent_weighted = -1;
	655	memset(cartesianbuffer_scratch,0,sizeof(double)sizeasizebm1);
	656	(this->*shell_block_function)(gc1, a, gc2, b-1,
	657	sizebm1, 0, 0,
	658	coef, cartesianbuffer_scratch);
	659	index1=0;
	660	FOR_CART(i1,j1,k1,a) {
	661	index2=0;
	662	FOR_CART(i2,j2,k2,b) {
	663	double ctmp = &cartesianbuffer[((index1+gcoff1)gcsize2
	664	+index2+gcoff2)*3];
	665	for (i=0; i<3; i++) {
	666	int sel = SELECT(i2,j2,k2,i);
	667	if (sel) {
	668	int ind = INT_CARTINDEX(b-1,i2-IN(i,0),j2-IN(i,1));
	669	ctmp[i] -= sel * cartesianbuffer_scratch[index1*sizebm1+ind];
	670	}
	671	}
	672	index2++;
	673	} END_FOR_CART;
	674	index1++;
	675	} END_FOR_CART;
	676	}
	677	}
	678	gcoff2 += INT_NCART_NN(b);
	679	}
	680	gcoff1 += INT_NCART_NN(a);
	681	}
	682
	683	accum_transform_1e_xyz(integral_,
	684	cartesianbuffer, buff, gshell1, gshell2);
	685
	686	#if DEBUG_NUC_SHELL_DER
	687	double *fastbuff = cartesianbuffer;
	688	cartesianbuffer = new double[gcsize1gcsize23];
	689	double junkbuff = new double[gcsize1gcsize2*3];
	690	memset(junkbuff,0,sizeof(double)gcsize1gcsize2*3);
	691	accum_shell_1der(junkbuff,ish,jsh,dercs,centernum,
	692	&Int1eV3::comp_shell_nuclear);
	693	delete[] junkbuff;
	694	int index = 0;
	695	for (i=0; i<gcsize1; i++) {
	696	for (int j=0; j<gcsize2; j++, index++) {
	697	ExEnv::outn() << scprintf(" (%d %d): % 18.15f % 18.15f",
	698	i,j,cartesianbuffer[index],fastbuff[index]);
	699	if (fabs(cartesianbuffer[index]-fastbuff[index])>1.0e-12) {
	700	ExEnv::outn() << " **";
	701	}
	702	ExEnv::outn() << endl;
	703	}
	704	}
	705	delete[] cartesianbuffer;
	706	cartesianbuffer = fastbuff;
	707	#endif
	708	}
	709
	710	/* Compute the kinetic energy for the shell. The arguments are the
	711	* cartesian exponents for centers 1 and 2. The shell1 and shell2
	712	* globals are used. */
	713	double
	714	Int1eV3::comp_shell_kinetic(int gc1, int i1, int j1, int k1,
	715	int gc2, int i2, int j2, int k2)
	716	{
	717	int i,j,xyz;
	718	double result;
	719	double Pi;
	720	double oozeta;
	721	double AmB,AmB2;
	722	double tmp;
	723
	724	/* Loop over the primitives in the shells. */
	725	result = 0.0;
	726	for (i=0; i<gshell1->nprimitive(); i++) {
	727	for (j=0; j<gshell2->nprimitive(); j++) {
	728
	729	/* Compute the intermediates. */
	730	oo2zeta_a = 0.5/gshell1->exponent(i);
	731	oo2zeta_b = 0.5/gshell2->exponent(j);
	732	oozeta = 1.0/(gshell1->exponent(i) + gshell2->exponent(j));
	733	oo2zeta = 0.5*oozeta;
	734	xi = oozeta * gshell1->exponent(i) * gshell2->exponent(j);
	735	AmB2 = 0.0;
	736	for (xyz=0; xyz<3; xyz++) {
	737	Pi = oozeta(gshell1->exponent(i) A[xyz]
	738	+ gshell2->exponent(j) * B[xyz]);
	739	PmA[xyz] = Pi - A[xyz];
	740	PmB[xyz] = Pi - B[xyz];
	741	AmB = A[xyz] - B[xyz];
	742	AmB2 += AmB*AmB;
	743	}
	744	/* The s integral kinetic energy. */
	745	ss = pow(M_PI/(gshell1->exponent(i)
	746	+gshell2->exponent(j)),1.5)
	747	* exp(- xi * AmB2);
	748	sTs = ss
	749	* xi
	750	* (3.0 - 2.0 * xi * AmB2);
	751	tmp = gshell1->coefficient_unnorm(gc1,i)
	752	* gshell2->coefficient_unnorm(gc2,j)
	753	* comp_prim_kinetic(i1,j1,k1,i2,j2,k2);
	754	if (exponent_weighted == 0) tmp *= gshell1->exponent(i);
	755	else if (exponent_weighted == 1) tmp *= gshell2->exponent(j);
	756	result += tmp;
	757	}
	758	}
	759
	760	return result;
	761	}
	762
	763	double
	764	Int1eV3::comp_prim_kinetic(int i1, int j1, int k1,
	765	int i2, int j2, int k2)
	766	{
	767	double tmp;
	768	double result;
	769
	770	if (i1) {
	771	result = PmA[0] * comp_prim_kinetic(i1-1,j1,k1,i2,j2,k2);
	772	if (i1>1) result += oo2zeta(i1-1)comp_prim_kinetic(i1-2,j1,k1,i2,j2,k2);
	773	if (i2) result += oo2zetai2comp_prim_kinetic(i1-1,j1,k1,i2-1,j2,k2);
	774	tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	775	if (i1>1) tmp -= oo2zeta_a(i1-1)comp_prim_overlap(i1-2,j1,k1,i2,j2,k2);
	776	result += 2.0 * xi * tmp;
	777	return result;
	778	}
	779	if (j1) {
	780	result = PmA[1] * comp_prim_kinetic(i1,j1-1,k1,i2,j2,k2);
	781	if (j1>1) result += oo2zeta(j1-1)comp_prim_kinetic(i1,j1-2,k1,i2,j2,k2);
	782	if (j2) result += oo2zetaj2comp_prim_kinetic(i1,j1-1,k1,i2,j2-1,k2);
	783	tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	784	if (j1>1) tmp -= oo2zeta_a(j1-1)comp_prim_overlap(i1,j1-2,k1,i2,j2,k2);
	785	result += 2.0 * xi * tmp;
	786	return result;
	787	}
	788	if (k1) {
	789	result = PmA[2] * comp_prim_kinetic(i1,j1,k1-1,i2,j2,k2);
	790	if (k1>1) result += oo2zeta(k1-1)comp_prim_kinetic(i1,j1,k1-2,i2,j2,k2);
	791	if (k2) result += oo2zetak2comp_prim_kinetic(i1,j1,k1-1,i2,j2,k2-1);
	792	tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	793	if (k1>1) tmp -= oo2zeta_a(k1-1)comp_prim_overlap(i1,j1,k1-2,i2,j2,k2);
	794	result += 2.0 * xi * tmp;
	795	return result;
	796	}
	797	if (i2) {
	798	result = PmB[0] * comp_prim_kinetic(i1,j1,k1,i2-1,j2,k2);
	799	if (i2>1) result += oo2zeta(i2-1)comp_prim_kinetic(i1,j1,k1,i2-2,j2,k2);
	800	if (i1) result += oo2zetai1comp_prim_kinetic(i1-1,j1,k1,i2-1,j2,k2);
	801	tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	802	if (i2>1) tmp -= oo2zeta_b(i2-1)comp_prim_overlap(i1,j1,k1,i2-2,j2,k2);
	803	result += 2.0 * xi * tmp;
	804	return result;
	805	}
	806	if (j2) {
	807	result = PmB[1] * comp_prim_kinetic(i1,j1,k1,i2,j2-1,k2);
	808	if (j2>1) result += oo2zeta(j2-1)comp_prim_kinetic(i1,j1,k1,i2,j2-2,k2);
	809	if (j1) result += oo2zetai1comp_prim_kinetic(i1,j1-1,k1,i2,j2-1,k2);
	810	tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	811	if (j2>1) tmp -= oo2zeta_b(j2-1)comp_prim_overlap(i1,j1,k1,i2,j2-2,k2);
	812	result += 2.0 * xi * tmp;
	813	return result;
	814	}
	815	if (k2) {
	816	result = PmB[2] * comp_prim_kinetic(i1,j1,k1,i2,j2,k2-1);
	817	if (k2>1) result += oo2zeta(k2-1)comp_prim_kinetic(i1,j1,k1,i2,j2,k2-2);
	818	if (k1) result += oo2zetai1comp_prim_kinetic(i1,j1,k1-1,i2,j2,k2-1);
	819	tmp = comp_prim_overlap(i1,j1,k1,i2,j2,k2);
	820	if (k2>1) tmp -= oo2zeta_b(k2-1)comp_prim_overlap(i1,j1,k1,i2,j2,k2-2);
	821	result += 2.0 * xi * tmp;
	822	return result;
	823	}
	824
	825	return sTs;
	826	}
	827
	828	/* --------------------------------------------------------------- */
	829	/* ------------- Routines for the nuclear attraction ------------- */
	830	/* --------------------------------------------------------------- */
	831
	832	/* This computes the nuclear attraction derivative integrals between functions
	833	* in two shells. The result is accumulated in the buffer which is ordered
	834	* atom 0 x, y, z, atom 1, ... .
	835	*/
	836	void
	837	Int1eV3::int_accum_shell_nuclear_1der(int ish, int jsh,
	838	Ref<GaussianBasisSet> dercs, int centernum)
	839	{
	840	int_accum_shell_nuclear_hf_1der(ish,jsh,dercs,centernum);
	841	int_accum_shell_nuclear_nonhf_1der(ish,jsh,dercs,centernum);
	842	}
	843
	844	/* A correction to the Hellman-Feynman part is computed which
	845	* is not included in the original HF routine. This is only needed
	846	* if the real Hellman-Feynman forces are desired, because the sum
	847	* of the hf_1der and nonhf_1der forces are still correct.
	848	*/
	849	void
	850	Int1eV3::int_accum_shell_nuclear_hfc_1der(int ish, int jsh,
	851	Ref<GaussianBasisSet> dercs,
	852	int centernum)
	853	{
	854	/* If both ish and jsh are not on the der center,
	855	* then there's no correction. */
	856	if (!( (bs1_==dercs)
	857	&&(bs2_==dercs)
	858	&&(bs1_->shell_to_center(ish)==centernum)
	859	&&(bs2_->shell_to_center(jsh)==centernum))) {
	860	return;
	861	}
	862
	863	/* Compute the nuc attr part of the nuclear derivative for three equal
	864	* centers. */
	865	scale_shell_result = 1;
	866	result_scale_factor = -bs1_->molecule()->charge(centernum);
	867	for (int xyz=0; xyz<3; xyz++) {
	868	C[xyz] = bs1_->r(centernum,xyz);
	869	}
	870	accum_shell_efield(buff,ish,jsh);
	871	scale_shell_result = 0;
	872
	873	}
	874
	875	/* This computes the nuclear attraction derivative integrals between functions
	876	* in two shells. The result is accumulated in the buffer which is ordered
	877	* atom 0 x, y, z, atom 1, ... . Only the Hellman-Feynman part is computed.
	878	*/
	879	void
	880	Int1eV3::int_accum_shell_nuclear_hf_1der(int ish, int jsh,
	881	Ref<GaussianBasisSet> dercs,
	882	int centernum)
	883	{
	884
	885	/* If both ish and jsh are on the der center, then the contrib is zero. */
	886	if ( (bs1_==dercs)
	887	&&(bs2_==dercs)
	888	&&(bs1_->shell_to_center(ish)==centernum)
	889	&&(bs2_->shell_to_center(jsh)==centernum)) {
	890	return;
	891	}
	892
	893	/* Compute the nuc attr part of the nuclear derivative. */
	894	if (bs1_ == dercs) {
	895	scale_shell_result = 1;
	896	result_scale_factor= -bs1_->molecule()->charge(centernum);
	897	for (int xyz=0; xyz<3; xyz++) {
	898	C[xyz] = bs1_->r(centernum,xyz);
	899	}
	900	//accum_shell_efield(buff,ish,jsh);
	901	accum_shell_block_efield(buff,ish,jsh);
	902	scale_shell_result = 0;
	903	}
	904	else if (bs2_ == dercs) {
	905	scale_shell_result = 1;
	906	result_scale_factor= -bs2_->molecule()->charge(centernum);
	907	for (int xyz=0; xyz<3; xyz++) {
	908	C[xyz] = bs2_->r(centernum,xyz);
	909	}
	910	//accum_shell_efield(buff,ish,jsh);
	911	accum_shell_block_efield(buff,ish,jsh);
	912	scale_shell_result = 0;
	913	}
	914
	915	}
	916
	917	/* This computes the nuclear attraction derivative integrals between functions
	918	* in two shells. The result is accumulated in the buffer which is ordered
	919	* atom 0 x, y, z, atom 1, ... . Only the non Hellman-Feynman part is computed.
	920	*/
	921	void
	922	Int1eV3::int_accum_shell_nuclear_nonhf_1der(int ish, int jsh,
	923	Ref<GaussianBasisSet> dercs,
	924	int centernum)
	925	{
	926	int i;
	927
	928	/* Get the basis function part of the nuclear derivative. */
	929	three_center = 1;
	930	third_centers = bs1_;
	931	for (i=0; i<bs1_->ncenter(); i++) {
	932	third_centernum = i;
	933	for (int xyz=0; xyz<3; xyz++) {
	934	C[xyz] = bs1_->r(i,xyz);
	935	}
	936	scale_shell_result = 1;
	937	result_scale_factor = -bs1_->molecule()->charge(i);
	938	//accum_shell_1der(buff,ish,jsh,dercs,centernum,
	939	// &Int1eV3::comp_shell_nuclear);
	940	accum_shell_block_1der(buff,ish,jsh,dercs,centernum,
	941	&Int1eV3::comp_shell_block_nuclear);
	942	scale_shell_result = 0;
	943	}
	944	if (bs2_!=bs1_) {
	945	third_centers = bs2_;
	946	for (i=0; i<bs2_->ncenter(); i++) {
	947	third_centernum = i;
	948	for (int xyz=0; xyz<3; xyz++) {
	949	C[xyz] = bs2_->r(i,xyz);
	950	}
	951	scale_shell_result = 1;
	952	result_scale_factor = -bs2_->molecule()->charge(i);
	953	//accum_shell_1der(buff,ish,jsh,dercs,centernum,
	954	// &Int1eV3::comp_shell_nuclear);
	955	accum_shell_block_1der(buff,ish,jsh,dercs,centernum,
	956	&Int1eV3::comp_shell_block_nuclear);
	957	scale_shell_result = 0;
	958	}
	959	}
	960	three_center = 0;
	961
	962	}
	963
	964	/* This computes the efield integrals between functions in two shells.
	965	* The result is accumulated in the buffer in the form bf1 x y z, bf2
	966	* x y z, etc.
	967	*/
	968	void
	969	Int1eV3::int_accum_shell_efield(int ish, int jsh,
	970	double *position)
	971	{
	972	scale_shell_result = 0;
	973	for (int xyz=0; xyz<3; xyz++) {
	974	C[xyz] = position[xyz];
	975	}
	976	accum_shell_efield(buff,ish,jsh);
	977	}
	978
	979	/* This computes the efield integrals between functions in two shells.
	980	* The result is accumulated in the buffer in the form bf1 x y z, bf2
	981	* x y z, etc. The globals scale_shell_result, result_scale_factor,
	982	* and C must be set before this is called.
	983	*/
	984	void
	985	Int1eV3::accum_shell_efield(double *buff, int ish, int jsh)
	986	{
	987	int i;
	988	int c1,i1,j1,k1,c2,i2,j2,k2;
	989	double efield[3];
	990	int gc1,gc2;
	991	int index1,index2;
	992	double *tmp = cartesianbuffer;
	993
	994	if (!(init_order >= 1)) {
	995	ExEnv::errn() << scprintf("accum_shell_efield: one electron routines are not init'ed\n");
	996	exit(1);
	997	}
	998
	999	c1 = bs1_->shell_to_center(ish);
	1000	c2 = bs2_->shell_to_center(jsh);
	1001	for (int xyz=0; xyz<3; xyz++) {
	1002	A[xyz] = bs1_->r(c1,xyz);
	1003	B[xyz] = bs2_->r(c2,xyz);
	1004	}
	1005	gshell1 = &bs1_->shell(ish);
	1006	gshell2 = &bs2_->shell(jsh);
	1007
	1008	FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
	1009	FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
	1010	comp_shell_efield(efield,gc1,i1,j1,k1,gc2,i2,j2,k2);
	1011	if (scale_shell_result) {
	1012	for (i=0; i<3; i++) efield[i] *= result_scale_factor;
	1013	}
	1014	for (i=0; i<3; i++) tmp[i] = efield[i];
	1015	tmp += 3;
	1016	END_FOR_GCCART_GS(index2)
	1017	END_FOR_GCCART_GS(index1)
	1018
	1019	accum_transform_1e_xyz(integral_,
	1020	cartesianbuffer, buff, gshell1, gshell2);
	1021	}
	1022
	1023	/* This computes the efield integrals between functions in two shells.
	1024	* The result is accumulated in the buffer in the form bf1 x y z, bf2
	1025	* x y z, etc. The globals scale_shell_result, result_scale_factor,
	1026	* and C must be set before this is called.
	1027	*/
	1028	void
	1029	Int1eV3::accum_shell_block_efield(double *buff, int ish, int jsh)
	1030	{
	1031	int c1,c2;
	1032	int gc1,gc2;
	1033
	1034	if (!(init_order >= 1)) {
	1035	ExEnv::errn() << scprintf("accum_shell_efield: one electron routines are not init'ed\n");
	1036	exit(1);
	1037	}
	1038
	1039	c1 = bs1_->shell_to_center(ish);
	1040	c2 = bs2_->shell_to_center(jsh);
	1041	for (int xyz=0; xyz<3; xyz++) {
	1042	A[xyz] = bs1_->r(c1,xyz);
	1043	B[xyz] = bs2_->r(c2,xyz);
	1044	}
	1045	gshell1 = &bs1_->shell(ish);
	1046	gshell2 = &bs2_->shell(jsh);
	1047
	1048	double coef;
	1049	if (scale_shell_result) coef = result_scale_factor;
	1050	else coef = 1.0;
	1051
	1052	int gcsize1 = gshell1->ncartesian();
	1053	int gcsize2 = gshell2->ncartesian();
	1054	memset(cartesianbuffer,0,sizeof(double)gcsize1gcsize2*3);
	1055	int gcoff1 = 0;
	1056	for (gc1=0; gc1<gshell1->ncontraction(); gc1++) {
	1057	int a = gshell1->am(gc1);
	1058	int sizea = INT_NCART_NN(a);
	1059	int gcoff2 = 0;
	1060	for (gc2=0; gc2<gshell2->ncontraction(); gc2++) {
	1061	int b = gshell2->am(gc2);
	1062	int sizeb = INT_NCART_NN(b);
	1063	comp_shell_block_efield(gc1,a,gc2,b,
	1064	gcsize2, gcoff1, gcoff2,
	1065	coef, cartesianbuffer);
	1066	gcoff2 += sizeb;
	1067	}
	1068	gcoff1 += sizea;
	1069	}
	1070
	1071	accum_transform_1e_xyz(integral_,
	1072	cartesianbuffer, buff, gshell1, gshell2);
	1073	}
	1074
	1075	/* This computes the efield integrals between functions in two shells.
	1076	* The result is placed in the buffer in the form bf1 x y z, bf2
	1077	* x y z, etc.
	1078	*/
	1079	void
	1080	Int1eV3::efield(int ish, int jsh, double *position)
	1081	{
	1082	scale_shell_result = 0;
	1083	int xyz;
	1084
	1085	for (xyz=0; xyz<3; xyz++) {
	1086	C[xyz] = position[xyz];
	1087	}
	1088
	1089	int i;
	1090	int c1,i1,j1,k1,c2,i2,j2,k2;
	1091	double efield[3];
	1092	int gc1,gc2;
	1093	int index1,index2;
	1094	double *tmp = cartesianbuffer;
	1095
	1096	if (!(init_order >= 1)) {
	1097	ExEnv::errn() << scprintf("Int1eV3::efield one electron routines are not ready\n");
	1098	exit(1);
	1099	}
	1100
	1101	c1 = bs1_->shell_to_center(ish);
	1102	c2 = bs2_->shell_to_center(jsh);
	1103	for (xyz=0; xyz<3; xyz++) {
	1104	A[xyz] = bs1_->r(c1,xyz);
	1105	B[xyz] = bs2_->r(c2,xyz);
	1106	}
	1107	gshell1 = &bs1_->shell(ish);
	1108	gshell2 = &bs2_->shell(jsh);
	1109
	1110	FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
	1111	FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
	1112	comp_shell_efield(efield,gc1,i1,j1,k1,gc2,i2,j2,k2);
	1113	if (scale_shell_result) {
	1114	for (i=0; i<3; i++) efield[i] *= result_scale_factor;
	1115	}
	1116	for (i=0; i<3; i++) tmp[i] = efield[i];
	1117	tmp += 3;
	1118	END_FOR_GCCART_GS(index2)
	1119	END_FOR_GCCART_GS(index1)
	1120
	1121	transform_1e_xyz(integral_,
	1122	cartesianbuffer, buff, gshell1, gshell2);
	1123	}
	1124
	1125	/* This slowly computes the nuc rep energy integrals between functions in
	1126	* two shells. The result is placed in the buffer. */
	1127	void
	1128	Int1eV3::nuclear_slow(int ish, int jsh)
	1129	{
	1130	int i;
	1131	int c1,i1,j1,k1,c2,i2,j2,k2;
	1132	int index;
	1133	int gc1,gc2;
	1134	int cart1,cart2;
	1135
	1136	if (!(init_order >= 0)) {
	1137	ExEnv::errn() << scprintf("int_shell_nuclear: one electron routines are not init'ed\n");
	1138	exit(1);
	1139	}
	1140
	1141	c1 = bs1_->shell_to_center(ish);
	1142	c2 = bs2_->shell_to_center(jsh);
	1143	for (int xyz=0; xyz<3; xyz++) {
	1144	A[xyz] = bs1_->r(c1,xyz);
	1145	B[xyz] = bs2_->r(c2,xyz);
	1146	}
	1147	gshell1 = &bs1_->shell(ish);
	1148	gshell2 = &bs2_->shell(jsh);
	1149	index = 0;
	1150
	1151	FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
	1152	FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
	1153	cartesianbuffer[index] = 0.0;
	1154	/* Loop thru the centers on bs1_. */
	1155	for (i=0; i<bs1_->ncenter(); i++) {
	1156	for (int xyz=0; xyz<3; xyz++) {
	1157	C[xyz] = bs1_->r(i,xyz);
	1158	}
	1159	cartesianbuffer[index] -= comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
	1160	* bs1_->molecule()->charge(i);
	1161	}
	1162	/* Loop thru the centers on bs2_ if necessary. */
	1163	if (bs2_ != bs1_) {
	1164	for (i=0; i<bs2_->ncenter(); i++) {
	1165	for (int xyz=0; xyz<3; xyz++) {
	1166	C[xyz] = bs2_->r(i,xyz);
	1167	}
	1168	cartesianbuffer[index]-=comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
	1169	* bs2_->molecule()->charge(i);
	1170	}
	1171	}
	1172	index++;
	1173	END_FOR_GCCART_GS(cart2)
	1174	END_FOR_GCCART_GS(cart1)
	1175
	1176	transform_1e(integral_,
	1177	cartesianbuffer, buff, gshell1, gshell2);
	1178	}
	1179
	1180	/* This computes the nuc rep energy integrals between functions in two
	1181	* shells. The result is placed in the buffer. */
	1182	void
	1183	Int1eV3::nuclear(int ish, int jsh)
	1184	{
	1185	int i;
	1186	int c1,c2;
	1187	int gc1,gc2;
	1188
	1189	if (!(init_order >= 0)) {
	1190	ExEnv::errn() << scprintf("int_shell_nuclear: one electron routines are not init'ed\n");
	1191	exit(1);
	1192	}
	1193
	1194	c1 = bs1_->shell_to_center(ish);
	1195	c2 = bs2_->shell_to_center(jsh);
	1196	for (int xyz=0; xyz<3; xyz++) {
	1197	A[xyz] = bs1_->r(c1,xyz);
	1198	B[xyz] = bs2_->r(c2,xyz);
	1199	}
	1200	gshell1 = &bs1_->shell(ish);
	1201	gshell2 = &bs2_->shell(jsh);
	1202
	1203	int ni = gshell1->ncartesian();
	1204	int nj = gshell2->ncartesian();
	1205	memset(cartesianbuffer,0,sizeof(double)ninj);
	1206
	1207	int offi = 0;
	1208	for (gc1=0; gc1<gshell1->ncontraction(); gc1++) {
	1209	int a = gshell1->am(gc1);
	1210	int offj = 0;
	1211	for (gc2=0; gc2<gshell2->ncontraction(); gc2++) {
	1212	int b = gshell2->am(gc2);
	1213	/* Loop thru the centers on bs1_. */
	1214	for (i=0; i<bs1_->ncenter(); i++) {
	1215	double charge = bs1_->molecule()->charge(i);
	1216	for (int xyz=0; xyz<3; xyz++) C[xyz] = bs1_->r(i,xyz);
	1217	comp_shell_block_nuclear(gc1, a, gc2, b,
	1218	nj, offi, offj,
	1219	-charge, cartesianbuffer);
	1220	}
	1221	/* Loop thru the centers on bs2_ if necessary. */
	1222	if (bs2_ != bs1_) {
	1223	for (i=0; i<bs2_->ncenter(); i++) {
	1224	double charge = bs2_->molecule()->charge(i);
	1225	for (int xyz=0; xyz<3; xyz++) C[xyz] = bs2_->r(i,xyz);
	1226	comp_shell_block_nuclear(gc1, a, gc2, b,
	1227	nj, offi, offj,
	1228	-charge, cartesianbuffer);
	1229	}
	1230	}
	1231	offj += INT_NCART_NN(b);
	1232	}
	1233	offi += INT_NCART_NN(a);
	1234	}
	1235
	1236	#if DEBUG_NUC_SHELL
	1237	double fastbuf = new double[ninj];
	1238	memcpy(fastbuf,cartesianbuffer,sizeof(double)ninj);
	1239	nuclear_slow(ish,jsh);
	1240
	1241	index = 0;
	1242	int cart1, cart2;
	1243	FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1) {
	1244	FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2) {
	1245	double fast = fastbuf[index];
	1246	double slow = cartesianbuffer[index];
	1247	if (fabs(fast-slow)>1.0e-12) {
	1248	ExEnv::outn() << scprintf("NUC SHELL FINAL: %d (%d %d %d) %d (%d %d %d): ",
	1249	gc1, i1,j1,k1, gc2, i2,j2,k2)
	1250	<< scprintf(" % 20.15f % 20.15f",
	1251	fast, slow)
	1252	<< endl;
	1253	}
	1254	index++;
	1255	} END_FOR_GCCART_GS(cart2);
	1256	} END_FOR_GCCART_GS(cart1);
	1257
	1258	memcpy(cartesianbuffer,fastbuf,sizeof(double)ninj);
	1259	delete[] fastbuf;
	1260	#endif
	1261
	1262	transform_1e(integral_,
	1263	cartesianbuffer, buff, gshell1, gshell2);
	1264	}
	1265
	1266	/* This computes the integrals between functions in two shells for
	1267	* a point charge interaction operator.
	1268	* The result is placed in the buffer.
	1269	*/
	1270	void
	1271	Int1eV3::int_accum_shell_point_charge(int ish, int jsh,
	1272	int ncharge, const double* charge,
	1273	const doubleconst position)
	1274	{
	1275	int i;
	1276	int c1,i1,j1,k1,c2,i2,j2,k2;
	1277	int index;
	1278	int gc1,gc2;
	1279	int cart1,cart2;
	1280	double tmp;
	1281
	1282	if (!(init_order >= 0)) {
	1283	ExEnv::errn() << scprintf("int_shell_pointcharge: one electron routines are not init'ed\n");
	1284	exit(1);
	1285	}
	1286
	1287	c1 = bs1_->shell_to_center(ish);
	1288	c2 = bs2_->shell_to_center(jsh);
	1289	for (int xyz=0; xyz<3; xyz++) {
	1290	A[xyz] = bs1_->r(c1,xyz);
	1291	B[xyz] = bs2_->r(c2,xyz);
	1292	}
	1293	gshell1 = &bs1_->shell(ish);
	1294	gshell2 = &bs2_->shell(jsh);
	1295	index = 0;
	1296
	1297	FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
	1298	FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
	1299	/* Loop thru the point charges. */
	1300	tmp = 0.0;
	1301	for (i=0; i<ncharge; i++) {
	1302	for (int xyz=0; xyz<3; xyz++) {
	1303	C[xyz] = position[i][xyz];
	1304	}
	1305	tmp -= comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
	1306	* charge[i];
	1307	}
	1308	cartesianbuffer[index] = tmp;
	1309	index++;
	1310	END_FOR_GCCART_GS(cart2)
	1311	END_FOR_GCCART_GS(cart1)
	1312
	1313	accum_transform_1e(integral_,
	1314	cartesianbuffer, buff, gshell1, gshell2);
	1315	}
	1316
	1317	/* This computes the integrals between functions in two shells for
	1318	* a point charge interaction operator.
	1319	* The result is placed in the buffer.
	1320	*/
	1321	void
	1322	Int1eV3::point_charge(int ish, int jsh,
	1323	int ncharge,
	1324	const double* charge, const doubleconst position)
	1325	{
	1326	int i;
	1327	int c1,i1,j1,k1,c2,i2,j2,k2;
	1328	int index;
	1329	int gc1,gc2;
	1330	int cart1,cart2;
	1331
	1332	if (!(init_order >= 0)) {
	1333	ExEnv::errn() << scprintf("Int1eV3::point_charge: one electron routines are not init'ed\n");
	1334	exit(1);
	1335	}
	1336
	1337	c1 = bs1_->shell_to_center(ish);
	1338	c2 = bs2_->shell_to_center(jsh);
	1339	for (int xyz=0; xyz<3; xyz++) {
	1340	A[xyz] = bs1_->r(c1,xyz);
	1341	B[xyz] = bs2_->r(c2,xyz);
	1342	}
	1343	gshell1 = &bs1_->shell(ish);
	1344	gshell2 = &bs2_->shell(jsh);
	1345	index = 0;
	1346
	1347	FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
	1348	FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
	1349	cartesianbuffer[index] = 0.0;
	1350	/* Loop thru the point charges. */
	1351	for (i=0; i<ncharge; i++) {
	1352	for (int xyz=0; xyz<3; xyz++) {
	1353	C[xyz] = position[i][xyz];
	1354	}
	1355	cartesianbuffer[index] -= comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
	1356	* charge[i];
	1357	}
	1358	index++;
	1359	END_FOR_GCCART_GS(cart2)
	1360	END_FOR_GCCART_GS(cart1)
	1361
	1362	transform_1e(integral_,
	1363	cartesianbuffer, buff, gshell1, gshell2);
	1364	}
	1365
	1366
	1367	/* This computes the 1e Hamiltonian integrals between functions in two shells.
	1368	* The result is placed in the buffer.
	1369	*/
	1370	void
	1371	Int1eV3::hcore(int ish, int jsh)
	1372	{
	1373	int i;
	1374	int c1,i1,j1,k1,c2,i2,j2,k2;
	1375	int index;
	1376	int cart1,cart2;
	1377	int gc1,gc2;
	1378
	1379	if (!(init_order >= 0)) {
	1380	ExEnv::errn() << scprintf("hcore: one electron routines are not init'ed\n");
	1381	exit(1);
	1382	}
	1383
	1384	c1 = bs1_->shell_to_center(ish);
	1385	c2 = bs2_->shell_to_center(jsh);
	1386	for (int xyz=0; xyz<3; xyz++) {
	1387	A[xyz] = bs1_->r(c1,xyz);
	1388	B[xyz] = bs2_->r(c2,xyz);
	1389	}
	1390	gshell1 = &bs1_->shell(ish);
	1391	gshell2 = &bs2_->shell(jsh);
	1392
	1393	index = 0;
	1394	FOR_GCCART_GS(gc1,cart1,i1,j1,k1,gshell1)
	1395	FOR_GCCART_GS(gc2,cart2,i2,j2,k2,gshell2)
	1396	cartesianbuffer[index] = comp_shell_kinetic(gc1,i1,j1,k1,gc2,i2,j2,k2);
	1397	/* Loop thru the centers on bs1_. */
	1398	for (i=0; i<bs1_->ncenter(); i++) {
	1399	for (int xyz=0; xyz<3; xyz++) {
	1400	C[xyz] = bs1_->r(i,xyz);
	1401	}
	1402	cartesianbuffer[index] -= comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
	1403	* bs1_->molecule()->charge(i);
	1404	}
	1405	/* Loop thru the centers on bs2_ if necessary. */
	1406	if (bs2_ != bs1_) {
	1407	for (i=0; i<bs2_->ncenter(); i++) {
	1408	for (int xyz=0; xyz<3; xyz++) {
	1409	C[xyz] = bs2_->r(i,xyz);
	1410	}
	1411	cartesianbuffer[index]-=comp_shell_nuclear(gc1,i1,j1,k1,gc2,i2,j2,k2)
	1412	* bs2_->molecule()->charge(i);
	1413	}
	1414	}
	1415	index++;
	1416	END_FOR_GCCART_GS(cart2)
	1417	END_FOR_GCCART_GS(cart1)
	1418
	1419	transform_1e(integral_,
	1420	cartesianbuffer, buff, gshell1, gshell2);
	1421	}
	1422
	1423	/* This computes the 1e Hamiltonian deriv ints between functions in two shells.
	1424	* The result is placed in the buffer.
	1425	*/
	1426	void
	1427	Int1eV3::hcore_1der(int ish, int jsh,
	1428	int idercs, int centernum)
	1429	{
	1430	int i;
	1431	int c1,c2;
	1432	int ni,nj;
	1433
	1434	if (!(init_order >= 0)) {
	1435	ExEnv::errn() << scprintf("int_shell_hcore: one electron routines are not init'ed\n");
	1436	exit(1);
	1437	}
	1438
	1439	Ref<GaussianBasisSet> dercs;
	1440	if (idercs == 0) dercs = bs1_;
	1441	else dercs = bs2_;
	1442
	1443	c1 = bs1_->shell_to_center(ish);
	1444	c2 = bs2_->shell_to_center(jsh);
	1445	for (int xyz=0; xyz<3; xyz++) {
	1446	A[xyz] = bs1_->r(c1,xyz);
	1447	B[xyz] = bs2_->r(c2,xyz);
	1448	}
	1449	gshell1 = &bs1_->shell(ish);
	1450	gshell2 = &bs2_->shell(jsh);
	1451
	1452	ni = gshell1->nfunction();
	1453	nj = gshell2->nfunction();
	1454
	1455	for (i=0; i<ninj3; i++) {
	1456	buff[i] = 0.0;
	1457	}
	1458
	1459	int_accum_shell_nuclear_1der(ish,jsh,dercs,centernum);
	1460	int_accum_shell_kinetic_1der(ish,jsh,dercs,centernum);
	1461	}
	1462
	1463	/* This computes the kinetic deriv ints between functions in two shells.
	1464	* The result is placed in the buffer.
	1465	*/
	1466	void
	1467	Int1eV3::kinetic_1der(int ish, int jsh,
	1468	int idercs, int centernum)
	1469	{
	1470	int i;
	1471	int c1,c2;
	1472	int ni,nj;
	1473
	1474	if (!(init_order >= 0)) {
	1475	ExEnv::errn() << scprintf("int_shell_kinetic: one electron routines are not init'ed\n");
	1476	exit(1);
	1477	}
	1478
	1479	Ref<GaussianBasisSet> dercs;
	1480	if (idercs == 0) dercs = bs1_;
	1481	else dercs = bs2_;
	1482
	1483	c1 = bs1_->shell_to_center(ish);
	1484	c2 = bs2_->shell_to_center(jsh);
	1485	for (int xyz=0; xyz<3; xyz++) {
	1486	A[xyz] = bs1_->r(c1,xyz);
	1487	B[xyz] = bs2_->r(c2,xyz);
	1488	}
	1489	gshell1 = &bs1_->shell(ish);
	1490	gshell2 = &bs2_->shell(jsh);
	1491
	1492	ni = gshell1->nfunction();
	1493	nj = gshell2->nfunction();
	1494
	1495	for (i=0; i<ninj3; i++) {
	1496	buff[i] = 0.0;
	1497	}
	1498
	1499	int_accum_shell_kinetic_1der(ish,jsh,dercs,centernum);
	1500	}
	1501
	1502	/* This computes the nuclear deriv ints between functions in two shells.
	1503	* The result is placed in the buffer.
	1504	*/
	1505	void
	1506	Int1eV3::nuclear_1der(int ish, int jsh, int idercs, int centernum)
	1507	{
	1508	int i;
	1509	int c1,c2;
	1510	int ni,nj;
	1511
	1512	if (!(init_order >= 0)) {
	1513	ExEnv::errn() << scprintf("int_shell_nuclear: one electron routines are not init'ed\n");
	1514	exit(1);
	1515	}
	1516
	1517	Ref<GaussianBasisSet> dercs;
	1518	if (idercs == 0) dercs = bs1_;
	1519	else dercs = bs2_;
	1520
	1521	c1 = bs1_->shell_to_center(ish);
	1522	c2 = bs2_->shell_to_center(jsh);
	1523	for (int xyz=0; xyz<3; xyz++) {
	1524	A[xyz] = bs1_->r(c1,xyz);
	1525	B[xyz] = bs2_->r(c2,xyz);
	1526	}
	1527	gshell1 = &bs1_->shell(ish);
	1528	gshell2 = &bs2_->shell(jsh);
	1529
	1530	ni = gshell1->nfunction();
	1531	nj = gshell2->nfunction();
	1532
	1533	for (i=0; i<ninj3; i++) {
	1534	buff[i] = 0.0;
	1535	}
	1536
	1537	int_accum_shell_nuclear_1der(ish,jsh,dercs,centernum);
	1538	}
	1539
	1540	/* This computes the nuclear deriv ints between functions in two shells.
	1541	* Only the Hellman-Feynman part is computed.
	1542	* The result is placed in the buffer.
	1543	*/
	1544	void
	1545	Int1eV3::int_shell_nuclear_hf_1der(int ish, int jsh,
	1546	Ref<GaussianBasisSet> dercs, int centernum)
	1547	{
	1548	int i;
	1549	int c1,c2;
	1550	int ni,nj;
	1551
	1552	if (!(init_order >= 0)) {
	1553	ExEnv::errn() << scprintf("int_shell_nuclear_hf_1der: one electron routines are not init'ed\n");
	1554	exit(1);
	1555	}
	1556
	1557	c1 = bs1_->shell_to_center(ish);
	1558	c2 = bs2_->shell_to_center(jsh);
	1559	for (int xyz=0; xyz<3; xyz++) {
	1560	A[xyz] = bs1_->r(c1,xyz);
	1561	B[xyz] = bs2_->r(c2,xyz);
	1562	}
	1563	gshell1 = &bs1_->shell(ish);
	1564	gshell2 = &bs2_->shell(jsh);
	1565
	1566	ni = gshell1->nfunction();
	1567	nj = gshell2->nfunction();
	1568
	1569	for (i=0; i<ninj3; i++) {
	1570	buff[i] = 0.0;
	1571	}
	1572
	1573	int_accum_shell_nuclear_hf_1der(ish,jsh,dercs,centernum);
	1574	}
	1575
	1576	/* This computes the nuclear deriv ints between functions in two shells.
	1577	* Only the non Hellman-Feynman part is computed.
	1578	* The result is placed in the buffer.
	1579	*/
	1580	void
	1581	Int1eV3::int_shell_nuclear_nonhf_1der(int ish, int jsh,
	1582	Ref<GaussianBasisSet> dercs, int centernum)
	1583	{
	1584	int i;
	1585	int c1,c2;
	1586	int ni,nj;
	1587
	1588	if (!(init_order >= 0)) {
	1589	ExEnv::errn() << scprintf("int_shell_nuclear_nonhf_1der: one electron routines are not init'ed\n");
	1590	exit(1);
	1591	}
	1592
	1593	c1 = bs1_->shell_to_center(ish);
	1594	c2 = bs2_->shell_to_center(jsh);
	1595	for (int xyz=0; xyz<3; xyz++) {
	1596	A[xyz] = bs1_->r(c1,xyz);
	1597	B[xyz] = bs2_->r(c2,xyz);
	1598	}
	1599	gshell1 = &bs1_->shell(ish);
	1600	gshell2 = &bs2_->shell(jsh);
	1601
	1602	ni = gshell1->nfunction();
	1603	nj = gshell2->nfunction();
	1604
	1605	#if 0
	1606	ExEnv::outn() << scprintf("int_shell_nuclear_nonhf_1der: zeroing %d doubles in buff\n",ninj3);
	1607	#endif
	1608	for (i=0; i<ninj3; i++) {
	1609	buff[i] = 0.0;
	1610	}
	1611
	1612	int_accum_shell_nuclear_nonhf_1der(ish,jsh,dercs,centernum);
	1613	}
	1614
	1615	/* Compute the nuclear attraction for the shell. The arguments are the
	1616	* cartesian exponents for centers 1 and 2. The shell1 and shell2
	1617	* globals are used. */
	1618	double
	1619	Int1eV3::comp_shell_nuclear(int gc1, int i1, int j1, int k1,
	1620	int gc2, int i2, int j2, int k2)
	1621	{
	1622	int i,j,k,xyz;
	1623	double result;
	1624	double Pi;
	1625	double oozeta;
	1626	double AmB,AmB2;
	1627	double PmC2;
	1628	double auxcoef;
	1629	int am;
	1630	double tmp;
	1631
	1632	am = i1+j1+k1+i2+j2+k2;
	1633
	1634	/* Loop over the primitives in the shells. */
	1635	result = 0.0;
	1636	for (i=0; i<gshell1->nprimitive(); i++) {
	1637	for (j=0; j<gshell2->nprimitive(); j++) {
	1638
	1639	/* Compute the intermediates. */
	1640	zeta = gshell1->exponent(i) + gshell2->exponent(j);
	1641	oozeta = 1.0/zeta;
	1642	oo2zeta = 0.5*oozeta;
	1643	AmB2 = 0.0;
	1644	PmC2 = 0.0;
	1645	for (xyz=0; xyz<3; xyz++) {
	1646	Pi = oozeta(gshell1->exponent(i) A[xyz]
	1647	+ gshell2->exponent(j) * B[xyz]);
	1648	PmA[xyz] = Pi - A[xyz];
	1649	PmB[xyz] = Pi - B[xyz];
	1650	PmC[xyz] = Pi - C[xyz];
	1651	AmB = A[xyz] - B[xyz];
	1652	AmB2 += AmB*AmB;
	1653	PmC2 += PmC[xyz]*PmC[xyz];
	1654	}
	1655
	1656	/* The auxillary integral coeficients. */
	1657	auxcoef = 2.0 * M_PI/(gshell1->exponent(i)
	1658	+gshell2->exponent(j))
	1659	* exp(- oozeta * gshell1->exponent(i)
	1660	* gshell2->exponent(j) * AmB2);
	1661
	1662	/* The Fm(U) intermediates. */
	1663	fjttable_ = fjt_->values(am,zeta*PmC2);
	1664
	1665	/* Convert the Fm(U) intermediates into the auxillary
	1666	* nuclear attraction integrals. */
	1667	for (k=0; k<=am; k++) {
	1668	fjttable_[k] *= auxcoef;
	1669	}
	1670
	1671	/* Compute the nuclear attraction integral. */
	1672	tmp = gshell1->coefficient_unnorm(gc1,i)
	1673	* gshell2->coefficient_unnorm(gc2,j)
	1674	* comp_prim_nuclear(i1,j1,k1,i2,j2,k2,0);
	1675
	1676	if (exponent_weighted == 0) tmp *= gshell1->exponent(i);
	1677	else if (exponent_weighted == 1) tmp *= gshell2->exponent(j);
	1678
	1679	result += tmp;
	1680	}
	1681	}
	1682
	1683	/* printf("comp_shell_nuclear(%d,%d,%d,%d,%d,%d): result = % 12.8lf\n",
	1684	* i1,j1,k1,i2,j2,k2,result);
	1685	*/
	1686	return result;
	1687	}
	1688
	1689	/* Compute the nuclear attraction for the shell. The arguments are the
	1690	* cartesian exponents for centers 1 and 2. The shell1 and shell2
	1691	* globals are used. */
	1692	void
	1693	Int1eV3::comp_shell_block_nuclear(int gc1, int a, int gc2, int b,
	1694	int gcsize2, int gcoff1, int gcoff2,
	1695	double coef, double *buffer)
	1696	{
	1697	int i,j,k,xyz;
	1698	double Pi;
	1699	double oozeta;
	1700	double AmB,AmB2;
	1701	double PmC2;
	1702	double auxcoef;
	1703	double tmp;
	1704	int am = a + b;
	1705	int size1 = INT_NCART_NN(a);
	1706	int size2 = INT_NCART_NN(b);
	1707
	1708	#if DEBUG_NUC_SHELL
	1709	double shellints = new double[size1size2];
	1710	memset(shellints,0,sizeof(double)size1size2);
	1711	#endif
	1712
	1713	/* Loop over the primitives in the shells. */
	1714	for (i=0; i<gshell1->nprimitive(); i++) {
	1715	for (j=0; j<gshell2->nprimitive(); j++) {
	1716
	1717	/* Compute the intermediates. */
	1718	zeta = gshell1->exponent(i) + gshell2->exponent(j);
	1719	oozeta = 1.0/zeta;
	1720	oo2zeta = 0.5*oozeta;
	1721	AmB2 = 0.0;
	1722	PmC2 = 0.0;
	1723	for (xyz=0; xyz<3; xyz++) {
	1724	Pi = oozeta(gshell1->exponent(i) A[xyz]
	1725	+ gshell2->exponent(j) * B[xyz]);
	1726	PmA[xyz] = Pi - A[xyz];
	1727	PmB[xyz] = Pi - B[xyz];
	1728	PmC[xyz] = Pi - C[xyz];
	1729	AmB = A[xyz] - B[xyz];
	1730	AmB2 += AmB*AmB;
	1731	PmC2 += PmC[xyz]*PmC[xyz];
	1732	}
	1733
	1734	/* The auxillary integral coeficients. */
	1735	auxcoef = 2.0 * M_PI/(gshell1->exponent(i)
	1736	+gshell2->exponent(j))
	1737	* exp(- oozeta * gshell1->exponent(i)
	1738	* gshell2->exponent(j) * AmB2);
	1739
	1740	/* The Fm(U) intermediates. */
	1741	fjttable_ = fjt_->values(am,zeta*PmC2);
	1742
	1743	/* Convert the Fm(U) intermediates into the auxillary
	1744	* nuclear attraction integrals. */
	1745	for (k=0; k<=am; k++) {
	1746	fjttable_[k] *= auxcoef;
	1747	}
	1748
	1749	/* Compute the primitive nuclear attraction integral. */
	1750	comp_prim_block_nuclear(a,b);
	1751
	1752	tmp = gshell1->coefficient_unnorm(gc1,i)
	1753	* gshell2->coefficient_unnorm(gc2,j)
	1754	* coef;
	1755
	1756	if (exponent_weighted == 0) tmp *= gshell1->exponent(i);
	1757	else if (exponent_weighted == 1) tmp *= gshell2->exponent(j);
	1758
	1759	#if DEBUG_NUC_SHELL
	1760	double *tprimbuffer = inter(a,b,0);
	1761	double *tmpshellints = shellints;
	1762	for (int ip=0; ip<size1; ip++) {
	1763	for (int jp=0; jp<size2; jp++) {
	1764	tmpshellints++ += tmp *tprimbuffer++ / coef;
	1765	}
	1766	}
	1767	#endif
	1768	double *primbuffer = inter(a,b,0);
	1769	for (int ip=0; ip<size1; ip++) {
	1770	for (int jp=0; jp<size2; jp++) {
	1771	//ExEnv::outn() << scprintf("buffer[%d] += %18.15f",
	1772	// (ip+gcoff1)*gcsize2+jp+gcoff2,
	1773	// tmp * *primbuffer)
	1774	// << endl;
	1775	buffer[(ip+gcoff1)gcsize2+jp+gcoff2] += tmp *primbuffer++;
	1776	}
	1777	}
	1778	}
	1779	}
	1780
	1781	#if DEBUG_NUC_SHELL
	1782	# if DEBUG_NUC_SHELL > 1
	1783	ExEnv::outn() << scprintf("GC = (%d %d), A = %d, B = %d", gc1, gc2, a, b)
	1784	<< endl;
	1785	# endif
	1786	int i1,j1,k1;
	1787	int i2,j2,k2;
	1788	int ip = 0;
	1789	double *tmpshellints = shellints;
	1790	FOR_CART(i1,j1,k1,a) {
	1791	int jp = 0;
	1792	FOR_CART(i2,j2,k2,b) {
	1793	double fast = *tmpshellints++;
	1794	double slow = comp_shell_nuclear(gc1, i1, j1, k1,
	1795	gc2, i2, j2, k2);
	1796	int bad = fabs(fast-slow)>1.0e-12;
	1797	if (DEBUG_NUC_SHELL > 1 \|\| bad) {
	1798	ExEnv::outn() << scprintf("NUC SHELL: (%d %d %d) (%d %d %d): ",
	1799	i1,j1,k1, i2,j2,k2)
	1800	<< scprintf(" % 20.15f % 20.15f",
	1801	fast, slow);
	1802	}
	1803	if (bad) {
	1804	ExEnv::outn() << " ****" << endl;
	1805	}
	1806	else if (DEBUG_NUC_SHELL > 1) {
	1807	ExEnv::outn() << endl;
	1808	}
	1809	jp++;
	1810	} END_FOR_CART;
	1811	ip++;
	1812	} END_FOR_CART;
	1813	delete[] shellints;
	1814	#endif
	1815	}
	1816
	1817	void
	1818	Int1eV3::comp_prim_block_nuclear(int a, int b)
	1819	{
	1820	int im, ia, ib;
	1821	int l = a + b;
	1822
	1823	// fill in the ia+ib=0 integrals
	1824	for (im=0; im<=l; im++) {
	1825	#if DEBUG_NUC_PRIM > 1
	1826	ExEnv::outn() << "BUILD: M = " << im
	1827	<< " A = " << 0
	1828	<< " B = " << 0
	1829	<< endl;
	1830	#endif
	1831	inter(0,0,im)[0] = fjttable_[im];
	1832	}
	1833
	1834	for (im=l-1; im>=0; im--) {
	1835	int lm = l-im;
	1836	// build the integrals for a = 0
	1837	for (ib=1; ib<=lm && ib<=b; ib++) {
	1838	#if DEBUG_NUC_PRIM > 1
	1839	ExEnv::outn() << "BUILD: M = " << im
	1840	<< " A = " << 0
	1841	<< " B = " << ib
	1842	<< endl;
	1843	#endif
	1844	comp_prim_block_nuclear_build_b(ib,im);
	1845	}
	1846	for (ia=1; ia<=lm && ia<=a; ia++) {
	1847	for (ib=0; ib<=lm-ia && ib<=b; ib++) {
	1848	#if DEBUG_NUC_PRIM > 1
	1849	ExEnv::outn() << "BUILD: M = " << im
	1850	<< " A = " << ia
	1851	<< " B = " << ib
	1852	<< endl;
	1853	#endif
	1854	comp_prim_block_nuclear_build_a(ia,ib,im);
	1855	}
	1856	}
	1857	}
	1858
	1859	#if DEBUG_NUC_PRIM
	1860	for (im=0; im<=l; im++) {
	1861	int lm = l-im;
	1862	for (ia=0; ia<=lm && ia<=a; ia++) {
	1863	int na = INT_NCART_NN(a);
	1864	for (ib=0; ib<=lm-ia && ib<=b; ib++) {
	1865	int nb = INT_NCART_NN(b);
	1866	#if DEBUG_NUC_PRIM > 1
	1867	ExEnv::outn() << "M = " << im
	1868	<< " A = " << ia
	1869	<< " B = " << ib
	1870	<< endl;
	1871	#endif
	1872	double *buf = inter(ia,ib,im);
	1873	int i1,j1,k1, i2,j2,k2;
	1874	FOR_CART(i1,j1,k1,ia) {
	1875	FOR_CART(i2,j2,k2,ib) {
	1876	double fast = *buf++;
	1877	double slow = comp_prim_nuclear(i1, j1, k1,
	1878	i2, j2, k2, im);
	1879	if (fast > 999.0) fast = 999.0;
	1880	if (fast < -999.0) fast = -999.0;
	1881	if (fabs(fast-slow)>1.0e-12) {
	1882	ExEnv::outn() << scprintf("(%d %d %d) (%d %d %d) (%d): ",
	1883	i1,j1,k1, i2,j2,k2, im)
	1884	<< scprintf(" % 20.15f % 20.15f",
	1885	fast, slow)
	1886	<< endl;
	1887	}
	1888	} END_FOR_CART;
	1889	} END_FOR_CART;
	1890	}
	1891	}
	1892	}
	1893	#endif
	1894	}
	1895
	1896	void
	1897	Int1eV3::comp_prim_block_nuclear_build_a(int a, int b, int m)
	1898	{
	1899	double *I000 = inter(a,b,m);
	1900	double *I100 = inter(a-1,b,m);
	1901	double *I101 = inter(a-1,b,m+1);
	1902	double *I200 = (a>1?inter(a-2,b,m):0);
	1903	double *I201 = (a>1?inter(a-2,b,m+1):0);
	1904	double *I110 = (b?inter(a-1,b-1,m):0);
	1905	double *I111 = (b?inter(a-1,b-1,m+1):0);
	1906	int i1,j1,k1;
	1907	int i2,j2,k2;
	1908	int carta=0;
	1909	int sizeb = INT_NCART_NN(b);
	1910	int sizebm1 = INT_NCART_DEC(b,sizeb);
	1911	FOR_CART(i1,j1,k1,a) {
	1912	int cartb=0;
	1913	FOR_CART(i2,j2,k2,b) {
	1914	double result = 0.0;
	1915	if (i1) {
	1916	int am1 = INT_CARTINDEX(a-1,i1-1,j1);
	1917	result = PmA[0] * I100[am1*sizeb+cartb];
	1918	result -= PmC[0] * I101[am1*sizeb+cartb];
	1919	if (i1>1) {
	1920	int am2 = INT_CARTINDEX(a-2,i1-2,j1);
	1921	result += oo2zeta * (i1-1)
	1922	(I200[am2sizeb+cartb] - I201[am2*sizeb+cartb]);
	1923	}
	1924	if (i2) {
	1925	int bm1 = INT_CARTINDEX(b-1,i2-1,j2);
	1926	result += oo2zeta * i2
	1927	(I110[am1sizebm1+bm1] - I111[am1*sizebm1+bm1]);
	1928	}
	1929	}
	1930	else if (j1) {
	1931	int am1 = INT_CARTINDEX(a-1,i1,j1-1);
	1932	result = PmA[1] * I100[am1*sizeb+cartb];
	1933	result -= PmC[1] * I101[am1*sizeb+cartb];
	1934	if (j1>1) {
	1935	int am2 = INT_CARTINDEX(a-2,i1,j1-2);
	1936	result += oo2zeta * (j1-1)
	1937	(I200[am2sizeb+cartb] - I201[am2*sizeb+cartb]);
	1938	}
	1939	if (j2) {
	1940	int bm1 = INT_CARTINDEX(b-1,i2,j2-1);
	1941	result += oo2zeta * j2
	1942	(I110[am1sizebm1+bm1] - I111[am1*sizebm1+bm1]);
	1943	}
	1944	}
	1945	else if (k1) {
	1946	int am1 = INT_CARTINDEX(a-1,i1,j1);
	1947	result = PmA[2] * I100[am1*sizeb+cartb];
	1948	result -= PmC[2] * I101[am1*sizeb+cartb];
	1949	if (k1>1) {
	1950	int am2 = INT_CARTINDEX(a-2,i1,j1);
	1951	result += oo2zeta * (k1-1)
	1952	(I200[am2sizeb+cartb] - I201[am2*sizeb+cartb]);
	1953	}
	1954	if (k2) {
	1955	int bm1 = INT_CARTINDEX(b-1,i2,j2);
	1956	result += oo2zeta * k2
	1957	(I110[am1sizebm1+bm1] - I111[am1*sizebm1+bm1]);
	1958	}
	1959	}
	1960	I000[carta*sizeb+cartb] = result;
	1961	cartb++;
	1962	} END_FOR_CART;
	1963	carta++;
	1964	} END_FOR_CART;
	1965	}
	1966
	1967	void
	1968	Int1eV3::comp_prim_block_nuclear_build_b(int b, int m)
	1969	{
	1970	double *I000 = inter(0,b,m);
	1971	double *I010 = inter(0,b-1,m);
	1972	double *I011 = inter(0,b-1,m+1);
	1973	double *I020 = (b>1?inter(0,b-2,m):0);
	1974	double *I021 = (b>1?inter(0,b-2,m+1):0);
	1975	int i2,j2,k2;
	1976	int cartb=0;
	1977	FOR_CART(i2,j2,k2,b) {
	1978	double result = 0.0;
	1979
	1980	if (i2) {
	1981	int bm1 = INT_CARTINDEX(b-1,i2-1,j2);
	1982	result = PmB[0] * I010[bm1];
	1983	result -= PmC[0] * I011[bm1];
	1984	if (i2>1) {
	1985	int bm2 = INT_CARTINDEX(b-2,i2-2,j2);
	1986	result += oo2zeta * (i2-1) * (I020[bm2] - I021[bm2]);
	1987	}
	1988	}
	1989	else if (j2) {
	1990	int bm1 = INT_CARTINDEX(b-1,i2,j2-1);
	1991	result = PmB[1] * I010[bm1];
	1992	result -= PmC[1] * I011[bm1];
	1993	if (j2>1) {
	1994	int bm2 = INT_CARTINDEX(b-2,i2,j2-2);
	1995	result += oo2zeta * (j2-1) * (I020[bm2] - I021[bm2]);
	1996	}
	1997	}
	1998	else if (k2) {
	1999	int bm1 = INT_CARTINDEX(b-1,i2,j2);
	2000	result = PmB[2] * I010[bm1];
	2001	result -= PmC[2] * I011[bm1];
	2002	if (k2>1) {
	2003	int bm2 = INT_CARTINDEX(b-2,i2,j2);
	2004	result += oo2zeta * (k2-1) * (I020[bm2] - I021[bm2]);
	2005	}
	2006	}
	2007
	2008	I000[cartb] = result;
	2009	cartb++;
	2010	} END_FOR_CART;
	2011	}
	2012
	2013	double
	2014	Int1eV3::comp_prim_nuclear(int i1, int j1, int k1,
	2015	int i2, int j2, int k2, int m)
	2016	{
	2017	double result;
	2018
	2019	if (i1) {
	2020	result = PmA[0] * comp_prim_nuclear(i1-1,j1,k1,i2,j2,k2,m);
	2021	result -= PmC[0] * comp_prim_nuclear(i1-1,j1,k1,i2,j2,k2,m+1);
	2022	if (i1>1) result += oo2zeta * (i1-1)
	2023	* ( comp_prim_nuclear(i1-2,j1,k1,i2,j2,k2,m)
	2024	- comp_prim_nuclear(i1-2,j1,k1,i2,j2,k2,m+1));
	2025	if (i2) result += oo2zeta * i2
	2026	* ( comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m)
	2027	- comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m+1));
	2028	}
	2029	else if (j1) {
	2030	result = PmA[1] * comp_prim_nuclear(i1,j1-1,k1,i2,j2,k2,m);
	2031	result -= PmC[1] * comp_prim_nuclear(i1,j1-1,k1,i2,j2,k2,m+1);
	2032	if (j1>1) result += oo2zeta * (j1-1)
	2033	* ( comp_prim_nuclear(i1,j1-2,k1,i2,j2,k2,m)
	2034	- comp_prim_nuclear(i1,j1-2,k1,i2,j2,k2,m+1));
	2035	if (j2) result += oo2zeta * j2
	2036	* ( comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m)
	2037	- comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m+1));
	2038	}
	2039	else if (k1) {
	2040	result = PmA[2] * comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2,m);
	2041	result -= PmC[2] * comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2,m+1);
	2042	if (k1>1) result += oo2zeta * (k1-1)
	2043	* ( comp_prim_nuclear(i1,j1,k1-2,i2,j2,k2,m)
	2044	- comp_prim_nuclear(i1,j1,k1-2,i2,j2,k2,m+1));
	2045	if (k2) result += oo2zeta * k2
	2046	* ( comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m)
	2047	- comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m+1));
	2048	}
	2049	else if (i2) {
	2050	result = PmB[0] * comp_prim_nuclear(i1,j1,k1,i2-1,j2,k2,m);
	2051	result -= PmC[0] * comp_prim_nuclear(i1,j1,k1,i2-1,j2,k2,m+1);
	2052	if (i2>1) result += oo2zeta * (i2-1)
	2053	* ( comp_prim_nuclear(i1,j1,k1,i2-2,j2,k2,m)
	2054	- comp_prim_nuclear(i1,j1,k1,i2-2,j2,k2,m+1));
	2055	if (i1) result += oo2zeta * i1
	2056	* ( comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m)
	2057	- comp_prim_nuclear(i1-1,j1,k1,i2-1,j2,k2,m+1));
	2058	}
	2059	else if (j2) {
	2060	result = PmB[1] * comp_prim_nuclear(i1,j1,k1,i2,j2-1,k2,m);
	2061	result -= PmC[1] * comp_prim_nuclear(i1,j1,k1,i2,j2-1,k2,m+1);
	2062	if (j2>1) result += oo2zeta * (j2-1)
	2063	* ( comp_prim_nuclear(i1,j1,k1,i2,j2-2,k2,m)
	2064	- comp_prim_nuclear(i1,j1,k1,i2,j2-2,k2,m+1));
	2065	if (j1) result += oo2zeta * j1
	2066	* ( comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m)
	2067	- comp_prim_nuclear(i1,j1-1,k1,i2,j2-1,k2,m+1));
	2068	}
	2069	else if (k2) {
	2070	result = PmB[2] * comp_prim_nuclear(i1,j1,k1,i2,j2,k2-1,m);
	2071	result -= PmC[2] * comp_prim_nuclear(i1,j1,k1,i2,j2,k2-1,m+1);
	2072	if (k2>1) result += oo2zeta * (k2-1)
	2073	* ( comp_prim_nuclear(i1,j1,k1,i2,j2,k2-2,m)
	2074	- comp_prim_nuclear(i1,j1,k1,i2,j2,k2-2,m+1));
	2075	if (k1) result += oo2zeta * k1
	2076	* ( comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m)
	2077	- comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2-1,m+1));
	2078	}
	2079	else result = fjttable_[m];
	2080
	2081	return result;
	2082	}
	2083
	2084	/* Compute the electric field integral for the shell. The arguments are the
	2085	* the electric field vector, the
	2086	* cartesian exponents for centers 1 and 2. The shell1 and shell2
	2087	* globals are used. */
	2088	void
	2089	Int1eV3::comp_shell_efield(double *efield,
	2090	int gc1, int i1, int j1, int k1,
	2091	int gc2, int i2, int j2, int k2)
	2092	{
	2093	int i,j,k,xyz;
	2094	double result[3];
	2095	double Pi;
	2096	double oozeta;
	2097	double AmB,AmB2;
	2098	double PmC2;
	2099	double auxcoef;
	2100	int am;
	2101
	2102	am = i1+j1+k1+i2+j2+k2;
	2103
	2104	/* Loop over the primitives in the shells. */
	2105	for (xyz=0; xyz<3; xyz++) result[xyz] = 0.0;
	2106	for (i=0; i<gshell1->nprimitive(); i++) {
	2107	for (j=0; j<gshell2->nprimitive(); j++) {
	2108
	2109	/* Compute the intermediates. */
	2110	zeta = gshell1->exponent(i) + gshell2->exponent(j);
	2111	oozeta = 1.0/zeta;
	2112	oo2zeta = 0.5*oozeta;
	2113	AmB2 = 0.0;
	2114	PmC2 = 0.0;
	2115	for (xyz=0; xyz<3; xyz++) {
	2116	Pi = oozeta(gshell1->exponent(i) A[xyz] + gshell2->exponent(j) * B[xyz]);
	2117	PmA[xyz] = Pi - A[xyz];
	2118	PmB[xyz] = Pi - B[xyz];
	2119	PmC[xyz] = Pi - C[xyz];
	2120	AmB = A[xyz] - B[xyz];
	2121	AmB2 += AmB*AmB;
	2122	PmC2 += PmC[xyz]*PmC[xyz];
	2123	}
	2124
	2125	/* The auxillary integral coeficients. */
	2126	auxcoef = 2.0 * M_PI/(gshell1->exponent(i)
	2127	+gshell2->exponent(j))
	2128	* exp(- oozeta * gshell1->exponent(i)
	2129	* gshell2->exponent(j) * AmB2);
	2130
	2131	/* The Fm(U) intermediates. */
	2132	fjttable_ = fjt_->values(am+1,zeta*PmC2);
	2133
	2134	/* Convert the Fm(U) intermediates into the auxillary
	2135	* nuclear attraction integrals. */
	2136	for (k=0; k<=am+1; k++) {
	2137	fjttable_[k] *= auxcoef;
	2138	}
	2139
	2140	/* Compute the nuclear attraction integral. */
	2141	for (xyz=0; xyz<3; xyz++) {
	2142	result[xyz] += gshell1->coefficient_unnorm(gc1,i)
	2143	* gshell2->coefficient_unnorm(gc2,j)
	2144	* comp_prim_efield(xyz,i1,j1,k1,i2,j2,k2,0);
	2145	}
	2146	}
	2147	}
	2148
	2149	for (xyz=0; xyz<3; xyz++) efield[xyz] = result[xyz];
	2150
	2151	}
	2152
	2153	/* Compute the electric field integral for the shell. The arguments are the
	2154	* the electric field vector, the
	2155	* cartesian exponents for centers 1 and 2. The shell1 and shell2
	2156	* globals are used. */
	2157	void
	2158	Int1eV3::comp_shell_block_efield(int gc1, int a, int gc2, int b,
	2159	int gcsize2, int gcoff1, int gcoff2,
	2160	double coef, double *buffer)
	2161	{
	2162	int i,j,k,xyz;
	2163	double Pi;
	2164	double oozeta;
	2165	double AmB,AmB2;
	2166	double PmC2;
	2167	double auxcoef;
	2168	int am = a + b;
	2169	int size1 = INT_NCART_NN(a);
	2170	int size2 = INT_NCART_NN(b);
	2171
	2172	/* Loop over the primitives in the shells. */
	2173	for (i=0; i<gshell1->nprimitive(); i++) {
	2174	for (j=0; j<gshell2->nprimitive(); j++) {
	2175
	2176	/* Compute the intermediates. */
	2177	zeta = gshell1->exponent(i) + gshell2->exponent(j);
	2178	oozeta = 1.0/zeta;
	2179	oo2zeta = 0.5*oozeta;
	2180	AmB2 = 0.0;
	2181	PmC2 = 0.0;
	2182	for (xyz=0; xyz<3; xyz++) {
	2183	Pi = oozeta(gshell1->exponent(i) A[xyz]
	2184	+ gshell2->exponent(j) * B[xyz]);
	2185	PmA[xyz] = Pi - A[xyz];
	2186	PmB[xyz] = Pi - B[xyz];
	2187	PmC[xyz] = Pi - C[xyz];
	2188	AmB = A[xyz] - B[xyz];
	2189	AmB2 += AmB*AmB;
	2190	PmC2 += PmC[xyz]*PmC[xyz];
	2191	}
	2192
	2193	/* The auxillary integral coeficients. */
	2194	auxcoef = 2.0 * M_PI/(gshell1->exponent(i)
	2195	+gshell2->exponent(j))
	2196	* exp(- oozeta * gshell1->exponent(i)
	2197	* gshell2->exponent(j) * AmB2);
	2198
	2199	/* The Fm(U) intermediates. */
	2200	fjttable_ = fjt_->values(am+1,zeta*PmC2);
	2201
	2202	/* Convert the Fm(U) intermediates into the auxillary
	2203	* nuclear attraction integrals. */
	2204	for (k=0; k<=am+1; k++) {
	2205	fjttable_[k] *= auxcoef;
	2206	}
	2207
	2208	double tmp = gshell1->coefficient_unnorm(gc1,i)
	2209	* gshell2->coefficient_unnorm(gc2,j)
	2210	* coef;
	2211
	2212	/* Compute the nuclear attraction integral. */
	2213	comp_prim_block_efield(a,b);
	2214
	2215	double *primbuffer = efield_inter(a,b,0);
	2216	for (int ip=0; ip<size1; ip++) {
	2217	for (int jp=0; jp<size2; jp++) {
	2218	for (xyz=0; xyz<3; xyz++) {
	2219	buffer[((ip+gcoff1)gcsize2+jp+gcoff2)3+xyz]
	2220	+= tmp * *primbuffer++;
	2221	}
	2222	}
	2223	}
	2224	}
	2225	}
	2226
	2227	}
	2228
	2229	void
	2230	Int1eV3::comp_prim_block_efield(int a, int b)
	2231	{
	2232	int xyz;
	2233	int im, ia, ib;
	2234	int l = a + b;
	2235
	2236	// Fill in the nuclear integrals which are needed as intermediates.
	2237	// m=0 is not needed.
	2238
	2239	// fill in the ia+ib=0 integrals, skipping m=0
	2240	for (im=1; im<=l; im++) {
	2241	#if DEBUG_EFIELD_PRIM > 1
	2242	ExEnv::outn() << "BUILD NUC: M = " << im
	2243	<< " A = " << 0
	2244	<< " B = " << 0
	2245	<< endl;
	2246	#endif
	2247	inter(0,0,im)[0] = fjttable_[im];
	2248	}
	2249
	2250	// skipping m=0
	2251	for (im=l-1; im>0; im--) {
	2252	int lm = l-im;
	2253	// build the integrals for a = 0
	2254	for (ib=1; ib<=lm && ib<=b; ib++) {
	2255	#if DEBUG_EFIELD_PRIM > 1
	2256	ExEnv::outn() << "BUILD NUC: M = " << im
	2257	<< " A = " << 0
	2258	<< " B = " << ib
	2259	<< endl;
	2260	#endif
	2261	comp_prim_block_nuclear_build_b(ib,im);
	2262	}
	2263	for (ia=1; ia<=lm && ia<=a; ia++) {
	2264	for (ib=0; ib<=lm-ia && ib<=b; ib++) {
	2265	#if DEBUG_EFIELD_PRIM > 1
	2266	ExEnv::outn() << "BUILD NUC: M = " << im
	2267	<< " A = " << ia
	2268	<< " B = " << ib
	2269	<< endl;
	2270	#endif
	2271	comp_prim_block_nuclear_build_a(ia,ib,im);
	2272	}
	2273	}
	2274	}
	2275
	2276	// now complete the efield integrals
	2277
	2278	// fill in the ia+ib=0 integrals
	2279	for (im=0; im<=l; im++) {
	2280	#if DEBUG_EFIELD_PRIM > 1
	2281	ExEnv::outn() << "BUILD EFIELD: M = " << im
	2282	<< " A = " << 0
	2283	<< " B = " << 0
	2284	<< endl;
	2285	#endif
	2286	double *tmp = efield_inter(0,0,im);
	2287	for (xyz=0; xyz<3; xyz++) {
	2288	tmp++ = 2.0 zeta * PmC[xyz] * fjttable_[im+1];
	2289	}
	2290	}
	2291
	2292	for (im=l-1; im>=0; im--) {
	2293	int lm = l-im;
	2294	// build the integrals for a = 0
	2295	for (ib=1; ib<=lm && ib<=b; ib++) {
	2296	#if DEBUG_EFIELD_PRIM > 1
	2297	ExEnv::outn() << "BUILD EFIELD: M = " << im
	2298	<< " A = " << 0
	2299	<< " B = " << ib
	2300	<< endl;
	2301	#endif
	2302	comp_prim_block_efield_build_b(ib,im);
	2303	}
	2304	for (ia=1; ia<=lm && ia<=a; ia++) {
	2305	for (ib=0; ib<=lm-ia && ib<=b; ib++) {
	2306	#if DEBUG_EFIELD_PRIM > 1
	2307	ExEnv::outn() << "BUILD EFIELD: M = " << im
	2308	<< " A = " << ia
	2309	<< " B = " << ib
	2310	<< endl;
	2311	#endif
	2312	comp_prim_block_efield_build_a(ia,ib,im);
	2313	}
	2314	}
	2315	}
	2316
	2317	#if DEBUG_EFIELD_PRIM
	2318	for (im=0; im<=l; im++) {
	2319	int lm = l-im;
	2320	for (ia=0; ia<=lm && ia<=a; ia++) {
	2321	int na = INT_NCART_NN(a);
	2322	for (ib=0; ib<=lm-ia && ib<=b; ib++) {
	2323	int nb = INT_NCART_NN(b);
	2324	#if DEBUG_EFIELD_PRIM > 1
	2325	ExEnv::outn() << "M = " << im
	2326	<< " A = " << ia
	2327	<< " B = " << ib
	2328	<< endl;
	2329	#endif
	2330	double *buf = efield_inter(ia,ib,im);
	2331	int i1,j1,k1, i2,j2,k2;
	2332	FOR_CART(i1,j1,k1,ia) {
	2333	FOR_CART(i2,j2,k2,ib) {
	2334	for (xyz=0; xyz<3; xyz++) {
	2335	double fast = *buf++;
	2336	double slow = comp_prim_efield(xyz, i1, j1, k1,
	2337	i2, j2, k2, im);
	2338	if (fast > 999.0) fast = 999.0;
	2339	if (fast < -999.0) fast = -999.0;
	2340	if (fabs(fast-slow)>1.0e-12) {
	2341	ExEnv::outn() << scprintf("(%d %d %d) (%d %d %d) (%d): ",
	2342	i1,j1,k1, i2,j2,k2, im)
	2343	<< scprintf(" % 20.15f % 20.15f",
	2344	fast, slow)
	2345	<< endl;
	2346	}
	2347	}
	2348	} END_FOR_CART;
	2349	} END_FOR_CART;
	2350	}
	2351	}
	2352	}
	2353	#endif
	2354	}
	2355
	2356	void
	2357	Int1eV3::comp_prim_block_efield_build_a(int a, int b, int m)
	2358	{
	2359	double *I000 = efield_inter(a,b,m);
	2360	double *I100 = efield_inter(a-1,b,m);
	2361	double *I101 = efield_inter(a-1,b,m+1);
	2362	double *I200 = (a>1?efield_inter(a-2,b,m):0);
	2363	double *I201 = (a>1?efield_inter(a-2,b,m+1):0);
	2364	double *I110 = (b?efield_inter(a-1,b-1,m):0);
	2365	double *I111 = (b?efield_inter(a-1,b-1,m+1):0);
	2366	double *nucI101 = inter(a-1,b,m+1);
	2367	int i1,j1,k1;
	2368	int i2,j2,k2;
	2369	int xyz;
	2370	int carta=0;
	2371	int sizeb = INT_NCART_NN(b);
	2372	int sizebm1 = INT_NCART_DEC(b,sizeb);
	2373	FOR_CART(i1,j1,k1,a) {
	2374	int cartb=0;
	2375	FOR_CART(i2,j2,k2,b) {
	2376	for (xyz=0; xyz<3; xyz++) {
	2377	double result = 0.0;
	2378	if (i1) {
	2379	int am1 = INT_CARTINDEX(a-1,i1-1,j1);
	2380	result = PmA[0] * I100[(am1sizeb+cartb)3+xyz];
	2381	result -= PmC[0] * I101[(am1sizeb+cartb)3+xyz];
	2382	if (i1>1) {
	2383	int am2 = INT_CARTINDEX(a-2,i1-2,j1);
	2384	result += oo2zeta * (i1-1)
	2385	(I200[(am2sizeb+cartb)*3+xyz]
	2386	- I201[(am2sizeb+cartb)3+xyz]);
	2387	}
	2388	if (i2) {
	2389	int bm1 = INT_CARTINDEX(b-1,i2-1,j2);
	2390	result += oo2zeta * i2
	2391	(I110[(am1sizebm1+bm1)*3+xyz]
	2392	- I111[(am1sizebm1+bm1)3+xyz]);
	2393	}
	2394	if (xyz==0) result += nucI101[am1*sizeb+cartb];
	2395	}
	2396	else if (j1) {
	2397	int am1 = INT_CARTINDEX(a-1,i1,j1-1);
	2398	result = PmA[1] * I100[(am1sizeb+cartb)3+xyz];
	2399	result -= PmC[1] * I101[(am1sizeb+cartb)3+xyz];
	2400	if (j1>1) {
	2401	int am2 = INT_CARTINDEX(a-2,i1,j1-2);
	2402	result += oo2zeta * (j1-1)
	2403	(I200[(am2sizeb+cartb)*3+xyz]
	2404	- I201[(am2sizeb+cartb)3+xyz]);
	2405	}
	2406	if (j2) {
	2407	int bm1 = INT_CARTINDEX(b-1,i2,j2-1);
	2408	result += oo2zeta * j2
	2409	(I110[(am1sizebm1+bm1)*3+xyz]
	2410	- I111[(am1sizebm1+bm1)3+xyz]);
	2411	}
	2412	if (xyz==1) result += nucI101[am1*sizeb+cartb];
	2413	}
	2414	else if (k1) {
	2415	int am1 = INT_CARTINDEX(a-1,i1,j1);
	2416	result = PmA[2] * I100[(am1sizeb+cartb)3+xyz];
	2417	result -= PmC[2] * I101[(am1sizeb+cartb)3+xyz];
	2418	if (k1>1) {
	2419	int am2 = INT_CARTINDEX(a-2,i1,j1);
	2420	result += oo2zeta * (k1-1)
	2421	(I200[(am2sizeb+cartb)*3+xyz]
	2422	- I201[(am2sizeb+cartb)3+xyz]);
	2423	}
	2424	if (k2) {
	2425	int bm1 = INT_CARTINDEX(b-1,i2,j2);
	2426	result += oo2zeta * k2
	2427	(I110[(am1sizebm1+bm1)*3+xyz]
	2428	- I111[(am1sizebm1+bm1)3+xyz]);
	2429	}
	2430	if (xyz==2) result += nucI101[am1*sizeb+cartb];
	2431	}
	2432	I000[(cartasizeb+cartb)3+xyz] = result;
	2433	}
	2434	cartb++;
	2435	} END_FOR_CART;
	2436	carta++;
	2437	} END_FOR_CART;
	2438	}
	2439
	2440	void
	2441	Int1eV3::comp_prim_block_efield_build_b(int b, int m)
	2442	{
	2443	double *I000 = efield_inter(0,b,m);
	2444	double *I010 = efield_inter(0,b-1,m);
	2445	double *I011 = efield_inter(0,b-1,m+1);
	2446	double *I020 = (b>1?efield_inter(0,b-2,m):0);
	2447	double *I021 = (b>1?efield_inter(0,b-2,m+1):0);
	2448	double *nucI011 = inter(0,b-1,m+1);
	2449	int xyz;
	2450	int i2,j2,k2;
	2451	int cartb=0;
	2452	FOR_CART(i2,j2,k2,b) {
	2453	for (xyz=0; xyz<3; xyz++) {
	2454	double result = 0.0;
	2455
	2456	if (i2) {
	2457	int bm1 = INT_CARTINDEX(b-1,i2-1,j2);
	2458	result = PmB[0] * I010[(bm1)*3+xyz];
	2459	result -= PmC[0] * I011[(bm1)*3+xyz];
	2460	if (i2>1) {
	2461	int bm2 = INT_CARTINDEX(b-2,i2-2,j2);
	2462	result += oo2zeta * (i2-1) * (I020[(bm2)*3+xyz]
	2463	- I021[(bm2)*3+xyz]);
	2464	}
	2465	if (xyz==0) result += nucI011[bm1];
	2466	}
	2467	else if (j2) {
	2468	int bm1 = INT_CARTINDEX(b-1,i2,j2-1);
	2469	result = PmB[1] * I010[(bm1)*3+xyz];
	2470	result -= PmC[1] * I011[(bm1)*3+xyz];
	2471	if (j2>1) {
	2472	int bm2 = INT_CARTINDEX(b-2,i2,j2-2);
	2473	result += oo2zeta * (j2-1) * (I020[(bm2)*3+xyz]
	2474	- I021[(bm2)*3+xyz]);
	2475	}
	2476	if (xyz==1) result += nucI011[bm1];
	2477	}
	2478	else if (k2) {
	2479	int bm1 = INT_CARTINDEX(b-1,i2,j2);
	2480	result = PmB[2] * I010[(bm1)*3+xyz];
	2481	result -= PmC[2] * I011[(bm1)*3+xyz];
	2482	if (k2>1) {
	2483	int bm2 = INT_CARTINDEX(b-2,i2,j2);
	2484	result += oo2zeta * (k2-1) * (I020[(bm2)*3+xyz]
	2485	- I021[(bm2)*3+xyz]);
	2486	}
	2487	if (xyz==2) result += nucI011[bm1];
	2488	}
	2489
	2490	I000[(cartb)*3+xyz] = result;
	2491	}
	2492	cartb++;
	2493	} END_FOR_CART;
	2494	}
	2495
	2496	double
	2497	Int1eV3::comp_prim_efield(int xyz, int i1, int j1, int k1,
	2498	int i2, int j2, int k2, int m)
	2499	{
	2500	double result;
	2501
	2502	/* if ((xyz != 0) \|\| (i1 != 1)) return 0.0; */
	2503
	2504	if (i1) {
	2505	result = PmA[0] * comp_prim_efield(xyz,i1-1,j1,k1,i2,j2,k2,m);
	2506	result -= PmC[0] * comp_prim_efield(xyz,i1-1,j1,k1,i2,j2,k2,m+1);
	2507	if (i1>1) result += oo2zeta * (i1-1)
	2508	* ( comp_prim_efield(xyz,i1-2,j1,k1,i2,j2,k2,m)
	2509	- comp_prim_efield(xyz,i1-2,j1,k1,i2,j2,k2,m+1));
	2510	if (i2) result += oo2zeta * i2
	2511	* ( comp_prim_efield(xyz,i1-1,j1,k1,i2-1,j2,k2,m)
	2512	- comp_prim_efield(xyz,i1-1,j1,k1,i2-1,j2,k2,m+1));
	2513	if (xyz==0) result += comp_prim_nuclear(i1-1,j1,k1,i2,j2,k2,m+1);
	2514	}
	2515	else if (j1) {
	2516	result = PmA[1] * comp_prim_efield(xyz,i1,j1-1,k1,i2,j2,k2,m);
	2517	result -= PmC[1] * comp_prim_efield(xyz,i1,j1-1,k1,i2,j2,k2,m+1);
	2518	if (j1>1) result += oo2zeta * (j1-1)
	2519	* ( comp_prim_efield(xyz,i1,j1-2,k1,i2,j2,k2,m)
	2520	- comp_prim_efield(xyz,i1,j1-2,k1,i2,j2,k2,m+1));
	2521	if (j2) result += oo2zeta * j2
	2522	* ( comp_prim_efield(xyz,i1,j1-1,k1,i2,j2-1,k2,m)
	2523	- comp_prim_efield(xyz,i1,j1-1,k1,i2,j2-1,k2,m+1));
	2524	if (xyz==1) result += comp_prim_nuclear(i1,j1-1,k1,i2,j2,k2,m+1);
	2525	}
	2526	else if (k1) {
	2527	result = PmA[2] * comp_prim_efield(xyz,i1,j1,k1-1,i2,j2,k2,m);
	2528	result -= PmC[2] * comp_prim_efield(xyz,i1,j1,k1-1,i2,j2,k2,m+1);
	2529	if (k1>1) result += oo2zeta * (k1-1)
	2530	* ( comp_prim_efield(xyz,i1,j1,k1-2,i2,j2,k2,m)
	2531	- comp_prim_efield(xyz,i1,j1,k1-2,i2,j2,k2,m+1));
	2532	if (k2) result += oo2zeta * k2
	2533	* ( comp_prim_efield(xyz,i1,j1,k1-1,i2,j2,k2-1,m)
	2534	- comp_prim_efield(xyz,i1,j1,k1-1,i2,j2,k2-1,m+1));
	2535	if (xyz==2) result += comp_prim_nuclear(i1,j1,k1-1,i2,j2,k2,m+1);
	2536	}
	2537	else if (i2) {
	2538	result = PmB[0] * comp_prim_efield(xyz,i1,j1,k1,i2-1,j2,k2,m);
	2539	result -= PmC[0] * comp_prim_efield(xyz,i1,j1,k1,i2-1,j2,k2,m+1);
	2540	if (i2>1) result += oo2zeta * (i2-1)
	2541	* ( comp_prim_efield(xyz,i1,j1,k1,i2-2,j2,k2,m)
	2542	- comp_prim_efield(xyz,i1,j1,k1,i2-2,j2,k2,m+1));
	2543	if (i1) result += oo2zeta * i1
	2544	* ( comp_prim_efield(xyz,i1-1,j1,k1,i2-1,j2,k2,m)
	2545	- comp_prim_efield(xyz,i1-1,j1,k1,i2-1,j2,k2,m+1));
	2546	if (xyz==0) result += comp_prim_nuclear(i1,j1,k1,i2-1,j2,k2,m+1);
	2547	}
	2548	else if (j2) {
	2549	result = PmB[1] * comp_prim_efield(xyz,i1,j1,k1,i2,j2-1,k2,m);
	2550	result -= PmC[1] * comp_prim_efield(xyz,i1,j1,k1,i2,j2-1,k2,m+1);
	2551	if (j2>1) result += oo2zeta * (j2-1)
	2552	* ( comp_prim_efield(xyz,i1,j1,k1,i2,j2-2,k2,m)
	2553	- comp_prim_efield(xyz,i1,j1,k1,i2,j2-2,k2,m+1));
	2554	if (j1) result += oo2zeta * j1
	2555	* ( comp_prim_efield(xyz,i1,j1-1,k1,i2,j2-1,k2,m)
	2556	- comp_prim_efield(xyz,i1,j1-1,k1,i2,j2-1,k2,m+1));
	2557	if (xyz==1) result += comp_prim_nuclear(i1,j1,k1,i2,j2-1,k2,m+1);
	2558	}
	2559	else if (k2) {
	2560	result = PmB[2] * comp_prim_efield(xyz,i1,j1,k1,i2,j2,k2-1,m);
	2561	result -= PmC[2] * comp_prim_efield(xyz,i1,j1,k1,i2,j2,k2-1,m+1);
	2562	if (k2>1) result += oo2zeta * (k2-1)
	2563	* ( comp_prim_efield(xyz,i1,j1,k1,i2,j2,k2-2,m)
	2564	- comp_prim_efield(xyz,i1,j1,k1,i2,j2,k2-2,m+1));
	2565	if (k1) result += oo2zeta * k1
	2566	* ( comp_prim_efield(xyz,i1,j1,k1-1,i2,j2,k2-1,m)
	2567	- comp_prim_efield(xyz,i1,j1,k1-1,i2,j2,k2-1,m+1));
	2568	if (xyz==2) result += comp_prim_nuclear(i1,j1,k1,i2,j2,k2-1,m+1);
	2569	}
	2570	else {
	2571	/* We arrive here if we have a (s\| \|s) type efield integral.
	2572	* The fjttable contains the standard (s\| \|s) nuc attr integrals.
	2573	*/
	2574	result = 2.0 * zeta * PmC[xyz] * fjttable_[m+1];
	2575	}
	2576
	2577	return result;
	2578	}
	2579
	2580
	2581	/* --------------------------------------------------------------- */
	2582	/* ------------- Routines for dipole moment integrals ------------ */
	2583	/* --------------------------------------------------------------- */
	2584
	2585	/* This computes the dipole integrals between functions in two shells.
	2586	* The result is accumulated in the buffer in the form bf1 x y z, bf2
	2587	* x y z, etc. The last arg, com, is the origin of the coordinate
	2588	* system used to compute the dipole moment.
	2589	*/
	2590	void
	2591	Int1eV3::int_accum_shell_dipole(int ish, int jsh,
	2592	double *com)
	2593	{
	2594	int c1,i1,j1,k1,c2,i2,j2,k2;
	2595	int gc1,gc2;
	2596	int index,index1,index2;
	2597	double dipole[3];
	2598	int xyz;
	2599
	2600	for (xyz=0; xyz<3; xyz++) C[xyz] = com[xyz];
	2601
	2602	c1 = bs1_->shell_to_center(ish);
	2603	c2 = bs2_->shell_to_center(jsh);
	2604	for (xyz=0; xyz<3; xyz++) {
	2605	A[xyz] = bs1_->r(c1,xyz);
	2606	B[xyz] = bs2_->r(c2,xyz);
	2607	}
	2608	gshell1 = &bs1_->shell(ish);
	2609	gshell2 = &bs2_->shell(jsh);
	2610	index = 0;
	2611	FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
	2612	FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
	2613	comp_shell_dipole(dipole,gc1,i1,j1,k1,gc2,i2,j2,k2);
	2614	for(mu=0; mu < 3; mu++) {
	2615	cartesianbuffer[index] = dipole[mu];
	2616	index++;
	2617	}
	2618	END_FOR_GCCART_GS(index2)
	2619	END_FOR_GCCART_GS(index1)
	2620	accum_transform_1e_xyz(integral_,
	2621	cartesianbuffer, buff, gshell1, gshell2);
	2622	}
	2623
	2624	/* This computes the dipole integrals between functions in two shells.
	2625	* The result is placed in the buffer in the form bf1 x y z, bf2
	2626	* x y z, etc. The last arg, com, is the origin of the coordinate
	2627	* system used to compute the dipole moment.
	2628	*/
	2629	void
	2630	Int1eV3::dipole(int ish, int jsh, double *com)
	2631	{
	2632	int c1,i1,j1,k1,c2,i2,j2,k2;
	2633	int gc1,gc2;
	2634	int index,index1,index2;
	2635	double dipole[3];
	2636	int xyz;
	2637
	2638	for (xyz=0; xyz<3; xyz++) C[xyz] = com[xyz];
	2639
	2640	c1 = bs1_->shell_to_center(ish);
	2641	c2 = bs2_->shell_to_center(jsh);
	2642	for (xyz=0; xyz<3; xyz++) {
	2643	A[xyz] = bs1_->r(c1,xyz);
	2644	B[xyz] = bs2_->r(c2,xyz);
	2645	}
	2646	gshell1 = &bs1_->shell(ish);
	2647	gshell2 = &bs2_->shell(jsh);
	2648	index = 0;
	2649	FOR_GCCART_GS(gc1,index1,i1,j1,k1,gshell1)
	2650	FOR_GCCART_GS(gc2,index2,i2,j2,k2,gshell2)
	2651	comp_shell_dipole(dipole,gc1,i1,j1,k1,gc2,i2,j2,k2);
	2652	for(mu=0; mu < 3; mu++) {
	2653	cartesianbuffer[index] = dipole[mu];
	2654	index++;
	2655	}
	2656	END_FOR_GCCART_GS(index2)
	2657	END_FOR_GCCART_GS(index1)
	2658	transform_1e_xyz(integral_,
	2659	cartesianbuffer, buff, gshell1, gshell2);
	2660	}
	2661
	2662	void
	2663	Int1eV3::comp_shell_dipole(double* dipole,
	2664	int gc1, int i1, int j1, int k1,
	2665	int gc2, int i2, int j2, int k2)
	2666	{
	2667	double exp1,exp2;
	2668	int i,j,xyz;
	2669	double Pi;
	2670	double oozeta;
	2671	double AmB,AmB2;
	2672	double tmp;
	2673
	2674	dipole[0] = dipole[1] = dipole[2] = 0.0;
	2675
	2676	if ((i1<0)\|\|(j1<0)\|\|(k1<0)\|\|(i2<0)\|\|(j2<0)\|\|(k2<0)) return;
	2677
	2678	/* Loop over the primitives in the shells. */
	2679	for (i=0; i<gshell1->nprimitive(); i++) {
	2680	for (j=0; j<gshell2->nprimitive(); j++) {
	2681
	2682	/* Compute the intermediates. */
	2683	exp1 = gshell1->exponent(i);
	2684	exp2 = gshell2->exponent(j);
	2685	oozeta = 1.0/(exp1 + exp2);
	2686	oo2zeta = 0.5*oozeta;
	2687	AmB2 = 0.0;
	2688	for (xyz=0; xyz<3; xyz++) {
	2689	Pi = oozeta(exp1 A[xyz] + exp2 * B[xyz]);
	2690	PmA[xyz] = Pi - A[xyz];
	2691	PmB[xyz] = Pi - B[xyz];
	2692	PmC[xyz] = Pi - C[xyz];
	2693	AmB = A[xyz] - B[xyz];
	2694	AmB2 += AmB*AmB;
	2695	}
	2696	ss = pow(M_PI/(exp1+exp2),1.5)
	2697	* exp(- oozeta * exp1 * exp2 * AmB2);
	2698	for (mu=0; mu<3; mu++) sMus[mu] = ss * PmC[mu];
	2699	tmp = gshell1->coefficient_unnorm(gc1,i)
	2700	* gshell2->coefficient_unnorm(gc2,j);
	2701	if (exponent_weighted == 0) tmp *= exp1;
	2702	else if (exponent_weighted == 1) tmp *= exp2;
	2703	dipole[0] += tmp * comp_prim_dipole(0,i1,j1,k1,i2,j2,k2);
	2704	dipole[1] += tmp * comp_prim_dipole(1,i1,j1,k1,i2,j2,k2);
	2705	dipole[2] += tmp * comp_prim_dipole(2,i1,j1,k1,i2,j2,k2);
	2706	}
	2707	}
	2708
	2709	}
	2710
	2711	double
	2712	Int1eV3::comp_prim_dipole(int axis,
	2713	int i1, int j1, int k1,
	2714	int i2, int j2, int k2)
	2715	{
	2716	double result;
	2717
	2718	if (i1) {
	2719	result = PmA[0] * comp_prim_dipole(axis,i1-1,j1,k1,i2,j2,k2);
	2720	if (i2)
	2721	result += oo2zetai2comp_prim_dipole(axis,i1-1,j1,k1,i2-1,j2,k2);
	2722	if (i1>1)
	2723	result += oo2zeta(i1-1)comp_prim_dipole(axis,i1-2,j1,k1,i2,j2,k2);
	2724	if(axis==0) result += oo2zeta*comp_prim_overlap(i1-1,j1,k1,i2,j2,k2);
	2725	return result;
	2726	}
	2727	if (j1) {
	2728	result = PmA[1] * comp_prim_dipole(axis,i1,j1-1,k1,i2,j2,k2);
	2729	if (j2)
	2730	result += oo2zetaj2comp_prim_dipole(axis,i1,j1-1,k1,i2,j2-1,k2);
	2731	if (j1>1)
	2732	result += oo2zeta(j1-1)comp_prim_dipole(axis,i1,j1-2,k1,i2,j2,k2);
	2733	if(axis==1) result += oo2zeta*comp_prim_overlap(i1,j1-1,k1,i2,j2,k2);
	2734	return result;
	2735	}
	2736	if (k1) {
	2737	result = PmA[2] * comp_prim_dipole(axis,i1,j1,k1-1,i2,j2,k2);
	2738	if (k2)
	2739	result += oo2zetak2comp_prim_dipole(axis,i1,j1,k1-1,i2,j2,k2-1);
	2740	if (k1>1)
	2741	result += oo2zeta(k1-1)comp_prim_dipole(axis,i1,j1,k1-2,i2,j2,k2);
	2742	if(axis==2) result += oo2zeta*comp_prim_overlap(i1,j1,k1-1,i2,j2,k2);
	2743	return result;
	2744	}
	2745	if (i2) {
	2746	result = PmB[0] * comp_prim_dipole(axis,i1,j1,k1,i2-1,j2,k2);
	2747	if (i1)
	2748	result += oo2zetai1comp_prim_dipole(axis,i1-1,j1,k1,i2-1,j2,k2);
	2749	if (i2>1)
	2750	result += oo2zeta(i2-1)comp_prim_dipole(axis,i1,j1,k1,i2-2,j2,k2);
	2751	if(axis==0) result += oo2zeta*comp_prim_overlap(i1,j1,k1,i2-1,j2,k2);
	2752	return result;
	2753	}
	2754	if (j2) {
	2755	result = PmB[1] * comp_prim_dipole(axis,i1,j1,k1,i2,j2-1,k2);
	2756	if (j1)
	2757	result += oo2zetai1comp_prim_dipole(axis,i1,j1-1,k1,i2,j2-1,k2);
	2758	if (j2>1)
	2759	result += oo2zeta(j2-1)comp_prim_dipole(axis,i1,j1,k1,i2,j2-2,k2);
	2760	if(axis==1) result += oo2zeta*comp_prim_overlap(i1,j1,k1,i2,j2-1,k2);
	2761	return result;
	2762	}
	2763	if (k2) {
	2764	result = PmB[2] * comp_prim_dipole(axis,i1,j1,k1,i2,j2,k2-1);
	2765	if (k1)
	2766	result += oo2zetai1comp_prim_dipole(axis,i1,j1,k1-1,i2,j2,k2-1);
	2767	if (k2>1)
	2768	result += oo2zeta(k2-1)comp_prim_dipole(axis,i1,j1,k1,i2,j2,k2-2);
	2769	if(axis==2) result += oo2zeta*comp_prim_overlap(i1,j1,k1,i2,j2,k2-1);
	2770	return result;
	2771	}
	2772
	2773	return sMus[axis];
	2774	}
	2775
	2776	/////////////////////////////////////////////////////////////////////////////
	2777
	2778	// Local Variables:
	2779	// mode: c++
	2780	// c-file-style: "CLJ-CONDENSED"
	2781	// End:

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

Download in other formats: