8
/* may no longer need #include <libc.h> */
9
#include <libciomr/libciomr.h>
11
#include <libiwl/iwl.h>
21
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
24
void tpdm_block(struct stringwr **alplist, struct stringwr **betlist,
25
int nbf, int nalpcodes, int nbetcodes,
26
double *twopdm, double **CJ, double **CI, int Ja_list,
27
int Jb_list, int Jnas, int Jnbs, int Ia_list, int Ib_list,
31
void tpdm(struct stringwr **alplist, struct stringwr **betlist,
32
int Inroots, int Iroot, int Inunits, int Ifirstunit,
33
int Jnroots, int Jroot, int Jnunits, int Jfirstunit,
34
int targetfile, int writeflag, int printflag)
39
int i, j, k, l, lmax, ij, kl, ijkl, ijksym;
40
int i2, j2, k2, l2, nfzc, populated_orbs;
42
int maxrows, maxcols, ntri, ntri2;
44
double **transp_tmp = NULL;
45
double **transp_tmp2 = NULL;
46
double *buffer1, *buffer2, *twopdm, **onepdm, value;
47
int Iblock, Iblock2, Ibuf, Iac, Ibc, Inas, Inbs, Iairr;
48
int Jblock, Jblock2, Jbuf, Jac, Jbc, Jnas, Jnbs, Jairr;
49
int do_Jblock, do_Jblock2;
52
nfzc = CalcInfo.num_fzc_orbs;
53
populated_orbs = CalcInfo.nmo - CalcInfo.num_fzv_orbs;
56
psio_open(Parameters.opdm_file, PSIO_OPEN_OLD);
57
onepdm = block_matrix(populated_orbs, populated_orbs);
60
Ivec.set(CIblks.vectlen, CIblks.num_blocks, Parameters.icore, 1,
61
CIblks.Ia_code, CIblks.Ib_code, CIblks.Ia_size, CIblks.Ib_size,
62
CIblks.offset, CIblks.num_alp_codes, CIblks.num_bet_codes,
63
CalcInfo.nirreps, AlphaG->subgr_per_irrep, Inroots, Inunits,
64
Ifirstunit, CIblks.first_iablk, CIblks.last_iablk, CIblks.decode);
66
Jvec.set(CIblks.vectlen, CIblks.num_blocks, Parameters.icore, 1,
67
CIblks.Ia_code, CIblks.Ib_code, CIblks.Ia_size, CIblks.Ib_size,
68
CIblks.offset, CIblks.num_alp_codes, CIblks.num_bet_codes,
69
CalcInfo.nirreps, AlphaG->subgr_per_irrep, Jnroots, Jnunits,
70
Jfirstunit, CIblks.first_iablk, CIblks.last_iablk, CIblks.decode);
72
buffer1 = Ivec.buf_malloc();
73
buffer2 = Jvec.buf_malloc();
74
Ivec.buf_lock(buffer1);
75
Jvec.buf_lock(buffer2);
77
ntri = CalcInfo.num_ci_orbs * CalcInfo.num_ci_orbs;
78
ntri2 = (ntri * (ntri + 1)) / 2;
79
twopdm = init_array(ntri2);
81
if ((Ivec.icore==2 && Ivec.Ms0 && CalcInfo.ref_sym != 0) ||
82
(Ivec.icore==0 && Ivec.Ms0)) {
83
for (i=0, maxrows=0, maxcols=0; i<Ivec.num_blocks; i++) {
84
if (Ivec.Ia_size[i] > maxrows) maxrows = Ivec.Ia_size[i];
85
if (Ivec.Ib_size[i] > maxcols) maxcols = Ivec.Ib_size[i];
87
if (maxcols > maxrows) maxrows = maxcols;
88
transp_tmp = (double **) malloc (maxrows * sizeof(double *));
89
transp_tmp2 = (double **) malloc (maxrows * sizeof(double *));
90
if (transp_tmp == NULL || transp_tmp2 == NULL) {
91
printf("(tpdm): Trouble with malloc'ing transp_tmp\n");
93
bufsz = Ivec.get_max_blk_size();
94
transp_tmp[0] = init_array(bufsz);
95
transp_tmp2[0] = init_array(bufsz);
96
if (transp_tmp[0] == NULL || transp_tmp2[0] == NULL) {
97
printf("(tpdm): Trouble with malloc'ing transp_tmp[0]\n");
102
if (Parameters.icore == 0) {
104
for (Ibuf=0; Ibuf<Ivec.buf_per_vect; Ibuf++) {
105
Ivec.read(Iroot, Ibuf);
106
Iblock = Ivec.buf2blk[Ibuf];
107
Iac = Ivec.Ia_code[Iblock];
108
Ibc = Ivec.Ib_code[Iblock];
109
Inas = Ivec.Ia_size[Iblock];
110
Inbs = Ivec.Ib_size[Iblock];
112
for (Jbuf=0; Jbuf<Jvec.buf_per_vect; Jbuf++) {
113
do_Jblock=0; do_Jblock2=0;
114
Jblock = Jvec.buf2blk[Jbuf];
116
Jac = Jvec.Ia_code[Jblock];
117
Jbc = Jvec.Ib_code[Jblock];
118
if (Jvec.Ms0) Jblock2 = Jvec.decode[Jbc][Jac];
119
Jnas = Jvec.Ia_size[Jblock];
120
Jnbs = Jvec.Ib_size[Jblock];
121
if (s1_contrib[Iblock][Jblock] || s2_contrib[Iblock][Jblock]
122
|| s3_contrib[Iblock][Jblock])
124
if (Jvec.buf_offdiag[Jbuf] && (s1_contrib[Iblock][Jblock2] ||
125
s2_contrib[Iblock][Jblock2] ||
126
s3_contrib[Iblock][Jblock2]))
128
if (!do_Jblock && !do_Jblock2) continue;
130
Jvec.read(Jroot, Jbuf);
133
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
134
Ivec.num_alpcodes, Ivec.num_betcodes, twopdm,
135
Jvec.blocks[Jblock], Ivec.blocks[Iblock],
136
Jac, Jbc, Jnas, Jnbs, Iac, Ibc, Inas, Inbs);
140
Jvec.transp_block(Jblock, transp_tmp);
141
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
142
Ivec.num_alpcodes, Ivec.num_betcodes,
143
twopdm, transp_tmp, Ivec.blocks[Iblock],
144
Jbc, Jac, Jnbs, Jnas, Iac, Ibc, Inas, Inbs);
147
} /* end loop over Jbuf */
149
if (Ivec.buf_offdiag[Ibuf]) { /* need to get contrib of transpose */
150
Iblock2 = Ivec.decode[Ibc][Iac];
151
Iac = Ivec.Ia_code[Iblock2];
152
Ibc = Ivec.Ib_code[Iblock2];
153
Inas = Ivec.Ia_size[Iblock2];
154
Inbs = Ivec.Ib_size[Iblock2];
156
Ivec.transp_block(Iblock, transp_tmp2);
158
for (Jbuf=0; Jbuf<Jvec.buf_per_vect; Jbuf++) {
159
do_Jblock=0; do_Jblock2=0;
160
Jblock = Jvec.buf2blk[Jbuf];
162
Jac = Jvec.Ia_code[Jblock];
163
Jbc = Jvec.Ib_code[Jblock];
164
if (Jvec.Ms0) Jblock2 = Jvec.decode[Jbc][Jac];
165
Jnas = Jvec.Ia_size[Jblock];
166
Jnbs = Jvec.Ib_size[Jblock];
167
if (s1_contrib[Iblock2][Jblock] || s2_contrib[Iblock2][Jblock] ||
168
s3_contrib[Iblock2][Jblock])
170
if (Jvec.buf_offdiag[Jbuf] && (s1_contrib[Iblock2][Jblock2] ||
171
s2_contrib[Iblock2][Jblock2] ||
172
s3_contrib[Iblock2][Jblock2]))
174
if (!do_Jblock && !do_Jblock2) continue;
176
Jvec.read(Jroot, Jbuf);
179
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
180
Ivec.num_alpcodes, Ivec.num_betcodes,
181
twopdm, Jvec.blocks[Jblock],
182
transp_tmp2, Jac, Jbc, Jnas,
183
Jnbs, Iac, Ibc, Inas, Inbs);
187
Jvec.transp_block(Jblock, transp_tmp);
188
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
189
Ivec.num_alpcodes, Ivec.num_betcodes,
190
twopdm, transp_tmp, transp_tmp2,
191
Jbc, Jac, Jnbs, Jnas, Iac, Ibc, Inas, Inbs);
193
} /* end loop over Jbuf */
195
} /* end loop over Ibuf transpose */
196
} /* end loop over Ibuf */
199
else if (Parameters.icore==1) { /* whole vectors in-core */
202
for (Iblock=0; Iblock<Ivec.num_blocks; Iblock++) {
203
Iac = Ivec.Ia_code[Iblock];
204
Ibc = Ivec.Ib_code[Iblock];
205
Inas = Ivec.Ia_size[Iblock];
206
Inbs = Ivec.Ib_size[Iblock];
207
if (Inas==0 || Inbs==0) continue;
208
for (Jblock=0; Jblock<Jvec.num_blocks; Jblock++) {
209
Jac = Jvec.Ia_code[Jblock];
210
Jbc = Jvec.Ib_code[Jblock];
211
Jnas = Jvec.Ia_size[Jblock];
212
Jnbs = Jvec.Ib_size[Jblock];
213
if (s1_contrib[Iblock][Jblock] || s2_contrib[Iblock][Jblock] ||
214
s3_contrib[Iblock][Jblock])
215
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
216
Ivec.num_alpcodes, Ivec.num_betcodes,
217
twopdm, Jvec.blocks[Jblock], Ivec.blocks[Iblock],
218
Jac, Jbc, Jnas, Jnbs, Iac, Ibc, Inas, Inbs);
221
} /* end loop over Iblock */
224
else if (Parameters.icore==2) { /* icore==2 */
225
for (Ibuf=0; Ibuf<Ivec.buf_per_vect; Ibuf++) {
226
Ivec.read(Iroot, Ibuf);
227
Iairr = Ivec.buf2blk[Ibuf];
229
for (Jbuf=0; Jbuf<Jvec.buf_per_vect; Jbuf++) {
230
Jvec.read(Jroot, Jbuf);
231
Jairr = Jvec.buf2blk[Jbuf];
233
for (Iblock=Ivec.first_ablk[Iairr]; Iblock<=Ivec.last_ablk[Iairr];
235
Iac = Ivec.Ia_code[Iblock];
236
Ibc = Ivec.Ib_code[Iblock];
237
Inas = Ivec.Ia_size[Iblock];
238
Inbs = Ivec.Ib_size[Iblock];
240
for (Jblock=Jvec.first_ablk[Jairr]; Jblock<=Jvec.last_ablk[Jairr];
242
Jac = Jvec.Ia_code[Jblock];
243
Jbc = Jvec.Ib_code[Jblock];
244
Jnas = Jvec.Ia_size[Jblock];
245
Jnbs = Jvec.Ib_size[Jblock];
247
if (s1_contrib[Iblock][Jblock] || s2_contrib[Iblock][Jblock] ||
248
s3_contrib[Iblock][Jblock])
249
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
250
Ivec.num_alpcodes, Ivec.num_betcodes,
251
twopdm, Jvec.blocks[Jblock], Ivec.blocks[Iblock],
252
Jac, Jbc, Jnas, Jnbs, Iac, Ibc, Inas, Inbs);
254
if (Jvec.buf_offdiag[Jbuf]) {
255
Jblock2 = Jvec.decode[Jbc][Jac];
256
if (s1_contrib[Iblock][Jblock2] ||
257
s2_contrib[Iblock][Jblock2] ||
258
s3_contrib[Iblock][Jblock2]) {
259
Jvec.transp_block(Jblock, transp_tmp);
260
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
261
Ivec.num_alpcodes, Ivec.num_betcodes,
262
twopdm, transp_tmp, Ivec.blocks[Iblock],
263
Jbc, Jac, Jnbs, Jnas, Iac, Ibc, Inas, Inbs);
267
} /* end loop over Jblock */
269
if (Ivec.buf_offdiag[Ibuf]) {
270
Iblock2 = Ivec.decode[Ibc][Iac];
271
Ivec.transp_block(Iblock, transp_tmp2);
272
Iac = Ivec.Ia_code[Iblock2];
273
Ibc = Ivec.Ib_code[Iblock2];
274
Inas = Ivec.Ia_size[Iblock2];
275
Inbs = Ivec.Ib_size[Iblock2];
277
for (Jblock=Jvec.first_ablk[Jairr]; Jblock<=Jvec.last_ablk[Jairr];
279
Jac = Jvec.Ia_code[Jblock];
280
Jbc = Jvec.Ib_code[Jblock];
281
Jnas = Jvec.Ia_size[Jblock];
282
Jnbs = Jvec.Ib_size[Jblock];
284
if (s1_contrib[Iblock2][Jblock] || s2_contrib[Iblock2][Jblock]
285
|| s3_contrib[Iblock2][Jblock])
286
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
287
Ivec.num_alpcodes, Ivec.num_betcodes,
288
twopdm, Jvec.blocks[Jblock], transp_tmp2,
289
Jac, Jbc, Jnas, Jnbs, Iac, Ibc, Inas, Inbs);
291
if (Jvec.buf_offdiag[Jbuf]) {
292
Jblock2 = Jvec.decode[Jbc][Jac];
293
if (s1_contrib[Iblock][Jblock2] ||
294
s2_contrib[Iblock][Jblock2] ||
295
s3_contrib[Iblock][Jblock2]) {
296
Jvec.transp_block(Jblock, transp_tmp);
297
tpdm_block(alplist, betlist, CalcInfo.num_ci_orbs,
298
Ivec.num_alpcodes, Ivec.num_betcodes,
299
twopdm, transp_tmp, transp_tmp2, Jbc, Jac,
300
Jnbs, Jnas, Iac, Ibc, Inas, Inbs);
304
} /* end loop over Jblock */
305
} /* end Ivec offdiag */
307
} /* end loop over Iblock */
308
} /* end loop over Jbuf */
309
} /* end loop over Ibuf */
313
printf("tpdm: unrecognized core option!\n");
317
/* write and/or print the tpdm */
321
if (printflag) fprintf(outfile, "\nTwo-particle density matrix\n\n");
322
iwl_buf_init(&TBuff, targetfile, 0.0, 0, 0);
323
orbsym = CalcInfo.orbsym + nfzc;
325
/* do the core-core and core-active part here */
327
sprintf(opdm_key, "MO-basis OPDM Root %d", Iroot);
328
psio_read_entry(Parameters.opdm_file, opdm_key, (char *) onepdm[0],
329
populated_orbs * populated_orbs * sizeof(double));
332
for (i=0; i<nfzc; i++) {
333
for (j=0; j<i; j++) {
334
iwl_buf_wrt_val(&TBuff,i,i,j,j, 2.00,printflag,outfile,0);
335
iwl_buf_wrt_val(&TBuff,i,j,j,i,-1.00,printflag,outfile,0);
337
iwl_buf_wrt_val(&TBuff,i,i,i,i,1.0,printflag,outfile,0);
340
/* core-active part */
341
for (i=nfzc; i<populated_orbs; i++) {
342
for (j=nfzc; j<populated_orbs; j++) {
343
value = onepdm[i][j];
344
for (k=0; k<nfzc; k++) {
345
iwl_buf_wrt_val(&TBuff,i,j,k,k,value,printflag,outfile,0);
346
iwl_buf_wrt_val(&TBuff,i,k,k,j,-0.5*value,printflag,outfile,0);
352
for (i=0; i<CalcInfo.num_ci_orbs; i++) {
354
for (j=0; j<CalcInfo.num_ci_orbs; j++) {
356
for (k=0; k<=i; k++) {
358
if (k==i) lmax = j+1;
359
else lmax = CalcInfo.num_ci_orbs;
360
ijksym = orbsym[i] ^ orbsym[j] ^ orbsym[k];
361
for (l=0; l<lmax; l++) {
363
if ((orbsym[l] ^ ijksym) != 0) continue;
364
ij = i * CalcInfo.num_ci_orbs + j;
365
kl = k * CalcInfo.num_ci_orbs + l;
368
// For PSI the factor of 1/2 is not pulled outside...so put
369
// it back inside now and then write out to the IWL buffer.
370
value = 0.5 * twopdm[ijkl];
371
iwl_buf_wrt_val(&TBuff,i2,j2,k2,l2,value,printflag,outfile,0);
376
iwl_buf_flush(&TBuff, 1);
377
iwl_buf_close(&TBuff, 1);
378
fprintf(outfile, "\n");
383
psio_close(Parameters.opdm_file, 1);
390
if (transp_tmp != NULL) free_block(transp_tmp);
391
if (transp_tmp2 != NULL) free_block(transp_tmp2);
398
void tpdm_block(struct stringwr **alplist, struct stringwr **betlist,
399
int nbf, int nalplists, int nbetlists,
400
double *twopdm, double **CJ, double **CI, int Ja_list,
401
int Jb_list, int Jnas, int Jnbs, int Ia_list, int Ib_list,
404
int Ia_idx, Ib_idx, Ja_idx, Jb_idx, Ja_ex, Jb_ex, Jbcnt, Jacnt;
405
int Kbcnt, Kacnt, Kb_ex, Ka_ex, Kb_list, Ka_list, Kb_idx, Ka_idx;
406
struct stringwr *Jb, *Ja, *Kb, *Ka;
407
signed char *Jbsgn, *Jasgn, *Kbsgn, *Kasgn;
408
unsigned int *Jbridx, *Jaridx, *Kbridx, *Karidx;
409
double C1, C2, Ib_sgn, Ia_sgn, Kb_sgn, Ka_sgn, tval;
410
int i, j, k, l, ij, kl, ijkl, oij, okl, *Jboij, *Jaoij, *Kboij, *Kaoij;
412
/* loop over Ia in Ia_list */
413
if (Ia_list == Ja_list) {
414
for (Ia_idx=0; Ia_idx<Inas; Ia_idx++) {
415
for (Jb=betlist[Jb_list], Jb_idx=0; Jb_idx<Jnbs; Jb_idx++, Jb++) {
416
C1 = CJ[Ia_idx][Jb_idx];
418
/* loop over excitations E^b_{kl} from |B(J_b)> */
419
for (Kb_list=0; Kb_list < nbetlists; Kb_list++) {
420
Jbcnt = Jb->cnt[Kb_list];
421
Jbridx = Jb->ridx[Kb_list];
422
Jbsgn = Jb->sgn[Kb_list];
423
Jboij = Jb->oij[Kb_list];
424
for (Jb_ex=0; Jb_ex < Jbcnt; Jb_ex++) {
427
Kb_sgn = (double) *Jbsgn++;
429
Kb = betlist[Kb_list] + Kb_idx;
430
if (Kb_list == Ib_list) {
431
C2 = CI[Ia_idx][Kb_idx];
434
for (j=0; j<nbf && j<=i; j++) {
439
twopdm[ijkl] -= Kb_sgn * C1 * C2;
444
/* loop over excitations E^b_{ij} from |B(K_b)> */
445
/* Ib_list pre-determined because of C blocking */
446
Kbcnt = Kb->cnt[Ib_list];
447
Kbridx = Kb->ridx[Ib_list];
448
Kbsgn = Kb->sgn[Ib_list];
449
Kboij = Kb->oij[Ib_list];
450
for (Kb_ex=0; Kb_ex<Kbcnt; Kb_ex++) {
452
Ib_sgn = (double) *Kbsgn++;
455
C2 = CI[Ia_idx][Ib_idx];
456
ijkl = INDEX(oij,okl);
457
twopdm[ijkl] += Ib_sgn * Kb_sgn * C1 * C2;
461
} /* end loop over Jb_ex */
462
} /* end loop over Kb_list */
463
} /* end loop over Jb_idx */
464
} /* end loop over Ia_idx */
465
} /* end case Ia_list == Ja_list */
467
/* loop over Ib in Ib_list */
468
if (Ib_list == Jb_list) {
469
for (Ib_idx=0; Ib_idx<Inbs; Ib_idx++) {
470
for (Ja=alplist[Ja_list], Ja_idx=0; Ja_idx<Jnas; Ja_idx++, Ja++) {
471
C1 = CJ[Ja_idx][Ib_idx];
473
/* loop over excitations E^a_{kl} from |A(J_a)> */
474
for (Ka_list=0; Ka_list < nalplists; Ka_list++) {
475
Jacnt = Ja->cnt[Ka_list];
476
Jaridx = Ja->ridx[Ka_list];
477
Jasgn = Ja->sgn[Ka_list];
478
Jaoij = Ja->oij[Ka_list];
479
for (Ja_ex=0; Ja_ex < Jacnt; Ja_ex++) {
482
Ka_sgn = (double) *Jasgn++;
484
Ka = alplist[Ka_list] + Ka_idx;
485
if (Ka_list == Ia_list) {
486
C2 = CI[Ka_idx][Ib_idx];
489
for (j=0; j<nbf && j<=i; j++) {
494
twopdm[ijkl] -= Ka_sgn * C1 * C2;
499
/* loop over excitations E^a_{ij} from |A(K_a)> */
500
/* Ia_list pre-determined because of C blocking */
501
Kacnt = Ka->cnt[Ia_list];
502
Karidx = Ka->ridx[Ia_list];
503
Kasgn = Ka->sgn[Ia_list];
504
Kaoij = Ka->oij[Ia_list];
505
for (Ka_ex=0; Ka_ex<Kacnt; Ka_ex++) {
507
Ia_sgn = (double) *Kasgn++;
510
C2 = CI[Ia_idx][Ib_idx];
511
ijkl = INDEX(oij,okl);
512
twopdm[ijkl] += Ia_sgn * Ka_sgn * C1 * C2;
516
} /* end loop over Ja_ex */
517
} /* end loop over Ka_list */
518
} /* end loop over Ja_idx */
519
} /* end loop over Ib_idx */
520
} /* end case Ib_list == Jb_list */
523
/* now do the sigma3 looking (alpha-beta) part */
525
for (Ja=alplist[Ja_list], Ja_idx=0; Ja_idx<Jnas; Ja_idx++, Ja++) {
527
/* loop over excitations E^a_{kl} from |A(I_a)> */
528
Jacnt = Ja->cnt[Ia_list];
529
Jaridx = Ja->ridx[Ia_list];
530
Jasgn = Ja->sgn[Ia_list];
531
Jaoij = Ja->oij[Ia_list];
532
for (Ja_ex=0; Ja_ex < Jacnt; Ja_ex++) {
535
Ia_sgn = (double) *Jasgn++;
538
for (Jb=betlist[Jb_list], Jb_idx=0; Jb_idx<Jnbs; Jb_idx++, Jb++) {
540
C1 = CJ[Ja_idx][Jb_idx];
542
/* loop over excitations E^b_{ij} from |B(J_b)> */
543
Jbcnt = Jb->cnt[Ib_list];
544
Jbridx = Jb->ridx[Ib_list];
545
Jbsgn = Jb->sgn[Ib_list];
546
Jboij = Jb->oij[Ib_list];
548
for (Jb_ex=0; Jb_ex < Jbcnt; Jb_ex++) {
551
Ib_sgn = (double) *Jbsgn++;
552
C2 = CI[Ia_idx][Ib_idx];
553
ijkl = INDEX(oij, okl);
554
tval = Ib_sgn * Ia_sgn * C1 * C2;
555
if (oij == okl) tval *= 2.0;
556
twopdm[ijkl] += tval;
558
} /* end loop over Jb */
559
} /* end loop over Ja_ex */
560
} /* end loop over Ja */
added
removed
src/bin/detci/tpdm.cc
