3
#include <libipv1/ip_lib.h>
4
#include <libciomr/libciomr.h>
6
#include "globaldefs.h"
10
** form_appx_diag_mo_hess
12
** Calculates the approximate diagonal MO Hessian according to
13
** G. Chaban, M. W. Schmidt, and M. S. Gordon, Theor. Chem. Acc.,
16
** Also uses active-active formulae from G. Chaban, personal communication
17
** (derived from H. J. Aa. Jensen and H. Agren, Chem. Phys. 104, 229 (1986))
19
** I am assuming that the pairs (p,q) are always given such that p>=q
24
** Updated with active-active parts, March 2004
26
void form_appx_diag_mo_hess(int npairs, int *ppair, int *qpair,
27
double *F_core, double *tei, double **opdm, double *tpdm, double *F_act,
28
int firstact, int lastact, double *hess)
31
int pair, p, q, pq, pp, qq;
33
int i,ii,a,aa,t,tt,u,tu,v,w,vw,uv,tuvw;
34
int qv,puvw,quvw,pv,pu,qu,pupv,quqv,ppuv,qquv,pvqu,puqv,pquv;
37
fprintf(outfile, "Forming approximate diagonal orbital Hessian\n");
39
/* loop over the independent pairs */
40
for (pair=0; pair<npairs; pair++) {
47
/* H_{ai,ai}, i.e., inactive virt/inactive occ */
48
if (p >= lastact && q < firstact) {
49
hess[pair] = 4.0 * (F_core[pp] + F_act[pp] - F_core[qq] - F_act[qq]);
52
/* H_{at,at}, i.e., inactive virt with active orb */
53
else if (p >= lastact && q >= firstact) {
57
hess[pair] = 2.0 * opdm[t][t] * (F_core[aa] + F_act[aa]);
60
for (u=firstact; u<lastact; u++) {
62
value += opdm[t][u] * F_core[tu];
64
hess[pair] -= 2.0 * value;
67
for (u=firstact; u<lastact; u++) {
70
/* loop over active v,w ... later restrict to symmetry allowed */
71
for (v=firstact; v<lastact; v++) {
72
for (w=firstact; w<lastact; w++) {
75
value += tpdm[tuvw] * tei[tuvw];
80
hess[pair] -= 2.0 * value;
84
/* H_{ti,ti}, i.e., active orb with inactive occ */
85
else if (p >= firstact && q < firstact) {
90
hess[pair] = 2.0 * opdm[t][t] * (F_core[ii] + F_act[ii]);
91
hess[pair] += 4.0 * (F_core[tt] + F_act[tt] - F_core[ii] - F_act[ii]);
94
for (u=firstact; u<lastact; u++) {
96
value += opdm[t][u] * F_core[tu];
98
hess[pair] -= 2.0 * value;
101
for (u=firstact; u<lastact; u++) {
104
/* loop over active v,w ... later restrict to symmetry allowed */
105
for (v=firstact; v<lastact; v++) {
106
for (w=firstact; w<lastact; w++) {
109
value += tpdm[tuvw] * tei[tuvw];
114
hess[pair] -= 2.0 * value;
118
H_{t1t2,t1t2}, i.e., active-active (happens for MCSCF, not CASSCF)
119
expression given by personal communication from Galina Chaban,
120
derived from H. J. Aa. Jensen, H. Agren, Chem. Phys. 104, 229 (1982),
123
else if (p >= firstact && q < lastact) {
124
hess[pair] = 2.0 * ( opdm[q][q] * F_core[pp]
125
+ opdm[p][p] * F_core[qq]
126
- 2.0 * opdm[p][q] * F_core[pq] );
128
for (u=firstact; u<lastact; u++) {
131
value += opdm[p][u] * F_core[pu] + opdm[q][u] * F_core[qu];
132
/* should be able to reduce work below by factor of 2 */
133
for (w=firstact; w<lastact; w++) {
134
for (v=firstact; v<lastact; v++) {
138
value += tpdm[puvw] * tei[puvw];
139
value += tpdm[quvw] * tei[quvw];
143
hess[pair] -= value * 2.0;
146
for (u=firstact; u<lastact; u++) {
149
for (v=firstact; v<lastact; v++) {
155
value += 2.0 * tpdm[pupv] * tei[quqv];
156
value += 2.0 * tpdm[quqv] * tei[pupv];
159
value += tpdm[ppuv] * tei[qquv] + tpdm[qquv] * tei[ppuv];
162
value -= 4.0 * tpdm[pvqu] * tei[puqv];
164
value -= 2.0 * tpdm[pquv] * tei[pquv];
167
hess[pair] += 2.0 * value;
168
if (Params.scale_act_act != 1.0) hess[pair] *= Params.scale_act_act;
169
} /* end case H_{t1t2,t1t2} */
174
"(form_diag_mo_hess): Error, unrecognized class of indep pair\n");
177
} /* end loop over pairs */
186
** Calculates the exact diagonal MO Hessian assuming CASSCF
188
** See Galina Chaban et al., Theor Chem Acc (1997) 97:88-95
189
** and references therein for mathematics of Newton-Raphson
190
** approach to MCSCF/CASSCF
192
** Supplemented with active-active parts from G. Chaban, personal
195
** I am assuming that the pairs (p,q) are always given such that p>=q
200
** Active-active parts added by C. D. Sherrill, March 2004
202
void form_diag_mo_hess(int npairs, int *ppair, int *qpair,
203
double *F_core, double *tei, double **opdm, double *tpdm, double *F_act,
204
int firstact, int lastact, double *hess)
207
int pair, p, q, pq, pp, qq, pqpq, ppqq;
208
int i,ii,a,aa,t,tt,u,tu,v,w,vw,tuvw;
209
int au, uv, av, tv, ttuv, aauv, tvtu, avau;
210
int ui, vi, ti, uvii, uivi, uiti, tuii;
211
int qv,puvw,quvw,pv,pu,qu,pupv,quqv,ppuv,qquv,pvqu,puqv,pquv;
215
fprintf(outfile, "Forming diagonal orbital Hessian\n");
217
/* loop over the independent pairs */
218
for (pair=0; pair<npairs; pair++) {
224
pqpq = ioff[pq] + pq;
225
ppqq = ioff[pp] + qq;
227
/* H_{ai,ai}, i.e., inactive virt/inactive occ */
228
if (p >= lastact && q < firstact) {
229
hess[pair] = 4.0 * (F_core[pp] + F_act[pp] - F_core[qq] - F_act[qq]
230
+ 3.0 * tei[pqpq] - tei[ppqq]);
233
/* H_{at,at}, i.e., inactive virt with active orb */
234
else if (p >= lastact && q >= firstact) {
239
hess[pair] = 2.0 * opdm[t][t] * (F_core[aa]);
242
for (u=firstact; u<lastact; u++) {
244
value += opdm[t][u] * F_core[tu];
246
hess[pair] -= 2.0 * value;
249
for (u=firstact; u<lastact; u++) {
252
/* loop over active v,w ... later restrict to symmetry allowed */
253
for (v=firstact; v<lastact; v++) {
254
for (w=firstact; w<lastact; w++) {
257
value += tpdm[tuvw] * tei[tuvw];
262
hess[pair] -= 2.0 * value;
265
for (u=firstact; u<lastact; u++) {
268
for (v=firstact; v<lastact; v++) {
276
value += ( ( tpdm[ttuv] * tei[aauv] ) +
277
( 2.0 * tpdm[tvtu] * tei[avau] ) );
280
hess[pair] += 2.0 * value;
284
/* H_{ti,ti}, i.e., active orb with inactive occ */
285
else if (p >= firstact && q < firstact) {
291
hess[pair] = 2.0 * opdm[t][t] * (F_core[ii]);
292
hess[pair] += 4.0 * (F_core[tt] + F_act[tt] - F_core[ii] - F_act[ii]);
295
for (u=firstact; u<lastact; u++) {
297
value += opdm[t][u] * F_core[tu];
299
hess[pair] -= 2.0 * value;
302
for (u=firstact; u<lastact; u++) {
305
/* loop over active v,w ... later restrict to symmetry allowed */
306
for (v=firstact; v<lastact; v++) {
307
for (w=firstact; w<lastact; w++) {
310
value += tpdm[tuvw] * tei[tuvw];
315
hess[pair] -= 2.0 * value;
318
for (u=firstact; u<lastact; u++) {
321
for (v=firstact; v<lastact; v++) {
329
value += ( ( tpdm[ttuv] * tei[uvii] ) +
330
( 2.0 * tpdm[tvtu] * tei[uivi] ) );
333
hess[pair] += 2.0 * value;
336
for (u=firstact; u<lastact; u++) {
341
/* Create delta(t,u) */
347
value += ( ( delta - opdm[t][u] ) *
348
( ( 3.0 * tei[uiti] ) - tei[tuii] ) );
350
hess[pair] += 4.0 * value;
354
H_{t1t2,t1t2}, i.e., active-active (happens for MCSCF, not CASSCF)
355
expression given by personal communication from Galina Chaban,
356
derived from H. J. Aa. Jensen, H. Agren, Chem. Phys. 104, 229 (1982),
359
else if (p >= firstact && q < lastact) {
360
hess[pair] = 2.0 * ( opdm[q][q] * F_core[pp]
361
+ opdm[p][p] * F_core[qq]
362
- 2.0 * opdm[p][q] * F_core[pq] );
364
for (u=firstact; u<lastact; u++) {
367
value += opdm[p][u] * F_core[pu] + opdm[q][u] * F_core[qu];
368
/* should be able to reduce work below by factor of 2 */
369
for (w=firstact; w<lastact; w++) {
370
for (v=firstact; v<lastact; v++) {
374
value += tpdm[puvw] * tei[puvw];
375
value += tpdm[quvw] * tei[quvw];
379
hess[pair] -= value * 2.0;
382
for (u=firstact; u<lastact; u++) {
385
for (v=firstact; v<lastact; v++) {
391
value += 2.0 * tpdm[pupv] * tei[quqv];
392
value += 2.0 * tpdm[quqv] * tei[pupv];
395
value += tpdm[ppuv] * tei[qquv] + tpdm[qquv] * tei[ppuv];
398
value -= 4.0 * tpdm[pvqu] * tei[puqv];
400
value -= 2.0 * tpdm[pquv] * tei[pquv];
403
hess[pair] += 2.0 * value;
404
if (Params.scale_act_act != 1.0) hess[pair] *= Params.scale_act_act;
405
} /* end case H_{t1t2,t1t2} */
409
"(form_diag_mo_hess): Error, unrecognized class of indep pair\n");
412
} /* end loop over pairs */
419
** Calculates the full MO Hessian d^2E / (Theta_{pq} Theta_{rs})
420
** for independent pairs (p,q) and (r,s). Assume independent pairs
421
** are coming in such that p>=q.
423
** d^2 E / (Theta_{pq} Theta_{rs}) =
424
** y_{pqrs} - y_{qprs} - y_{pqsr} + y_{qpsr}
425
** + 1/2 delta_{ps} ( x_{qr} + x_{rq} )
426
** - 1/2 delta_{qs} ( x_{pr} + x_{rp} )
427
** - 1/2 delta_{pr} ( x_{qs} + x_{sq} )
428
** + 1/2 delta_{qr} ( x_{ps} + x_{sp} )
430
** x is the Lagrangian
431
** y_{pqrs} = gamma_{qs}
432
** + \sum_{mn} ( 2 Gamma_{qsmn} (pr|mn) + 4 Gamma_{qmsn} (pm|rn) )
434
** note: indices q and s must be populated for y to be nonzero
436
** Based on notes by Yukio Yamaguchi (MCSCF eq. 22)
437
** We need twice his value to be consistent with the "actual" gradient
438
** (his eq. 21, not his eq. 27) or to match G. Chaban, M. W. Schmidt,
439
** and M. S. Gordon, Theor. Chem. Acc. 97, 88 (1997).
444
void form_full_mo_hess(int npairs, int *ppair, int *qpair, double *oei,
445
double *tei, double **opdm, double *tpdm, double **lag,
449
int pair1, pair2, p, q, r, s, m, n;
450
int qs, ps, pr, qr, mn, qm, sn, pm, rn;
451
int qsmn, prmn, qmsn, pmrn, psmn, pmsn, qrmn, qmrn;
454
fprintf(outfile, "Forming full MCSCF orbital Hessian\n");
457
npop = CalcInfo.npop;
459
/* loop over the pairs of independent pairs */
460
for (pair1=0; pair1<npairs; pair1++) {
464
for (pair2=0; pair2<=pair1; pair2++) {
468
/* compute the hessian contribution for p,q,r,s */
470
if (p == s) sum += 0.5 * (lag[q][r] + lag[r][q]);
471
if (q == s) sum -= 0.5 * (lag[p][r] + lag[r][p]);
472
if (p == r) sum -= 0.5 * (lag[q][s] + lag[s][q]);
473
if (q == r) sum += 0.5 * (lag[p][s] + lag[s][p]);
480
/* compute y_{pqrs} */
482
if (q < npop && s < npop) {
483
ysum += oei[pr] * opdm[q][s];
484
for (m=0; m<npop; m++) {
487
for (n=0; n<npop; n++) {
495
ysum += tpdm[qsmn] * tei[prmn];
496
ysum += 2.0 * tpdm[qmsn] * tei[pmrn];
502
/* compute y_{qprs} */
504
if (p < npop && s < npop) {
505
ysum += oei[qr] * opdm[p][s];
506
for (m=0; m<npop; m++) {
509
for (n=0; n<npop; n++) {
517
ysum += tpdm[psmn] * tei[qrmn];
518
ysum += 2.0 * tpdm[pmsn] * tei[qmrn];
524
/* compute y_{pqsr} */
526
if (q < npop && r < npop) {
527
ysum += oei[ps] * opdm[q][r];
528
for (m=0; m<npop; m++) {
531
for (n=0; n<npop; n++) {
539
ysum += tpdm[qrmn] * tei[psmn];
540
ysum += 2.0 * tpdm[qmrn] * tei[pmsn];
546
/* compute y_{qpsr} */
548
if (p < npop && r < npop) {
549
ysum += oei[qs] * opdm[p][r];
550
for (m=0; m<npop; m++) {
553
for (n=0; n<npop; n++) {
561
ysum += tpdm[prmn] * tei[qsmn];
562
ysum += 2.0 * tpdm[pmrn] * tei[qmsn];
568
hess[pair1][pair2] = 2.0 * sum;
569
hess[pair2][pair1] = 2.0 * sum;
570
} /* end loop over pair2 */
571
} /* end loop over pair1 */
576
** form_diag_mo_hess_yy
578
** Calculates the diagonal MO Hessian d^2E / Theta_{pq}^2
579
** for independent pair (p,q). Assume independent pairs
580
** are coming in such that p>=q.
582
** d^2 E / Theta_{pq}^2
583
** y_{pqpq} - y_{qppq} - y_{pqqp} + y_{qpqp}
584
** - 1/2 delta_{qq} ( x_{pp} + x_{pp} )
585
** - 1/2 delta_{pp} ( x_{qq} + x_{qq} )
587
** x is the Lagrangian
588
** y_{pqrs} = gamma_{qs}
589
** + \sum_{mn} ( 2 Gamma_{qsmn} (pr|mn) + 4 Gamma_{qmsn} (pm|rn) )
591
** note: indices q and p must be populated for y to be nonzero
593
** Based on notes by Yukio Yamaguchi (MCSCF eq. 22)
598
void form_diag_mo_hess_yy(int npairs, int *ppair, int *qpair, double *oei,
599
double *tei, double **opdm, double *tpdm, double **lag,
603
int pair1, pair2, p, q, m, n;
604
int pq, pp, qq, mn, pn, qn, pm, qm;
605
int qqmn, ppmn, qmqn, pmpn, pqmn, qpmn, pmqn, qmpn;
608
fprintf(outfile, "Forming diagonal MCSCF orbital Hessian (YY)\n");
611
npop = CalcInfo.npop;
613
/* loop over the pairs of independent pairs */
614
for (pair1=0; pair1<npairs; pair1++) {
622
/* compute the hessian contribution for p,q,p,q */
623
sum = - lag[p][p] - lag[q][q];
625
/* compute y_{pqpq} */
628
ysum += oei[pp] * opdm[q][q];
630
for (m=0; m<npop; m++) {
633
for (n=0; n<npop; n++) {
641
ysum += tpdm[qqmn] * tei[ppmn];
642
ysum += 2.0 * tpdm[qmqn] * tei[pmpn];
648
/* compute y_{qppq} = y_{pqqp} */
650
if (p < npop && q < npop) {
651
ysum += oei[pq] * opdm[p][q];
652
for (m=0; m<npop; m++) {
655
for (n=0; n<npop; n++) {
663
ysum += tpdm[pqmn] * tei[pqmn];
664
ysum += 2.0 * tpdm[pmqn] * tei[qmpn];
670
/* compute y_{qpqp} */
673
ysum += oei[qq] * opdm[p][p];
674
for (m=0; m<npop; m++) {
677
for (n=0; n<npop; n++) {
686
ysum += tpdm[ppmn] * tei[qqmn];
687
ysum += 2.0 * tpdm[pmpn] * tei[qmqn];
693
hess[pair1] = 2.0 * sum;
694
} /* end loop over pair1 */
added
removed
src/bin/detcas/hessian.c
