3
\brief One-electron integral second derivatives
10
#include <libciomr/libciomr.h>
11
#include <libint/libint.h>
12
#include <libderiv/libderiv.h>
13
#include <Tools/prints.h>
18
#include "oe_deriv2_osrr.h"
19
#include"symmetrize.h"
21
#include"taylor_fm_eval.h"
23
#include "small_fns.h"
25
#define INCLUDE_OVERLAP 1
26
#define INCLUDE_KINETIC 1
27
#define INCLUDE_POTENTIAL 1
32
namespace psi { namespace CINTS {
33
/*--- These frequently used numbers are to avoid costs of passing parameters ---*/
34
static double oo2g, oog, gam;
36
/*-------------------------------
37
Explicit function declarations
38
-------------------------------*/
39
static double overlap_int(double a1, int l1, int m1, int n1, double norm1,
40
double a2, int l2, int m2, int n2, double norm2,
41
struct coordinates AB,
42
struct coordinates PA,
43
struct coordinates PB);
44
static double ke_int(double a1, int l1, int m1, int n1, double norm1,
45
double a2, int l2, int m2, int n2, double norm2,
46
struct coordinates AB,
47
struct coordinates PA,
48
struct coordinates PB);
49
static double f_n(int k, int l1, int l2, double A, double B);
50
static double int_pow(double a, int p);
52
/*!-------------------------------------------------------------
53
This function computes derivatives of one-electron integrals
54
-------------------------------------------------------------*/
57
/* Computing up to second-order derivatives here */
58
const int deriv_lvl = 2;
59
const int deriv0_lvl = 0;
60
const int deriv1_lvl = 1;
61
const int deriv2_lvl = 2;
63
double **hess_ov, **hess_oe;
65
struct coordinates PA, PB, AB, PC;
66
struct shell_pair *sp;
67
int i, j, k, l, ii, jj, kk, ll;
74
int l1, l2, m1, m2, n1, n2;
75
int ioffset, joffset ;
81
int ni,li,nj,lj,ai,aj;
83
int ixm, iym, izm, jxm, jym, jzm;
84
int ixm1, iym1, izm1, jxm1, jym1, jzm1;
85
int indmax, iind,jind;
86
int atom, atom1, atom2;
87
int coord_ax, coord_ay, coord_az,
88
coord_bx, coord_by, coord_bz,
89
coord_cx, coord_cy, coord_cz;
95
double inorm, jnorm, over_pf;
96
double norm_pf, normover_pf, dens_pf, wdens_pf, zdens_pf, znormover_pf, hds_norm_pf;
97
double twozeta_a, twozeta_b, upuppfac, updownpfac, s_int, t_int, v_int;
98
double *ptr1, *ptr2, norm1, norm12;
100
double ***AIX, ***AIY, ***AIZ;
101
double ***AIXX, ***AIXY, ***AIXZ,
102
***AIYY, ***AIYZ, ***AIZZ;
105
/*--- +4*deriv_lvl because of the way we invoke AI_Deriv2_OSrecurs ---*/
106
init_Taylor_Fm_Eval(BasisSet.max_am*4-4+4*deriv_lvl,UserOptions.cutoff);
109
hess_ov = block_matrix(Molecule.num_atoms*3,Molecule.num_atoms*3);
110
hess_oe = block_matrix(Molecule.num_atoms*3,Molecule.num_atoms*3);
112
indmax = (BasisSet.max_am+deriv_lvl-1)*(BasisSet.max_am+deriv_lvl)*(BasisSet.max_am+deriv_lvl)+1;
113
AI0 = init_box(indmax,indmax,2*(BasisSet.max_am+deriv2_lvl)+1);
114
AIX = init_box(indmax,indmax,2*(BasisSet.max_am+deriv1_lvl)+1);
115
AIY = init_box(indmax,indmax,2*(BasisSet.max_am+deriv1_lvl)+1);
116
AIZ = init_box(indmax,indmax,2*(BasisSet.max_am+deriv1_lvl)+1);
117
AIXX = init_box(indmax,indmax,2*(BasisSet.max_am+deriv0_lvl)+1);
118
AIXY = init_box(indmax,indmax,2*(BasisSet.max_am+deriv0_lvl)+1);
119
AIXZ = init_box(indmax,indmax,2*(BasisSet.max_am+deriv0_lvl)+1);
120
AIYY = init_box(indmax,indmax,2*(BasisSet.max_am+deriv0_lvl)+1);
121
AIYZ = init_box(indmax,indmax,2*(BasisSet.max_am+deriv0_lvl)+1);
122
AIZZ = init_box(indmax,indmax,2*(BasisSet.max_am+deriv0_lvl)+1);
124
for (si=0; si<BasisSet.num_shells; si++){
125
ni = ioff[BasisSet.shells[si].am];
126
am_i = BasisSet.shells[si].am-1;
131
iym1 = am_i+deriv2_lvl+1;
133
atom1 = BasisSet.shells[si].center-1;
134
si_fao = BasisSet.shells[si].fao-1;
137
coord_ay = coord_ax+1;
138
coord_az = coord_ax+2;
140
for (sj=0; sj<=si; sj++){
141
nj = ioff[BasisSet.shells[sj].am];
142
am_j = BasisSet.shells[sj].am-1;
147
jym1 = am_j+deriv2_lvl+1;
149
atom2 = BasisSet.shells[sj].center-1;
150
sj_fao = BasisSet.shells[sj].fao - 1;
153
coord_by = coord_bx+1;
154
coord_bz = coord_bx+2;
156
sp = &(BasisSet.shell_pairs[si][sj]);
164
/*--- contract by primitives here ---*/
165
for (i = 0; i < BasisSet.shells[si].n_prims; i++) {
167
twozeta_a = 2.0 * a1;
168
inorm = sp->inorm[i];
169
for (j = 0; j < BasisSet.shells[sj].n_prims; j++) {
171
twozeta_b = 2.0 * a2;
172
gam = sp->gamma[i][j];
173
jnorm = sp->jnorm[j];
174
PA.x = sp->PA[i][j][0];
175
PA.y = sp->PA[i][j][1];
176
PA.z = sp->PA[i][j][2];
177
PB.x = sp->PB[i][j][0];
178
PB.y = sp->PB[i][j][1];
179
PB.z = sp->PB[i][j][2];
181
over_pf = exp(-a1*a2*ab2*oog)*sqrt(M_PI*oog)*M_PI*oog*inorm*jnorm;
183
/*--- create all am components of si ---*/
185
for(ii = 0; ii <= am_i; ii++){
187
for(jj = 0; jj <= ii; jj++){
192
/*--- create all am components of sj ---*/
194
for(kk = 0; kk <= am_j; kk++){
196
for(ll = 0; ll <= kk; ll++){
201
if (si == sj && ai < aj)
204
norm_pf = (GTOs.bf_norm[am_i][ai] * GTOs.bf_norm[am_j][aj]);
205
dens_pf = Dens[I][J];
207
wdens_pf = (-1.0)*Lagr[I][J];
209
hds_norm_pf = norm_pf;
215
/*--- A factor of 1/2 for correlated Lagrangians ---*/
216
if (strcmp(UserOptions.wfn,"SCF"))
219
/*----------------------------------
220
Half-derivative overap integrals
221
----------------------------------*/
223
tmp = 2.0*a1*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2, m2,
224
n2, jnorm, AB, PA, PB);
226
tmp -= l1*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2, m2,
227
n2, jnorm, AB, PA, PB);
228
HDS[coord_ax][I][J] += tmp*hds_norm_pf;
231
tmp = 2.0*a1*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2, m2,
232
n2, jnorm, AB, PA, PB);
234
tmp -= m1*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2, m2,
235
n2, jnorm, AB, PA, PB);
236
HDS[coord_ay][I][J] += tmp*hds_norm_pf;
239
tmp = 2.0*a1*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2, m2,
240
n2, jnorm, AB, PA, PB);
242
tmp -= n1*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2, m2,
243
n2, jnorm, AB, PA, PB);
244
HDS[coord_az][I][J] += tmp*hds_norm_pf;
248
tmp = 2.0*a2*overlap_int(a2, l2+1, m2, n2, jnorm, a1, l1, m1,
249
n1, inorm, AB, PB, PA);
251
tmp -= l2*overlap_int(a2, l2-1, m2, n2, jnorm, a1, l1, m1,
252
n1, inorm, AB, PB, PA);
253
HDS[coord_bx][J][I] += tmp*hds_norm_pf;
256
tmp = 2.0*a2*overlap_int(a2, l2, m2+1, n2, jnorm, a1, l1, m1,
257
n1, inorm, AB, PB, PA);
259
tmp -= m2*overlap_int(a2, l2, m2-1, n2, jnorm, a1, l1, m1,
260
n1, inorm, AB, PB, PA);
261
HDS[coord_by][J][I] += tmp*hds_norm_pf;
264
tmp = 2.0*a2*overlap_int(a2, l2, m2, n2+1, jnorm, a1, l1, m1,
265
n1, inorm, AB, PB, PA);
267
tmp -= n2*overlap_int(a2, l2, m2, n2-1, jnorm, a1, l1, m1,
268
n1, inorm, AB, PB, PA);
269
HDS[coord_bz][J][I] += tmp*hds_norm_pf;
272
/*----------------------------------
273
First derivative overap integrals
274
----------------------------------*/
276
tmp = 2.0*a1*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2, m2,
277
n2, jnorm, AB, PA, PB);
279
tmp -= l1*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2, m2,
280
n2, jnorm, AB, PA, PB);
281
S[coord_ax][I][J] += tmp*norm_pf;
284
tmp = 2.0*a1*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2, m2,
285
n2, jnorm, AB, PA, PB);
287
tmp -= m1*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2, m2,
288
n2, jnorm, AB, PA, PB);
289
S[coord_ay][I][J] += tmp*norm_pf;
292
tmp = 2.0*a1*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2, m2,
293
n2, jnorm, AB, PA, PB);
295
tmp -= n1*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2, m2,
296
n2, jnorm, AB, PA, PB);
297
S[coord_az][I][J] += tmp*norm_pf;
300
tmp = 2.0*a2*overlap_int(a1, l1, m1, n1, inorm, a2, l2+1, m2,
301
n2, jnorm, AB, PA, PB);
303
tmp -= l2*overlap_int(a1, l1, m1, n1, inorm, a2, l2-1, m2,
304
n2, jnorm, AB, PA, PB);
305
S[coord_bx][I][J] += tmp*norm_pf;
308
tmp = 2.0*a2*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2+1,
309
n2, jnorm, AB, PA, PB);
311
tmp -= m2*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2-1,
312
n2, jnorm, AB, PA, PB);
313
S[coord_by][I][J] += tmp*norm_pf;
316
tmp = 2.0*a2*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2,
317
n2+1, jnorm, AB, PA, PB);
319
tmp -= n2*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2,
320
n2-1, jnorm, AB, PA, PB);
321
S[coord_bz][I][J] += tmp*norm_pf;
324
/*------------------------------------------
325
First derivative kinetic energy integrals
326
------------------------------------------*/
328
tmp = 2.0*a1*ke_int(a1, l1+1, m1, n1, inorm, a2, l2, m2,
329
n2, jnorm, AB, PA, PB);
331
tmp -= l1*ke_int(a1, l1-1, m1, n1, inorm, a2, l2, m2,
332
n2, jnorm, AB, PA, PB);
333
F[coord_ax][I][J] += tmp*norm_pf;
336
tmp = 2.0*a1*ke_int(a1, l1, m1+1, n1, inorm, a2, l2, m2,
337
n2, jnorm, AB, PA, PB);
339
tmp -= m1*ke_int(a1, l1, m1-1, n1, inorm, a2, l2, m2,
340
n2, jnorm, AB, PA, PB);
341
F[coord_ay][I][J] += tmp*norm_pf;
344
tmp = 2.0*a1*ke_int(a1, l1, m1, n1+1, inorm, a2, l2, m2,
345
n2, jnorm, AB, PA, PB);
347
tmp -= n1*ke_int(a1, l1, m1, n1-1, inorm, a2, l2, m2,
348
n2, jnorm, AB, PA, PB);
349
F[coord_az][I][J] += tmp*norm_pf;
352
tmp = 2.0*a2*ke_int(a1, l1, m1, n1, inorm, a2, l2+1, m2,
353
n2, jnorm, AB, PA, PB);
355
tmp -= l2*ke_int(a1, l1, m1, n1, inorm, a2, l2-1, m2,
356
n2, jnorm, AB, PA, PB);
357
F[coord_bx][I][J] += tmp*norm_pf;
360
tmp = 2.0*a2*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2+1,
361
n2, jnorm, AB, PA, PB);
363
tmp -= m2*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2-1,
364
n2, jnorm, AB, PA, PB);
365
F[coord_by][I][J] += tmp*norm_pf;
368
tmp = 2.0*a2*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2,
369
n2+1, jnorm, AB, PA, PB);
371
tmp -= n2*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2,
372
n2-1, jnorm, AB, PA, PB);
373
F[coord_bz][I][J] += tmp*norm_pf;
376
/*------------------------------------
377
Second derivative overlap integrals
378
------------------------------------*/
380
s_int = overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2,
381
n2, jnorm, AB, PA, PB);
383
upuppfac = twozeta_a*twozeta_a;
385
tmp = upuppfac*overlap_int(a1, l1+2, m1, n1, inorm, a2, l2, m2,
386
n2, jnorm, AB, PA, PB);
387
updownpfac = l1 + 1; if (l1) updownpfac += l1;
388
tmp -= twozeta_a*updownpfac*s_int;
390
tmp += l1*(l1-1)*overlap_int(a1, l1-2, m1, n1, inorm, a2, l2, m2,
391
n2, jnorm, AB, PA, PB);
393
hess_ov[coord_ax][coord_ax] += tmp*wdens_pf;
396
tmp = upuppfac*overlap_int(a1, l1, m1+2, n1, inorm, a2, l2, m2,
397
n2, jnorm, AB, PA, PB);
398
updownpfac = m1 + 1; if (m1) updownpfac += m1;
399
tmp -= twozeta_a*updownpfac*s_int;
401
tmp += m1*(m1-1)*overlap_int(a1, l1, m1-2, n1, inorm, a2, l2, m2,
402
n2, jnorm, AB, PA, PB);
404
hess_ov[coord_ay][coord_ay] += tmp*wdens_pf;
407
tmp = upuppfac*overlap_int(a1, l1, m1, n1+2, inorm, a2, l2, m2,
408
n2, jnorm, AB, PA, PB);
409
updownpfac = n1 + 1; if (n1) updownpfac += n1;
410
tmp -= twozeta_a*updownpfac*s_int;
412
tmp += n1*(n1-1)*overlap_int(a1, l1, m1, n1-2, inorm, a2, l2, m2,
413
n2, jnorm, AB, PA, PB);
415
hess_ov[coord_az][coord_az] += tmp*wdens_pf;
417
/*--- d2/dAxdAy ---*/
418
tmp = upuppfac*overlap_int(a1, l1+1, m1+1, n1, inorm, a2, l2, m2,
419
n2, jnorm, AB, PA, PB);
421
tmp -= twozeta_a*l1*overlap_int(a1, l1-1, m1+1, n1, inorm, a2, l2, m2,
422
n2, jnorm, AB, PA, PB);
424
tmp -= twozeta_a*m1*overlap_int(a1, l1+1, m1-1, n1, inorm, a2, l2, m2,
425
n2, jnorm, AB, PA, PB);
426
if (l1 > 0 && m1 > 0)
427
tmp += l1*m1*overlap_int(a1, l1-1, m1-1, n1, inorm, a2, l2, m2,
428
n2, jnorm, AB, PA, PB);
430
hess_ov[coord_ax][coord_ay] += tmp*wdens_pf;
432
/*--- d2/dAxdAz ---*/
433
tmp = upuppfac*overlap_int(a1, l1+1, m1, n1+1, inorm, a2, l2, m2,
434
n2, jnorm, AB, PA, PB);
436
tmp -= twozeta_a*l1*overlap_int(a1, l1-1, m1, n1+1, inorm, a2, l2, m2,
437
n2, jnorm, AB, PA, PB);
439
tmp -= twozeta_a*n1*overlap_int(a1, l1+1, m1, n1-1, inorm, a2, l2, m2,
440
n2, jnorm, AB, PA, PB);
441
if (l1 > 0 && n1 > 0)
442
tmp += l1*n1*overlap_int(a1, l1-1, m1, n1-1, inorm, a2, l2, m2,
443
n2, jnorm, AB, PA, PB);
445
hess_ov[coord_ax][coord_az] += tmp*wdens_pf;
447
/*--- d2/dAydAz ---*/
448
tmp = upuppfac*overlap_int(a1, l1, m1+1, n1+1, inorm, a2, l2, m2,
449
n2, jnorm, AB, PA, PB);
451
tmp -= twozeta_a*m1*overlap_int(a1, l1, m1-1, n1+1, inorm, a2, l2, m2,
452
n2, jnorm, AB, PA, PB);
454
tmp -= twozeta_a*n1*overlap_int(a1, l1, m1+1, n1-1, inorm, a2, l2, m2,
455
n2, jnorm, AB, PA, PB);
456
if (m1 > 0 && n1 > 0)
457
tmp += m1*n1*overlap_int(a1, l1, m1-1, n1-1, inorm, a2, l2, m2,
458
n2, jnorm, AB, PA, PB);
460
hess_ov[coord_ay][coord_az] += tmp*wdens_pf;
463
upuppfac = twozeta_b*twozeta_b;
465
tmp = upuppfac*overlap_int(a1, l1, m1, n1, inorm, a2, l2+2, m2,
466
n2, jnorm, AB, PA, PB);
467
updownpfac = l2 + 1; if (l2) updownpfac += l2;
468
tmp -= twozeta_b*updownpfac*s_int;
470
tmp += l2*(l2-1)*overlap_int(a1, l1, m1, n1, inorm, a2, l2-2, m2,
471
n2, jnorm, AB, PA, PB);
473
hess_ov[coord_bx][coord_bx] += tmp*wdens_pf;
476
tmp = upuppfac*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2+2,
477
n2, jnorm, AB, PA, PB);
478
updownpfac = m2 + 1; if (m2) updownpfac += m2;
479
tmp -= twozeta_b*updownpfac*s_int;
481
tmp += m2*(m2-1)*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2-2,
482
n2, jnorm, AB, PA, PB);
484
hess_ov[coord_by][coord_by] += tmp*wdens_pf;
487
tmp = upuppfac*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2,
488
n2+2, jnorm, AB, PA, PB);
489
updownpfac = n2 + 1; if (n2) updownpfac += n2;
490
tmp -= twozeta_b*updownpfac*s_int;
492
tmp += n2*(n2-1)*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2,
493
n2-2, jnorm, AB, PA, PB);
495
hess_ov[coord_bz][coord_bz] += tmp*wdens_pf;
497
/*--- d2/dBxdBy ---*/
498
tmp = upuppfac*overlap_int(a1, l1, m1, n1, inorm, a2, l2+1, m2+1,
499
n2, jnorm, AB, PA, PB);
501
tmp -= twozeta_b*l2*overlap_int(a1, l1, m1, n1, inorm, a2, l2-1, m2+1,
502
n2, jnorm, AB, PA, PB);
504
tmp -= twozeta_b*m2*overlap_int(a1, l1, m1, n1, inorm, a2, l2+1, m2-1,
505
n2, jnorm, AB, PA, PB);
506
if (l2 > 0 && m2 > 0)
507
tmp += l2*m2*overlap_int(a1, l1, m1, n1, inorm, a2, l2-1, m2-1,
508
n2, jnorm, AB, PA, PB);
510
hess_ov[coord_bx][coord_by] += tmp*wdens_pf;
512
/*--- d2/dBxdBz ---*/
513
tmp = upuppfac*overlap_int(a1, l1, m1, n1, inorm, a2, l2+1, m2,
514
n2+1, jnorm, AB, PA, PB);
516
tmp -= twozeta_b*l2*overlap_int(a1, l1, m1, n1, inorm, a2, l2-1, m2,
517
n2+1, jnorm, AB, PA, PB);
519
tmp -= twozeta_b*n2*overlap_int(a1, l1, m1, n1, inorm, a2, l2+1, m2,
520
n2-1, jnorm, AB, PA, PB);
521
if (l2 > 0 && n2 > 0)
522
tmp += l2*n2*overlap_int(a1, l1, m1, n1, inorm, a2, l2-1, m2,
523
n2-1, jnorm, AB, PA, PB);
525
hess_ov[coord_bx][coord_bz] += tmp*wdens_pf;
527
/*--- d2/dBydBz ---*/
528
tmp = upuppfac*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2+1,
529
n2+1, jnorm, AB, PA, PB);
531
tmp -= twozeta_b*m2*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2-1,
532
n2+1, jnorm, AB, PA, PB);
534
tmp -= twozeta_b*n2*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2+1,
535
n2-1, jnorm, AB, PA, PB);
536
if (m2 > 0 && n2 > 0)
537
tmp += m2*n2*overlap_int(a1, l1, m1, n1, inorm, a2, l2, m2-1,
538
n2-1, jnorm, AB, PA, PB);
540
hess_ov[coord_by][coord_bz] += tmp*wdens_pf;
543
upuppfac = twozeta_a*twozeta_b;
544
/*--- d2/dAxdBx ---*/
545
tmp = upuppfac*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2+1, m2,
546
n2, jnorm, AB, PA, PB);
548
tmp -= twozeta_b*l1*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2+1, m2,
549
n2, jnorm, AB, PA, PB);
551
tmp -= twozeta_a*l2*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2-1, m2,
552
n2, jnorm, AB, PA, PB);
553
if (l1 > 0 && l2 > 0)
554
tmp += l1*l2*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2-1, m2,
555
n2, jnorm, AB, PA, PB);
556
if (coord_ax == coord_bx)
558
hess_ov[coord_ax][coord_bx] += tmp*wdens_pf;
560
/*--- d2/dAxdBy ---*/
561
tmp = upuppfac*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2, m2+1,
562
n2, jnorm, AB, PA, PB);
564
tmp -= twozeta_b*l1*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2, m2+1,
565
n2, jnorm, AB, PA, PB);
567
tmp -= twozeta_a*m2*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2, m2-1,
568
n2, jnorm, AB, PA, PB);
569
if (l1 > 0 && m2 > 0)
570
tmp += l1*m2*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2, m2-1,
571
n2, jnorm, AB, PA, PB);
573
hess_ov[coord_ax][coord_by] += tmp*wdens_pf;
575
/*--- d2/dAxdBz ---*/
576
tmp = upuppfac*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2, m2,
577
n2+1, jnorm, AB, PA, PB);
579
tmp -= twozeta_b*l1*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2, m2,
580
n2+1, jnorm, AB, PA, PB);
582
tmp -= twozeta_a*n2*overlap_int(a1, l1+1, m1, n1, inorm, a2, l2, m2,
583
n2-1, jnorm, AB, PA, PB);
584
if (l1 > 0 && n2 > 0)
585
tmp += l1*n2*overlap_int(a1, l1-1, m1, n1, inorm, a2, l2, m2,
586
n2-1, jnorm, AB, PA, PB);
588
hess_ov[coord_ax][coord_bz] += tmp*wdens_pf;
590
/*--- d2/dAydBx ---*/
591
tmp = upuppfac*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2+1, m2,
592
n2, jnorm, AB, PA, PB);
594
tmp -= twozeta_b*m1*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2+1, m2,
595
n2, jnorm, AB, PA, PB);
597
tmp -= twozeta_a*l2*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2-1, m2,
598
n2, jnorm, AB, PA, PB);
599
if (m1 > 0 && l2 > 0)
600
tmp += m1*l2*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2-1, m2,
601
n2, jnorm, AB, PA, PB);
603
hess_ov[coord_ay][coord_bx] += tmp*wdens_pf;
605
/*--- d2/dAydBy ---*/
606
tmp = upuppfac*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2, m2+1,
607
n2, jnorm, AB, PA, PB);
609
tmp -= twozeta_b*m1*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2, m2+1,
610
n2, jnorm, AB, PA, PB);
612
tmp -= twozeta_a*m2*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2, m2-1,
613
n2, jnorm, AB, PA, PB);
614
if (m1 > 0 && m2 > 0)
615
tmp += m1*m2*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2, m2-1,
616
n2, jnorm, AB, PA, PB);
617
if (coord_ay == coord_by)
619
hess_ov[coord_ay][coord_by] += tmp*wdens_pf;
621
/*--- d2/dAydBz ---*/
622
tmp = upuppfac*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2, m2,
623
n2+1, jnorm, AB, PA, PB);
625
tmp -= twozeta_b*m1*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2, m2,
626
n2+1, jnorm, AB, PA, PB);
628
tmp -= twozeta_a*n2*overlap_int(a1, l1, m1+1, n1, inorm, a2, l2, m2,
629
n2-1, jnorm, AB, PA, PB);
630
if (m1 > 0 && n2 > 0)
631
tmp += m1*n2*overlap_int(a1, l1, m1-1, n1, inorm, a2, l2, m2,
632
n2-1, jnorm, AB, PA, PB);
634
hess_ov[coord_ay][coord_bz] += tmp*wdens_pf;
636
/*--- d2/dAzdBx ---*/
637
tmp = upuppfac*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2+1, m2,
638
n2, jnorm, AB, PA, PB);
640
tmp -= twozeta_b*n1*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2+1, m2,
641
n2, jnorm, AB, PA, PB);
643
tmp -= twozeta_a*l2*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2-1, m2,
644
n2, jnorm, AB, PA, PB);
645
if (n1 > 0 && l2 > 0)
646
tmp += n1*l2*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2-1, m2,
647
n2, jnorm, AB, PA, PB);
649
hess_ov[coord_az][coord_bx] += tmp*wdens_pf;
651
/*--- d2/dAzdBy ---*/
652
tmp = upuppfac*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2, m2+1,
653
n2, jnorm, AB, PA, PB);
655
tmp -= twozeta_b*n1*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2, m2+1,
656
n2, jnorm, AB, PA, PB);
658
tmp -= twozeta_a*m2*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2, m2-1,
659
n2, jnorm, AB, PA, PB);
660
if (n1 > 0 && m2 > 0)
661
tmp += n1*m2*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2, m2-1,
662
n2, jnorm, AB, PA, PB);
664
hess_ov[coord_az][coord_by] += tmp*wdens_pf;
666
/*--- d2/dAzdBz ---*/
667
tmp = upuppfac*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2, m2,
668
n2+1, jnorm, AB, PA, PB);
670
tmp -= twozeta_b*n1*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2, m2,
671
n2+1, jnorm, AB, PA, PB);
673
tmp -= twozeta_a*n2*overlap_int(a1, l1, m1, n1+1, inorm, a2, l2, m2,
674
n2-1, jnorm, AB, PA, PB);
675
if (n1 > 0 && n2 > 0)
676
tmp += n1*n2*overlap_int(a1, l1, m1, n1-1, inorm, a2, l2, m2,
677
n2-1, jnorm, AB, PA, PB);
678
if (coord_az == coord_bz)
680
hess_ov[coord_az][coord_bz] += tmp*wdens_pf;
684
/*-------------------------------------------
685
Second derivative kinetic energy integrals
686
-------------------------------------------*/
688
t_int = ke_int(a1, l1, m1, n1, inorm, a2, l2, m2,
689
n2, jnorm, AB, PA, PB);
691
upuppfac = twozeta_a*twozeta_a;
693
tmp = upuppfac*ke_int(a1, l1+2, m1, n1, inorm, a2, l2, m2,
694
n2, jnorm, AB, PA, PB);
695
updownpfac = l1 + 1; if (l1) updownpfac += l1;
696
tmp -= twozeta_a*updownpfac*t_int;
698
tmp += l1*(l1-1)*ke_int(a1, l1-2, m1, n1, inorm, a2, l2, m2,
699
n2, jnorm, AB, PA, PB);
701
hess_oe[coord_ax][coord_ax] += tmp*dens_pf;
704
tmp = upuppfac*ke_int(a1, l1, m1+2, n1, inorm, a2, l2, m2,
705
n2, jnorm, AB, PA, PB);
706
updownpfac = m1 + 1; if (m1) updownpfac += m1;
707
tmp -= twozeta_a*updownpfac*t_int;
709
tmp += m1*(m1-1)*ke_int(a1, l1, m1-2, n1, inorm, a2, l2, m2,
710
n2, jnorm, AB, PA, PB);
712
hess_oe[coord_ay][coord_ay] += tmp*dens_pf;
715
tmp = upuppfac*ke_int(a1, l1, m1, n1+2, inorm, a2, l2, m2,
716
n2, jnorm, AB, PA, PB);
717
updownpfac = n1 + 1; if (n1) updownpfac += n1;
718
tmp -= twozeta_a*updownpfac*t_int;
720
tmp += n1*(n1-1)*ke_int(a1, l1, m1, n1-2, inorm, a2, l2, m2,
721
n2, jnorm, AB, PA, PB);
723
hess_oe[coord_az][coord_az] += tmp*dens_pf;
725
/*--- d2/dAxdAy ---*/
726
tmp = upuppfac*ke_int(a1, l1+1, m1+1, n1, inorm, a2, l2, m2,
727
n2, jnorm, AB, PA, PB);
729
tmp -= twozeta_a*l1*ke_int(a1, l1-1, m1+1, n1, inorm, a2, l2, m2,
730
n2, jnorm, AB, PA, PB);
732
tmp -= twozeta_a*m1*ke_int(a1, l1+1, m1-1, n1, inorm, a2, l2, m2,
733
n2, jnorm, AB, PA, PB);
734
if (l1 > 0 && m1 > 0)
735
tmp += l1*m1*ke_int(a1, l1-1, m1-1, n1, inorm, a2, l2, m2,
736
n2, jnorm, AB, PA, PB);
738
hess_oe[coord_ax][coord_ay] += tmp*dens_pf;
740
/*--- d2/dAxdAz ---*/
741
tmp = upuppfac*ke_int(a1, l1+1, m1, n1+1, inorm, a2, l2, m2,
742
n2, jnorm, AB, PA, PB);
744
tmp -= twozeta_a*l1*ke_int(a1, l1-1, m1, n1+1, inorm, a2, l2, m2,
745
n2, jnorm, AB, PA, PB);
747
tmp -= twozeta_a*n1*ke_int(a1, l1+1, m1, n1-1, inorm, a2, l2, m2,
748
n2, jnorm, AB, PA, PB);
749
if (l1 > 0 && n1 > 0)
750
tmp += l1*n1*ke_int(a1, l1-1, m1, n1-1, inorm, a2, l2, m2,
751
n2, jnorm, AB, PA, PB);
753
hess_oe[coord_ax][coord_az] += tmp*dens_pf;
755
/*--- d2/dAydAz ---*/
756
tmp = upuppfac*ke_int(a1, l1, m1+1, n1+1, inorm, a2, l2, m2,
757
n2, jnorm, AB, PA, PB);
759
tmp -= twozeta_a*m1*ke_int(a1, l1, m1-1, n1+1, inorm, a2, l2, m2,
760
n2, jnorm, AB, PA, PB);
762
tmp -= twozeta_a*n1*ke_int(a1, l1, m1+1, n1-1, inorm, a2, l2, m2,
763
n2, jnorm, AB, PA, PB);
764
if (m1 > 0 && n1 > 0)
765
tmp += m1*n1*ke_int(a1, l1, m1-1, n1-1, inorm, a2, l2, m2,
766
n2, jnorm, AB, PA, PB);
768
hess_oe[coord_ay][coord_az] += tmp*dens_pf;
770
upuppfac = twozeta_b*twozeta_b;
772
tmp = upuppfac*ke_int(a1, l1, m1, n1, inorm, a2, l2+2, m2,
773
n2, jnorm, AB, PA, PB);
774
updownpfac = l2 + 1; if (l2) updownpfac += l2;
775
tmp -= twozeta_b*updownpfac*t_int;
777
tmp += l2*(l2-1)*ke_int(a1, l1, m1, n1, inorm, a2, l2-2, m2,
778
n2, jnorm, AB, PA, PB);
780
hess_oe[coord_bx][coord_bx] += tmp*dens_pf;
783
tmp = upuppfac*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2+2,
784
n2, jnorm, AB, PA, PB);
785
updownpfac = m2 + 1; if (m2) updownpfac += m2;
786
tmp -= twozeta_b*updownpfac*t_int;
788
tmp += m2*(m2-1)*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2-2,
789
n2, jnorm, AB, PA, PB);
791
hess_oe[coord_by][coord_by] += tmp*dens_pf;
794
tmp = upuppfac*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2,
795
n2+2, jnorm, AB, PA, PB);
796
updownpfac = n2 + 1; if (n2) updownpfac += n2;
797
tmp -= twozeta_b*updownpfac*t_int;
799
tmp += n2*(n2-1)*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2,
800
n2-2, jnorm, AB, PA, PB);
802
hess_oe[coord_bz][coord_bz] += tmp*dens_pf;
804
/*--- d2/dBxdBy ---*/
805
tmp = upuppfac*ke_int(a1, l1, m1, n1, inorm, a2, l2+1, m2+1,
806
n2, jnorm, AB, PA, PB);
808
tmp -= twozeta_b*l2*ke_int(a1, l1, m1, n1, inorm, a2, l2-1, m2+1,
809
n2, jnorm, AB, PA, PB);
811
tmp -= twozeta_b*m2*ke_int(a1, l1, m1, n1, inorm, a2, l2+1, m2-1,
812
n2, jnorm, AB, PA, PB);
813
if (l2 > 0 && m2 > 0)
814
tmp += l2*m2*ke_int(a1, l1, m1, n1, inorm, a2, l2-1, m2-1,
815
n2, jnorm, AB, PA, PB);
817
hess_oe[coord_bx][coord_by] += tmp*dens_pf;
819
/*--- d2/dBxdBz ---*/
820
tmp = upuppfac*ke_int(a1, l1, m1, n1, inorm, a2, l2+1, m2,
821
n2+1, jnorm, AB, PA, PB);
823
tmp -= twozeta_b*l2*ke_int(a1, l1, m1, n1, inorm, a2, l2-1, m2,
824
n2+1, jnorm, AB, PA, PB);
826
tmp -= twozeta_b*n2*ke_int(a1, l1, m1, n1, inorm, a2, l2+1, m2,
827
n2-1, jnorm, AB, PA, PB);
828
if (l2 > 0 && n2 > 0)
829
tmp += l2*n2*ke_int(a1, l1, m1, n1, inorm, a2, l2-1, m2,
830
n2-1, jnorm, AB, PA, PB);
832
hess_oe[coord_bx][coord_bz] += tmp*dens_pf;
834
/*--- d2/dBydBz ---*/
835
tmp = upuppfac*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2+1,
836
n2+1, jnorm, AB, PA, PB);
838
tmp -= twozeta_b*m2*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2-1,
839
n2+1, jnorm, AB, PA, PB);
841
tmp -= twozeta_b*n2*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2+1,
842
n2-1, jnorm, AB, PA, PB);
843
if (m2 > 0 && n2 > 0)
844
tmp += m2*n2*ke_int(a1, l1, m1, n1, inorm, a2, l2, m2-1,
845
n2-1, jnorm, AB, PA, PB);
847
hess_oe[coord_by][coord_bz] += tmp*dens_pf;
850
upuppfac = twozeta_a*twozeta_b;
851
/*--- d2/dAxdBx ---*/
852
tmp = upuppfac*ke_int(a1, l1+1, m1, n1, inorm, a2, l2+1, m2,
853
n2, jnorm, AB, PA, PB);
855
tmp -= twozeta_b*l1*ke_int(a1, l1-1, m1, n1, inorm, a2, l2+1, m2,
856
n2, jnorm, AB, PA, PB);
858
tmp -= twozeta_a*l2*ke_int(a1, l1+1, m1, n1, inorm, a2, l2-1, m2,
859
n2, jnorm, AB, PA, PB);
860
if (l1 > 0 && l2 > 0)
861
tmp += l1*l2*ke_int(a1, l1-1, m1, n1, inorm, a2, l2-1, m2,
862
n2, jnorm, AB, PA, PB);
863
if (coord_ax == coord_bx)
865
hess_oe[coord_ax][coord_bx] += tmp*dens_pf;
867
/*--- d2/dAxdBy ---*/
868
tmp = upuppfac*ke_int(a1, l1+1, m1, n1, inorm, a2, l2, m2+1,
869
n2, jnorm, AB, PA, PB);
871
tmp -= twozeta_b*l1*ke_int(a1, l1-1, m1, n1, inorm, a2, l2, m2+1,
872
n2, jnorm, AB, PA, PB);
874
tmp -= twozeta_a*m2*ke_int(a1, l1+1, m1, n1, inorm, a2, l2, m2-1,
875
n2, jnorm, AB, PA, PB);
876
if (l1 > 0 && m2 > 0)
877
tmp += l1*m2*ke_int(a1, l1-1, m1, n1, inorm, a2, l2, m2-1,
878
n2, jnorm, AB, PA, PB);
880
hess_oe[coord_ax][coord_by] += tmp*dens_pf;
882
/*--- d2/dAxdBz ---*/
883
tmp = upuppfac*ke_int(a1, l1+1, m1, n1, inorm, a2, l2, m2,
884
n2+1, jnorm, AB, PA, PB);
886
tmp -= twozeta_b*l1*ke_int(a1, l1-1, m1, n1, inorm, a2, l2, m2,
887
n2+1, jnorm, AB, PA, PB);
889
tmp -= twozeta_a*n2*ke_int(a1, l1+1, m1, n1, inorm, a2, l2, m2,
890
n2-1, jnorm, AB, PA, PB);
891
if (l1 > 0 && n2 > 0)
892
tmp += l1*n2*ke_int(a1, l1-1, m1, n1, inorm, a2, l2, m2,
893
n2-1, jnorm, AB, PA, PB);
895
hess_oe[coord_ax][coord_bz] += tmp*dens_pf;
897
/*--- d2/dAydBx ---*/
898
tmp = upuppfac*ke_int(a1, l1, m1+1, n1, inorm, a2, l2+1, m2,
899
n2, jnorm, AB, PA, PB);
901
tmp -= twozeta_b*m1*ke_int(a1, l1, m1-1, n1, inorm, a2, l2+1, m2,
902
n2, jnorm, AB, PA, PB);
904
tmp -= twozeta_a*l2*ke_int(a1, l1, m1+1, n1, inorm, a2, l2-1, m2,
905
n2, jnorm, AB, PA, PB);
906
if (m1 > 0 && l2 > 0)
907
tmp += m1*l2*ke_int(a1, l1, m1-1, n1, inorm, a2, l2-1, m2,
908
n2, jnorm, AB, PA, PB);
910
hess_oe[coord_ay][coord_bx] += tmp*dens_pf;
912
/*--- d2/dAydBy ---*/
913
tmp = upuppfac*ke_int(a1, l1, m1+1, n1, inorm, a2, l2, m2+1,
914
n2, jnorm, AB, PA, PB);
916
tmp -= twozeta_b*m1*ke_int(a1, l1, m1-1, n1, inorm, a2, l2, m2+1,
917
n2, jnorm, AB, PA, PB);
919
tmp -= twozeta_a*m2*ke_int(a1, l1, m1+1, n1, inorm, a2, l2, m2-1,
920
n2, jnorm, AB, PA, PB);
921
if (m1 > 0 && m2 > 0)
922
tmp += m1*m2*ke_int(a1, l1, m1-1, n1, inorm, a2, l2, m2-1,
923
n2, jnorm, AB, PA, PB);
924
if (coord_ay == coord_by)
926
hess_oe[coord_ay][coord_by] += tmp*dens_pf;
928
/*--- d2/dAydBz ---*/
929
tmp = upuppfac*ke_int(a1, l1, m1+1, n1, inorm, a2, l2, m2,
930
n2+1, jnorm, AB, PA, PB);
932
tmp -= twozeta_b*m1*ke_int(a1, l1, m1-1, n1, inorm, a2, l2, m2,
933
n2+1, jnorm, AB, PA, PB);
935
tmp -= twozeta_a*n2*ke_int(a1, l1, m1+1, n1, inorm, a2, l2, m2,
936
n2-1, jnorm, AB, PA, PB);
937
if (m1 > 0 && n2 > 0)
938
tmp += m1*n2*ke_int(a1, l1, m1-1, n1, inorm, a2, l2, m2,
939
n2-1, jnorm, AB, PA, PB);
941
hess_oe[coord_ay][coord_bz] += tmp*dens_pf;
943
/*--- d2/dAzdBx ---*/
944
tmp = upuppfac*ke_int(a1, l1, m1, n1+1, inorm, a2, l2+1, m2,
945
n2, jnorm, AB, PA, PB);
947
tmp -= twozeta_b*n1*ke_int(a1, l1, m1, n1-1, inorm, a2, l2+1, m2,
948
n2, jnorm, AB, PA, PB);
950
tmp -= twozeta_a*l2*ke_int(a1, l1, m1, n1+1, inorm, a2, l2-1, m2,
951
n2, jnorm, AB, PA, PB);
952
if (n1 > 0 && l2 > 0)
953
tmp += n1*l2*ke_int(a1, l1, m1, n1-1, inorm, a2, l2-1, m2,
954
n2, jnorm, AB, PA, PB);
956
hess_oe[coord_az][coord_bx] += tmp*dens_pf;
958
/*--- d2/dAzdBy ---*/
959
tmp = upuppfac*ke_int(a1, l1, m1, n1+1, inorm, a2, l2, m2+1,
960
n2, jnorm, AB, PA, PB);
962
tmp -= twozeta_b*n1*ke_int(a1, l1, m1, n1-1, inorm, a2, l2, m2+1,
963
n2, jnorm, AB, PA, PB);
965
tmp -= twozeta_a*m2*ke_int(a1, l1, m1, n1+1, inorm, a2, l2, m2-1,
966
n2, jnorm, AB, PA, PB);
967
if (n1 > 0 && m2 > 0)
968
tmp += n1*m2*ke_int(a1, l1, m1, n1-1, inorm, a2, l2, m2-1,
969
n2, jnorm, AB, PA, PB);
971
hess_oe[coord_az][coord_by] += tmp*dens_pf;
973
/*--- d2/dAzdBz ---*/
974
tmp = upuppfac*ke_int(a1, l1, m1, n1+1, inorm, a2, l2, m2,
975
n2+1, jnorm, AB, PA, PB);
977
tmp -= twozeta_b*n1*ke_int(a1, l1, m1, n1-1, inorm, a2, l2, m2,
978
n2+1, jnorm, AB, PA, PB);
980
tmp -= twozeta_a*n2*ke_int(a1, l1, m1, n1+1, inorm, a2, l2, m2,
981
n2-1, jnorm, AB, PA, PB);
982
if (n1 > 0 && n2 > 0)
983
tmp += n1*n2*ke_int(a1, l1, m1, n1-1, inorm, a2, l2, m2,
984
n2-1, jnorm, AB, PA, PB);
985
if (coord_az == coord_bz)
987
hess_oe[coord_az][coord_bz] += tmp*dens_pf;
995
} /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
998
/*---------------------------------------------------------
999
First and second derivative nuclear attraction integrals
1000
---------------------------------------------------------*/
1001
for(atom=0;atom<Molecule.num_atoms;atom++) {
1003
coord_cy = coord_cx + 1;
1004
coord_cz = coord_cx + 2;
1005
PC.x = sp->P[i][j][0] - Molecule.centers[atom].x;
1006
PC.y = sp->P[i][j][1] - Molecule.centers[atom].y;
1007
PC.z = sp->P[i][j][2] - Molecule.centers[atom].z;
1008
AI_Deriv2_OSrecurs(AI0,AIX,AIY,AIZ,AIXX,AIXY,AIXZ,AIYY,AIYZ,AIZZ,PA,PB,PC,gam,am_i+deriv2_lvl,am_j+deriv2_lvl);
1010
/*--- create all am components of si ---*/
1012
for(ii = 0; ii <= am_i; ii++){
1014
for(jj = 0; jj <= ii; jj++){
1017
iind = n1*izm1 + m1*iym1 + l1*ixm1;
1020
/*--- create all am components of sj ---*/
1022
for(kk = 0; kk <= am_j; kk++){
1024
for(ll = 0; ll <= kk; ll++){
1027
jind = n2*jzm1 + m2*jym1 + l2*jxm1;
1030
if (si == sj && ai < aj)
1033
norm_pf = (GTOs.bf_norm[am_i][ai] * GTOs.bf_norm[am_j][aj]);
1034
normover_pf = norm_pf * over_pf;
1035
dens_pf = Dens[I][J];
1043
zdens_pf = over_pf * dens_pf * Molecule.centers[atom].Z_nuc;
1044
znormover_pf = normover_pf * Molecule.centers[atom].Z_nuc;
1046
tmp = 2.0*a1*AI0[iind+ixm1][jind][0];
1048
tmp -= l1*AI0[iind-ixm1][jind][0];
1049
F[coord_ax][I][J] -= tmp * Molecule.centers[atom].Z_nuc * (normover_pf);
1051
tmp = 2.0*a2*AI0[iind][jind+jxm1][0];
1053
tmp -= l2*AI0[iind][jind-jxm1][0];
1054
F[coord_bx][I][J] -= tmp * Molecule.centers[atom].Z_nuc * (normover_pf);
1056
F[coord_cx][I][J] -= AIX[iind][jind][0] * Molecule.centers[atom].Z_nuc * normover_pf;
1058
tmp = 2.0*a1*AI0[iind+iym1][jind][0];
1060
tmp -= m1*AI0[iind-iym1][jind][0];
1061
F[coord_ay][I][J] -= tmp * Molecule.centers[atom].Z_nuc * (normover_pf);
1063
tmp = 2.0*a2*AI0[iind][jind+jym1][0];
1065
tmp -= m2*AI0[iind][jind-jym1][0];
1066
F[coord_by][I][J] -= tmp * Molecule.centers[atom].Z_nuc * (normover_pf);
1068
F[coord_cy][I][J] -= AIY[iind][jind][0] * Molecule.centers[atom].Z_nuc * normover_pf;
1070
tmp = 2.0*a1*AI0[iind+izm1][jind][0];
1072
tmp -= n1*AI0[iind-izm1][jind][0];
1073
F[coord_az][I][J] -= tmp * Molecule.centers[atom].Z_nuc * (normover_pf);
1075
tmp = 2.0*a2*AI0[iind][jind+jzm1][0];
1077
tmp -= n2*AI0[iind][jind-jzm1][0];
1078
F[coord_bz][I][J] -= tmp * Molecule.centers[atom].Z_nuc * (normover_pf);
1080
F[coord_cz][I][J] -= AIZ[iind][jind][0] * Molecule.centers[atom].Z_nuc * normover_pf;
1082
#if INCLUDE_POTENTIAL
1083
v_int = AI0[iind][jind][0];
1085
upuppfac = twozeta_a*twozeta_a;
1087
tmp = upuppfac*AI0[iind+ixm1+ixm1][jind][0];
1088
updownpfac = l1 + 1; if (l1) updownpfac += l1;
1089
tmp -= twozeta_a*updownpfac*v_int;
1091
tmp += l1*(l1-1)*AI0[iind-ixm1-ixm1][jind][0];
1093
hess_oe[coord_ax][coord_ax] -= tmp * zdens_pf;
1096
tmp = upuppfac*AI0[iind+iym1+iym1][jind][0];
1097
updownpfac = m1 + 1; if (m1) updownpfac += m1;
1098
tmp -= twozeta_a*updownpfac*v_int;
1100
tmp += m1*(m1-1)*AI0[iind-iym1-iym1][jind][0];
1102
hess_oe[coord_ay][coord_ay] -= tmp * zdens_pf;
1105
tmp = upuppfac*AI0[iind+izm1+izm1][jind][0];
1106
updownpfac = n1 + 1; if (n1) updownpfac += n1;
1107
tmp -= twozeta_a*updownpfac*v_int;
1109
tmp += n1*(n1-1)*AI0[iind-izm1-izm1][jind][0];
1111
hess_oe[coord_az][coord_az] -= tmp * zdens_pf;
1113
/*--- d2/dAxdAy ---*/
1114
tmp = upuppfac*AI0[iind+ixm1+iym1][jind][0];
1116
tmp -= twozeta_a*l1*AI0[iind-ixm1+iym1][jind][0];
1118
tmp -= twozeta_a*m1*AI0[iind+ixm1-iym1][jind][0];
1119
if (l1 > 0 && m1 > 0)
1120
tmp += l1*m1*AI0[iind-ixm1-iym1][jind][0];
1122
hess_oe[coord_ax][coord_ay] -= tmp*zdens_pf;
1124
/*--- d2/dAxdAz ---*/
1125
tmp = upuppfac*AI0[iind+ixm1+izm1][jind][0];
1127
tmp -= twozeta_a*l1*AI0[iind-ixm1+izm1][jind][0];
1129
tmp -= twozeta_a*n1*AI0[iind+ixm1-izm1][jind][0];
1130
if (l1 > 0 && n1 > 0)
1131
tmp += l1*n1*AI0[iind-ixm1-izm1][jind][0];
1133
hess_oe[coord_ax][coord_az] -= tmp*zdens_pf;
1135
/*--- d2/dAydAz ---*/
1136
tmp = upuppfac*AI0[iind+iym1+izm1][jind][0];
1138
tmp -= twozeta_a*m1*AI0[iind-iym1+izm1][jind][0];
1140
tmp -= twozeta_a*n1*AI0[iind+iym1-izm1][jind][0];
1141
if (m1 > 0 && n1 > 0)
1142
tmp += m1*n1*AI0[iind-iym1-izm1][jind][0];
1144
hess_oe[coord_ay][coord_az] -= tmp*zdens_pf;
1147
upuppfac = twozeta_b*twozeta_b;
1149
tmp = upuppfac*AI0[iind][jind+jxm1+jxm1][0];
1150
updownpfac = l2 + 1; if (l2) updownpfac += l2;
1151
tmp -= twozeta_b*updownpfac*v_int;
1153
tmp += l2*(l2-1)*AI0[iind][jind-jxm1-jxm1][0];
1155
hess_oe[coord_bx][coord_bx] -= tmp * zdens_pf;
1158
tmp = upuppfac*AI0[iind][jind+jym1+jym1][0];
1159
updownpfac = m2 + 1; if (m2) updownpfac += m2;
1160
tmp -= twozeta_b*updownpfac*v_int;
1162
tmp += m2*(m2-1)*AI0[iind][jind-jym1-jym1][0];
1164
hess_oe[coord_by][coord_by] -= tmp * zdens_pf;
1167
tmp = upuppfac*AI0[iind][jind+jzm1+jzm1][0];
1168
updownpfac = n2 + 1; if (n2) updownpfac += n2;
1169
tmp -= twozeta_b*updownpfac*v_int;
1171
tmp += n2*(n2-1)*AI0[iind][jind-jzm1-jzm1][0];
1173
hess_oe[coord_bz][coord_bz] -= tmp * zdens_pf;
1175
/*--- d2/dBxdBy ---*/
1176
tmp = upuppfac*AI0[iind][jind+jxm1+jym1][0];
1178
tmp -= twozeta_b*l2*AI0[iind][jind-jxm1+jym1][0];
1180
tmp -= twozeta_b*m2*AI0[iind][jind+jxm1-jym1][0];
1181
if (l2 > 0 && m2 > 0)
1182
tmp += l2*m2*AI0[iind][jind-jxm1-jym1][0];
1184
hess_oe[coord_bx][coord_by] -= tmp*zdens_pf;
1186
/*--- d2/dBxdBz ---*/
1187
tmp = upuppfac*AI0[iind][jind+jxm1+jzm1][0];
1189
tmp -= twozeta_b*l2*AI0[iind][jind-jxm1+jzm1][0];
1191
tmp -= twozeta_b*n2*AI0[iind][jind+jxm1-jzm1][0];
1192
if (l2 > 0 && n2 > 0)
1193
tmp += l2*n2*AI0[iind][jind-jxm1-jzm1][0];
1195
hess_oe[coord_bx][coord_bz] -= tmp*zdens_pf;
1197
/*--- d2/dBydBz ---*/
1198
tmp = upuppfac*AI0[iind][jind+jym1+jzm1][0];
1200
tmp -= twozeta_b*m2*AI0[iind][jind-jym1+jzm1][0];
1202
tmp -= twozeta_b*n2*AI0[iind][jind+jym1-jzm1][0];
1203
if (m2 > 0 && n2 > 0)
1204
tmp += m2*n2*AI0[iind][jind-jym1-jzm1][0];
1206
hess_oe[coord_by][coord_bz] -= tmp*zdens_pf;
1209
upuppfac = twozeta_a*twozeta_b;
1210
/*--- d2/dAxdBx ---*/
1211
tmp = upuppfac*AI0[iind+ixm1][jind+jxm1][0];
1213
tmp -= twozeta_b*l1*AI0[iind-ixm1][jind+jxm1][0];
1215
tmp -= twozeta_a*l2*AI0[iind+ixm1][jind-jxm1][0];
1216
if (l1 > 0 && l2 > 0)
1217
tmp += l1*l2*AI0[iind-ixm1][jind-jxm1][0];
1218
if (coord_ax == coord_bx)
1220
hess_oe[coord_ax][coord_bx] -= tmp*zdens_pf;
1222
/*--- d2/dAxdBy ---*/
1223
tmp = upuppfac*AI0[iind+ixm1][jind+jym1][0];
1225
tmp -= twozeta_b*l1*AI0[iind-ixm1][jind+jym1][0];
1227
tmp -= twozeta_a*m2*AI0[iind+ixm1][jind-jym1][0];
1228
if (l1 > 0 && m2 > 0)
1229
tmp += l1*m2*AI0[iind-ixm1][jind-jym1][0];
1230
hess_oe[coord_ax][coord_by] -= tmp*zdens_pf;
1232
/*--- d2/dAxdBz ---*/
1233
tmp = upuppfac*AI0[iind+ixm1][jind+jzm1][0];
1235
tmp -= twozeta_b*l1*AI0[iind-ixm1][jind+jzm1][0];
1237
tmp -= twozeta_a*n2*AI0[iind+ixm1][jind-jzm1][0];
1238
if (l1 > 0 && n2 > 0)
1239
tmp += l1*n2*AI0[iind-ixm1][jind-jzm1][0];
1240
hess_oe[coord_ax][coord_bz] -= tmp*zdens_pf;
1242
/*--- d2/dAydBx ---*/
1243
tmp = upuppfac*AI0[iind+iym1][jind+jxm1][0];
1245
tmp -= twozeta_b*m1*AI0[iind-iym1][jind+jxm1][0];
1247
tmp -= twozeta_a*l2*AI0[iind+iym1][jind-jxm1][0];
1248
if (m1 > 0 && l2 > 0)
1249
tmp += m1*l2*AI0[iind-iym1][jind-jxm1][0];
1250
hess_oe[coord_ay][coord_bx] -= tmp*zdens_pf;
1252
/*--- d2/dAydBy ---*/
1253
tmp = upuppfac*AI0[iind+iym1][jind+jym1][0];
1255
tmp -= twozeta_b*m1*AI0[iind-iym1][jind+jym1][0];
1257
tmp -= twozeta_a*m2*AI0[iind+iym1][jind-jym1][0];
1258
if (m1 > 0 && m2 > 0)
1259
tmp += m1*m2*AI0[iind-iym1][jind-jym1][0];
1260
if (coord_ay == coord_by)
1262
hess_oe[coord_ay][coord_by] -= tmp*zdens_pf;
1264
/*--- d2/dAydBz ---*/
1265
tmp = upuppfac*AI0[iind+iym1][jind+jzm1][0];
1267
tmp -= twozeta_b*m1*AI0[iind-iym1][jind+jzm1][0];
1269
tmp -= twozeta_a*n2*AI0[iind+iym1][jind-jzm1][0];
1270
if (m1 > 0 && n2 > 0)
1271
tmp += m1*n2*AI0[iind-iym1][jind-jzm1][0];
1272
hess_oe[coord_ay][coord_bz] -= tmp*zdens_pf;
1274
/*--- d2/dAzdBx ---*/
1275
tmp = upuppfac*AI0[iind+izm1][jind+jxm1][0];
1277
tmp -= twozeta_b*n1*AI0[iind-izm1][jind+jxm1][0];
1279
tmp -= twozeta_a*l2*AI0[iind+izm1][jind-jxm1][0];
1280
if (n1 > 0 && l2 > 0)
1281
tmp += n1*l2*AI0[iind-izm1][jind-jxm1][0];
1282
hess_oe[coord_az][coord_bx] -= tmp*zdens_pf;
1284
/*--- d2/dAzdBy ---*/
1285
tmp = upuppfac*AI0[iind+izm1][jind+jym1][0];
1287
tmp -= twozeta_b*n1*AI0[iind-izm1][jind+jym1][0];
1289
tmp -= twozeta_a*m2*AI0[iind+izm1][jind-jym1][0];
1290
if (n1 > 0 && m2 > 0)
1291
tmp += n1*m2*AI0[iind-izm1][jind-jym1][0];
1292
hess_oe[coord_az][coord_by] -= tmp*zdens_pf;
1294
/*--- d2/dAzdBz ---*/
1295
tmp = upuppfac*AI0[iind+izm1][jind+jzm1][0];
1297
tmp -= twozeta_b*n1*AI0[iind-izm1][jind+jzm1][0];
1299
tmp -= twozeta_a*n2*AI0[iind+izm1][jind-jzm1][0];
1300
if (n1 > 0 && n2 > 0)
1301
tmp += n1*n2*AI0[iind-izm1][jind-jzm1][0];
1302
if (coord_az == coord_bz)
1304
hess_oe[coord_az][coord_bz] -= tmp*zdens_pf;
1308
/*--- d2/dCxdAx ---*/
1309
tmp = twozeta_a*AIX[iind+ixm1][jind][0];
1311
tmp -= l1*AIX[iind-ixm1][jind][0];
1312
if (coord_ax == coord_cx) tmp *= 2.0;
1313
hess_oe[coord_ax][coord_cx] -= tmp*zdens_pf;
1315
/*--- d2/dCxdAy ---*/
1316
tmp = twozeta_a*AIX[iind+iym1][jind][0];
1318
tmp -= m1*AIX[iind-iym1][jind][0];
1319
hess_oe[coord_ay][coord_cx] -= tmp*zdens_pf;
1321
/*--- d2/dCxdAz ---*/
1322
tmp = twozeta_a*AIX[iind+izm1][jind][0];
1324
tmp -= n1*AIX[iind-izm1][jind][0];
1325
hess_oe[coord_az][coord_cx] -= tmp*zdens_pf;
1327
/*--- d2/dCydAx ---*/
1328
tmp = twozeta_a*AIY[iind+ixm1][jind][0];
1330
tmp -= l1*AIY[iind-ixm1][jind][0];
1331
hess_oe[coord_ax][coord_cy] -= tmp*zdens_pf;
1333
/*--- d2/dCydAy ---*/
1334
tmp = twozeta_a*AIY[iind+iym1][jind][0];
1336
tmp -= m1*AIY[iind-iym1][jind][0];
1337
if (coord_ay == coord_cy) tmp *= 2.0;
1338
hess_oe[coord_ay][coord_cy] -= tmp*zdens_pf;
1340
/*--- d2/dCydAz ---*/
1341
tmp = twozeta_a*AIY[iind+izm1][jind][0];
1343
tmp -= n1*AIY[iind-izm1][jind][0];
1344
hess_oe[coord_az][coord_cy] -= tmp*zdens_pf;
1346
/*--- d2/dCzdAx ---*/
1347
tmp = twozeta_a*AIZ[iind+ixm1][jind][0];
1349
tmp -= l1*AIZ[iind-ixm1][jind][0];
1350
hess_oe[coord_ax][coord_cz] -= tmp*zdens_pf;
1352
/*--- d2/dCzdAy ---*/
1353
tmp = twozeta_a*AIZ[iind+iym1][jind][0];
1355
tmp -= m1*AIZ[iind-iym1][jind][0];
1356
hess_oe[coord_ay][coord_cz] -= tmp*zdens_pf;
1358
/*--- d2/dCzdAz ---*/
1359
tmp = twozeta_a*AIZ[iind+izm1][jind][0];
1361
tmp -= n1*AIZ[iind-izm1][jind][0];
1362
if (coord_az == coord_cz) tmp *= 2.0;
1363
hess_oe[coord_az][coord_cz] -= tmp*zdens_pf;
1366
/*--- d2/dCxdBx ---*/
1367
tmp = twozeta_b*AIX[iind][jind+jxm1][0];
1369
tmp -= l2*AIX[iind][jind-jxm1][0];
1370
if (coord_bx == coord_cx) tmp *= 2.0;
1371
hess_oe[coord_bx][coord_cx] -= tmp*zdens_pf;
1373
/*--- d2/dCxdBy ---*/
1374
tmp = twozeta_b*AIX[iind][jind+jym1][0];
1376
tmp -= m2*AIX[iind][jind-jym1][0];
1377
hess_oe[coord_by][coord_cx] -= tmp*zdens_pf;
1379
/*--- d2/dCxdBz ---*/
1380
tmp = twozeta_b*AIX[iind][jind+jzm1][0];
1382
tmp -= n2*AIX[iind][jind-jzm1][0];
1383
hess_oe[coord_bz][coord_cx] -= tmp*zdens_pf;
1385
/*--- d2/dCydBx ---*/
1386
tmp = twozeta_b*AIY[iind][jind+jxm1][0];
1388
tmp -= l2*AIY[iind][jind-jxm1][0];
1389
hess_oe[coord_bx][coord_cy] -= tmp*zdens_pf;
1391
/*--- d2/dCydBy ---*/
1392
tmp = twozeta_b*AIY[iind][jind+jym1][0];
1394
tmp -= m2*AIY[iind][jind-jym1][0];
1395
if (coord_by == coord_cy) tmp *= 2.0;
1396
hess_oe[coord_by][coord_cy] -= tmp*zdens_pf;
1398
/*--- d2/dCydBz ---*/
1399
tmp = twozeta_b*AIY[iind][jind+jzm1][0];
1401
tmp -= n2*AIY[iind][jind-jzm1][0];
1402
hess_oe[coord_bz][coord_cy] -= tmp*zdens_pf;
1404
/*--- d2/dCzdBx ---*/
1405
tmp = twozeta_b*AIZ[iind][jind+jxm1][0];
1407
tmp -= l2*AIZ[iind][jind-jxm1][0];
1408
hess_oe[coord_bx][coord_cz] -= tmp*zdens_pf;
1410
/*--- d2/dCzdBy ---*/
1411
tmp = twozeta_b*AIZ[iind][jind+jym1][0];
1413
tmp -= m2*AIZ[iind][jind-jym1][0];
1414
hess_oe[coord_by][coord_cz] -= tmp*zdens_pf;
1416
/*--- d2/dCzdBz ---*/
1417
tmp = twozeta_b*AIZ[iind][jind+jzm1][0];
1419
tmp -= n2*AIZ[iind][jind-jzm1][0];
1420
if (coord_bz == coord_cz) tmp *= 2.0;
1421
hess_oe[coord_bz][coord_cz] -= tmp*zdens_pf;
1425
/*---- d2/dCx2 ---*/
1426
hess_oe[coord_cx][coord_cx] -= AIXX[iind][jind][0]*zdens_pf;
1427
/*---- d2/dCy2 ---*/
1428
hess_oe[coord_cy][coord_cy] -= AIYY[iind][jind][0]*zdens_pf;
1429
/*---- d2/dCz2 ---*/
1430
hess_oe[coord_cz][coord_cz] -= AIZZ[iind][jind][0]*zdens_pf;
1431
/*---- d2/dCxdCy ---*/
1432
hess_oe[coord_cx][coord_cy] -= AIXY[iind][jind][0]*zdens_pf;
1433
/*---- d2/dCxdCz ---*/
1434
hess_oe[coord_cx][coord_cz] -= AIXZ[iind][jind][0]*zdens_pf;
1435
/*---- d2/dCydCz ---*/
1436
hess_oe[coord_cy][coord_cz] -= AIYZ[iind][jind][0]*zdens_pf;
1444
} /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
1447
} /*--- end primitive contraction ---*/
1453
symmetrize_hessian(hess_ov);
1454
symmetrize_hessian(hess_oe);
1455
if (UserOptions.print_lvl >= PRINT_OEDERIV) {
1456
print_atommat("Overlap component of the molecular Hessian (a.u.)",hess_ov);
1457
print_atommat("Core Hamiltonian component of the molecular Hessian (a.u.)",hess_oe);
1460
add_mat(Hess,hess_ov,Hess,Molecule.num_atoms*3,Molecule.num_atoms*3);
1461
add_mat(Hess,hess_oe,Hess,Molecule.num_atoms*3,Molecule.num_atoms*3);
1463
/*--- flush it all away ---*/
1465
free_box(AI0,indmax,indmax);
1466
free_box(AIX,indmax,indmax);
1467
free_box(AIY,indmax,indmax);
1468
free_box(AIZ,indmax,indmax);
1469
free_box(AIXX,indmax,indmax);
1470
free_box(AIXY,indmax,indmax);
1471
free_box(AIXZ,indmax,indmax);
1472
free_box(AIYY,indmax,indmax);
1473
free_box(AIYZ,indmax,indmax);
1474
free_box(AIZZ,indmax,indmax);
1475
free_block(hess_ov);
1476
free_block(hess_oe);
1480
/*!-----------------------------------------------
1481
This computes overlap of 2 primitive gaussians
1482
-----------------------------------------------*/
1483
double overlap_int(double a1, int l1, int m1, int n1, double norm1,
1484
double a2, int l2, int m2, int n2, double norm2,
1485
struct coordinates AB,
1486
struct coordinates PA,
1487
struct coordinates PB)
1490
int imax, jmax, kmax;
1496
double tval, tval1, tval2 ;
1500
AB2 = AB.x*AB.x+AB.y*AB.y+AB.z*AB.z;
1502
norm_fact = norm1*norm2;
1507
for (i = 0; i<= imax; i++){
1508
tval = f_n(i*2, l1, l2, PA.x, PB.x);
1509
tval2 = int_pow(tval1, i);
1510
Ix += tval*(num_ser[i])/(tval2);
1514
for (j = 0; j<= jmax; j++){
1515
tval = f_n(j*2, m1, m2, PA.y, PB.y);
1516
tval2 = int_pow(tval1, j);
1517
Iy += tval*num_ser[j]/(tval2);
1521
for (k = 0; k<= kmax; k++){
1522
tval = f_n(k*2, n1, n2, PA.z, PB.z);
1523
tval2 = int_pow(tval1, k);
1524
Iz += tval*num_ser[k]/(tval2);
1527
I=exp(-1*a1*a2*AB2/gam)*Ix*Iy*Iz*sqrt(M_PI/gam)*(M_PI/gam);
1533
/*!-----------------------------------------------------------------------------
1534
This computes matrix element of kinetic energy between 2 primitive gaussians
1535
-----------------------------------------------------------------------------*/
1536
double ke_int(double a1, int l1, int m1, int n1, double norm1,
1537
double a2, int l2, int m2, int n2, double norm2,
1538
struct coordinates AB,
1539
struct coordinates PA,
1540
struct coordinates PB)
1542
int localprint = 0 ;
1544
double I1 = 0.0, I2 = 0.0, I3 = 0.0, I4 = 0.0 ;
1548
I2 = overlap_int(a1, l1+1, m1, n1, norm1, a2, l2+1, m2, n2, norm2,
1550
I1 = overlap_int(a1, l1-1, m1, n1, norm1, a2, l2-1, m2, n2, norm2,
1552
I3 = overlap_int(a1, l1+1, m1, n1, norm1, a2, l2-1, m2, n2, norm2,
1554
I4 = overlap_int(a1, l1-1, m1, n1, norm1, a2, l2+1, m2, n2, norm2,
1557
Ix = (l1*l2*I1/2.0 + 2*a1*a2*I2 - a1*l2*I3 - a2*l1*I4);
1564
I2 = overlap_int(a1, l1, m1+1, n1, norm1, a2, l2, m2+1, n2, norm2,
1566
I1 = overlap_int(a1, l1, m1-1, n1, norm1, a2, l2, m2-1, n2, norm2,
1568
I3 = overlap_int(a1, l1, m1+1, n1, norm1, a2, l2, m2-1, n2, norm2,
1570
I4 = overlap_int(a1, l1, m1-1, n1, norm1, a2, l2, m2+1, n2, norm2,
1573
Iy = (m1*m2*I1/2.0 + 2*a1*a2*I2 - a1*m2*I3 - a2*m1*I4);
1580
I2 = overlap_int(a1, l1, m1, n1+1, norm1, a2, l2, m2, n2+1, norm2,
1582
I1 = overlap_int(a1, l1, m1, n1-1, norm1, a2, l2, m2, n2-1, norm2,
1584
I3 = overlap_int(a1, l1, m1, n1+1, norm1, a2, l2, m2, n2-1, norm2,
1586
I4 = overlap_int(a1, l1, m1, n1-1, norm1, a2, l2, m2, n2+1, norm2,
1589
Iz = (n1*n2*I1/2.0 + 2*a1*a2*I2 - a1*n2*I3 - a2*n1*I4);
1595
double f_n(int k, int l1, int l2, double A, double B)
1600
double tmp1, tmp2, tmp3, tmp4 ;
1602
for(i=0; i<=l1; i++){
1603
for(j=0; j<=l2; j++){
1604
tmp1 = tmp2 = tmp3 = tmp4 = 0.0 ;
1606
tmp1 = int_pow(A, (l1-i));
1607
tmp2 = int_pow(B, (l2-j));
1610
tval = tmp1 * tmp2 * tmp3 * tmp4 ;
1618
double int_pow(double a, int p)
1623
for(i=0; i<p; i++) b = b*a;
added
removed
src/bin/cints/Default_Deriv2/oe_deriv2.cc
