3
#include <libipv1/ip_lib.h>
4
#include <libciomr/libciomr.h>
5
#include <libpsio/psio.h>
7
#include <libiwl/iwl.h>
8
#include <libchkpt/chkpt.h>
15
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
16
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
17
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
19
void trans_one_forwards(void);
20
void trans_one_backwards(void);
21
void tran_one(int nirreps, double ***C, int src_orbs, int *src_first, int *src_last,
22
int *src_orbspi, int dst_orbs, int *dst_first, int *dst_last,
23
int *dst_orbspi, double *src_ints, double *dst_ints, int *order,
24
char *label, int backtran, int nfzc, int printflg,FILE *outfile);
25
double *** construct_evects(char *spin, int nirreps, int *active, int *sopi, int *orbspi,
26
int *first_so, int *last_so, int *first, int *last,
27
int *fstact, int *lstact, int printflag);
30
** TRANSFORM_ONE(): This function does the SO to MO transformation of
31
** the one-electron integrals.
33
** This function was getting too complicated to read, so I split
34
** it up into two functions, one for forwards transforms, and one
35
** for backwards transforms.
44
trans_one_backwards();
55
** TRANS_ONE_FORWARDS()
57
** This does a forwards transform of one-electron integrals from the SO
58
** to the MO basis. All the real work is done in a function
59
** tran_one(); the rest is just setting up to call that function.
60
** I also simplified matters by changing libiwl to allow filtering
61
** out of frozen orbitals. Hence, for forwards transforms, we now
62
** only write (up to) two files: the bare one-electron integrals,
63
** and the frozen core operator. Each file runs over all indices
64
** (frozen or not). This makes this code much simpler than it was.
69
void trans_one_forwards(void)
71
int itap, nirreps, nmo, nso, nfzc, nfzv;
72
int *src_first, *dst_first, *src_last, *dst_last, *src_orbspi, *dst_orbspi;
73
int src_orbs, dst_orbs, src_ntri, dst_ntri;
75
double *oe_ints, *new_oe_ints;
77
if (params.print_lvl) {
78
fprintf(outfile, "\n\tTransforming one-electron integrals...\n");
82
nirreps = moinfo.nirreps;
85
print_integrals = params.print_oe_ints;
88
src_first = moinfo.first_so;
89
dst_first = moinfo.first;
90
src_last = moinfo.last_so;
91
dst_last = moinfo.last;
92
src_orbspi = moinfo.sopi;
93
dst_orbspi = moinfo.orbspi;
96
dst_orbs = nmo - nfzv - nfzc;
97
dst_ntri = (dst_orbs * (dst_orbs + 1)) / 2;
98
src_ntri = (src_orbs * (src_orbs + 1)) / 2;
100
new_oe_ints = init_array(dst_ntri);
102
/* We do the two-electron ints first, for which our target ints may
103
* not have allowed frozen indices. However, here we now want to run
104
* over all MOs. So, we need to re-build the C matrix to include
105
* frozen orbital columns. This annoying stuff could be avoided by
106
* a generalized mmult routine
108
if (!params.do_all_tei) {
109
if(!strcmp(params.ref,"UHF")) {
110
destruct_evects(nirreps, moinfo.evects_alpha);
111
destruct_evects(nirreps, moinfo.evects_beta);
113
else destruct_evects(nirreps, moinfo.evects);
115
if(!strcmp(params.ref,"UHF")) {
116
moinfo.evects_alpha = construct_evects("alpha", moinfo.nirreps,
118
moinfo.sopi, moinfo.orbspi,
119
moinfo.first_so, moinfo.last_so,
120
moinfo.first, moinfo.last,
121
moinfo.first, moinfo.last,
123
moinfo.evects_beta = construct_evects("beta", moinfo.nirreps,
125
moinfo.sopi, moinfo.orbspi,
126
moinfo.first_so, moinfo.last_so,
127
moinfo.first, moinfo.last,
128
moinfo.first, moinfo.last,
132
moinfo.evects = construct_evects("RHF", moinfo.nirreps, moinfo.orbspi,
133
moinfo.sopi, moinfo.orbspi,
134
moinfo.first_so, moinfo.last_so,
135
moinfo.first, moinfo.last,
136
moinfo.first, moinfo.last,
141
if (params.do_h_bare) {
143
oe_ints = moinfo.oe_ints;
145
if(!strcmp(params.ref,"UHF")) {
147
itap = params.h_bare_a_file;
149
tran_one(nirreps, moinfo.evects_alpha, src_orbs, src_first, src_last,
150
src_orbspi, dst_orbs, dst_first, dst_last, dst_orbspi, oe_ints,
151
new_oe_ints, moinfo.order_alpha,
152
"\n\n\tOne-Electron Alpha Integrals (MO Basis):\n",
153
params.backtr, nfzc, print_integrals, outfile);
155
iwl_wrtone(itap, PSIF_MO_A_OEI, dst_ntri, new_oe_ints);
157
if (params.print_lvl) {
158
fprintf(outfile, "\tOne-electron A integrals written to file %d.\n",
163
itap = params.h_bare_b_file;
165
/* Clear the new_oe_ints */
166
zero_arr(new_oe_ints,dst_ntri);
168
tran_one(nirreps, moinfo.evects_beta, src_orbs, src_first, src_last,
169
src_orbspi, dst_orbs, dst_first, dst_last, dst_orbspi, oe_ints,
170
new_oe_ints, moinfo.order_beta,
171
"\n\n\tOne-Electron Beta Integrals (MO Basis):\n",
172
params.backtr, nfzc, print_integrals, outfile);
174
iwl_wrtone(itap, PSIF_MO_B_OEI, dst_ntri, new_oe_ints);
175
if (params.print_lvl) {
176
fprintf(outfile, "\tOne-electron B integrals written to file %d.\n",
183
itap = params.h_bare_file;
185
tran_one(nirreps, moinfo.evects, src_orbs, src_first, src_last,
186
src_orbspi, dst_orbs, dst_first, dst_last, dst_orbspi, oe_ints,
187
new_oe_ints, moinfo.order,
188
"\n\n\tOne-Electron Integrals (MO Basis):\n",
189
params.backtr, nfzc, print_integrals, outfile);
192
iwl_wrtone(itap, PSIF_MO_OEI, dst_ntri, new_oe_ints);
193
if (params.print_lvl) {
194
fprintf(outfile, "\tOne-electron integrals written to file %d.\n",itap);
203
if (params.do_h_fzc) {
205
if(!strcmp(params.ref,"UHF")) {
207
itap = params.h_fzc_a_file;
208
oe_ints = moinfo.fzc_operator_alpha;
210
tran_one(nirreps, moinfo.evects_alpha, src_orbs, src_first, src_last,
211
src_orbspi, dst_orbs, dst_first, dst_last, dst_orbspi, oe_ints,
212
new_oe_ints, moinfo.order_alpha,
213
"\n\n\tAlpha Frozen-Core Operator (MO Basis):\n",
214
params.backtr, nfzc, print_integrals, outfile);
216
iwl_wrtone(itap, PSIF_MO_A_FZC, dst_ntri, new_oe_ints);
218
if (params.print_lvl) {
219
fprintf(outfile, "\tAlpha frozen-core operator written to file %d.\n",
224
/* Clear the new_oe_ints */
225
zero_arr(new_oe_ints,dst_ntri);
227
itap = params.h_fzc_b_file;
228
oe_ints = moinfo.fzc_operator_beta;
230
tran_one(nirreps, moinfo.evects_beta, src_orbs, src_first, src_last,
231
src_orbspi, dst_orbs, dst_first, dst_last, dst_orbspi, oe_ints,
232
new_oe_ints, moinfo.order_beta,
233
"\n\n\tBeta Frozen-Core Operator (MO Basis):\n",
234
params.backtr, nfzc, print_integrals, outfile);
236
iwl_wrtone(itap, PSIF_MO_B_FZC, dst_ntri, new_oe_ints);
238
if (params.print_lvl) {
239
fprintf(outfile, "\tBeta frozen-core operator written to file %d.\n",
246
itap = params.h_fzc_file;
247
oe_ints = moinfo.fzc_operator;
249
tran_one(nirreps, moinfo.evects, src_orbs, src_first, src_last,
250
src_orbspi, dst_orbs, dst_first, dst_last, dst_orbspi, oe_ints,
251
new_oe_ints, moinfo.order,
252
"\n\n\tFrozen-Core Operator (MO Basis):\n",
253
params.backtr, nfzc, print_integrals, outfile);
255
iwl_wrtone(itap, PSIF_MO_FZC, dst_ntri, new_oe_ints);
257
if (params.print_lvl) {
258
fprintf(outfile, "\tFrozen-core operator written to file %d.\n", itap);
272
** TRANS_ONE_BACKWARDS()
274
** This does a backwards transform of the one-particle density matrix
275
** and the MO Lagrangian from the MO to the AO basis. All the real
276
** work is done in a function tran_one()---the same one called in the
277
** forwards transformations; the code in this function just sets up
278
** to call tran_one().
283
void trans_one_backwards(void)
285
int itap, nirreps, nmo, nfzc, nfzv;
286
int *src_first, *dst_first, *src_last, *dst_last, *src_orbspi, *dst_orbspi;
287
int src_orbs, dst_orbs, src_ntri, dst_ntri;
289
double *oe_ints, *new_oe_ints;
290
double **A,**B,*opdm,**tmat,**tmat2,**so2ao,tval;
291
double *opdm_a, *opdm_b, *lag_a, *lag_b;
292
double *new_opdm_a, *new_opdm_b, *new_lag_a, *new_lag_b;
295
PSI_FPTR pdm_fp_in=0, pdm_fp_out=0;
297
if (params.print_lvl) {
298
fprintf(outfile, "\n\tTransforming one-electron integrals...\n");
302
itap = params.opdm_out_file;
303
nirreps = moinfo.backtr_nirreps;
305
print_integrals = params.print_oe_ints;
309
src_first = moinfo.backtr_mo_first;
310
dst_first = moinfo.backtr_ao_first;
311
src_last = moinfo.backtr_mo_lstact;
312
dst_last = moinfo.backtr_ao_last;
313
src_orbspi = moinfo.backtr_mo_active;
314
dst_orbspi = moinfo.backtr_ao_orbspi;
316
src_orbs = nmo - nfzv;
317
dst_orbs = moinfo.nao;
318
dst_ntri = (dst_orbs * (dst_orbs + 1)) / 2;
319
src_ntri = (src_orbs * (src_orbs + 1)) / 2;
321
if(!strcmp(params.ref, "UHF")) {
323
tmat = block_matrix(src_orbs, src_orbs);
324
psio_open(PSIF_MO_OPDM, PSIO_OPEN_OLD);
326
psio_read_entry(PSIF_MO_OPDM, "MO-basis Alpha OPDM", (char *) tmat[0],
327
sizeof(double)*src_orbs*src_orbs);
329
opdm_a = init_array(src_ntri);
330
for (p=0; p<src_orbs; p++) {
331
for (q=0; q<=p; q++) {
332
P = moinfo.corr2pitz_nofzv_a[p];
333
Q = moinfo.corr2pitz_nofzv_a[q];
335
opdm_a[PQ] = 0.5 * (tmat[p][q] + tmat[q][p]);
339
new_opdm_a = init_array(dst_ntri);
340
tran_one(nirreps, moinfo.evects_alpha, src_orbs, src_first, src_last, src_orbspi, dst_orbs,
341
dst_first, dst_last, dst_orbspi, opdm_a, new_opdm_a, moinfo.order,
342
"\n\n\tAlpha Contribution to One-Particle Density Matrix (AO Basis):\n",
343
params.backtr, nfzc, print_integrals, outfile);
346
psio_read_entry(PSIF_MO_OPDM, "MO-basis Beta OPDM", (char *) tmat[0],
347
sizeof(double)*src_orbs*src_orbs);
349
opdm_b = init_array(src_ntri);
350
for (p=0; p<src_orbs; p++) {
351
for (q=0; q<=p; q++) {
352
P = moinfo.corr2pitz_nofzv_b[p];
353
Q = moinfo.corr2pitz_nofzv_b[q];
355
opdm_b[PQ] = 0.5 * (tmat[p][q] + tmat[q][p]);
359
new_opdm_b = init_array(dst_ntri);
360
tran_one(nirreps, moinfo.evects_beta, src_orbs, src_first, src_last, src_orbspi, dst_orbs,
361
dst_first, dst_last, dst_orbspi, opdm_b, new_opdm_b, moinfo.order,
362
"\n\n\tBeta Contribution to One-Particle Density Matrix (AO Basis):\n",
363
params.backtr, nfzc, print_integrals, outfile);
366
psio_close(PSIF_MO_OPDM, 1);
369
if(print_integrals) {
370
fprintf(outfile, "\n\tAlpha AO-basis OPDM:\n");
371
print_array(new_opdm_a, dst_orbs, outfile);
372
fprintf(outfile, "\n\tBeta AO-basis OPDM:\n");
373
print_array(new_opdm_b, dst_orbs, outfile);
376
/* combine the alpha and beta contributions */
377
for(p=0; p < dst_ntri; p++)
378
new_opdm_a[p] += new_opdm_b[p];
380
if(print_integrals) {
381
fprintf(outfile, "\n\tTotal AO-basis OPDM:\n");
382
print_array(new_opdm_a, dst_orbs, outfile);
385
psio_open(itap, PSIO_OPEN_OLD);
386
psio_write_entry(itap, "AO-basis OPDM", (char *) new_opdm_a, dst_ntri*sizeof(double));
394
/* read in the MO one-particle density matrix */
395
opdm = init_array(src_ntri);
397
tmat = block_matrix(src_orbs, src_orbs);
399
psio_open(params.opdm_in_file, PSIO_OPEN_OLD);
400
if (strcmp(params.wfn, "OOCCD")!=0)
401
psio_read_entry(params.opdm_in_file, "MO-basis OPDM", (char *) tmat[0],
402
src_orbs*src_orbs*sizeof(double));
404
tmat2 = block_matrix(nmo, nmo);
405
psio_read_entry(params.opdm_in_file, "MO-basis OPDM", (char *) tmat2[0],
406
nmo*nmo*sizeof(double));
407
for (p=0; p<src_orbs; p++)
408
for (q=0; q<src_orbs; q++)
409
tmat[p][q] = tmat2[p][q];
412
psio_close(params.opdm_in_file, 1);
413
if (params.print_lvl > 3) {
414
fprintf(outfile, "One-particle density matrix\n");
415
print_mat(tmat, src_orbs, src_orbs, outfile);
418
/* reorder the opdm into the backtransform order */
419
for (p=0; p<src_orbs; p++) {
420
for (q=0; q<=p; q++) {
421
P = moinfo.corr2pitz_nofzv[p];
422
Q = moinfo.corr2pitz_nofzv[q];
424
opdm[PQ] = 0.5 * (tmat[p][q] + tmat[q][p]);
429
new_oe_ints = init_array(dst_ntri);
431
tran_one(nirreps, moinfo.evects, src_orbs, src_first, src_last, src_orbspi, dst_orbs,
432
dst_first, dst_last, dst_orbspi, oe_ints, new_oe_ints, moinfo.order,
433
"\n\n\tOne-Particle Density Matrix (AO Basis):\n",
434
params.backtr, nfzc, print_integrals, outfile);
436
psio_open(itap, PSIO_OPEN_OLD);
437
psio_write_entry(itap, "AO-basis OPDM", (char *) new_oe_ints, dst_ntri*sizeof(double));
445
/* If this is a back-transform, we need to do the Lagrangian also */
446
/* Assume the lagrangian has indices over all orbitals, incl fzv */
447
src_last = moinfo.backtr_mo_last;
448
src_orbspi = moinfo.backtr_mo_orbspi;
450
src_ntri = (src_orbs * (src_orbs + 1)) / 2;
452
if(!strcmp(params.ref,"UHF")) {
454
/* we need to re-construct the C matrices because we now need the */
455
/* columns corresponding to frozen virtual orbitals */
456
destruct_evects(moinfo.backtr_nirreps, moinfo.evects_alpha);
457
moinfo.evects_alpha = (double ***) malloc (1 * sizeof(double **));
458
moinfo.evects_alpha[0] = block_matrix(moinfo.nao, moinfo.nmo);
460
destruct_evects(moinfo.backtr_nirreps, moinfo.evects_beta);
461
moinfo.evects_beta = (double ***) malloc (1 * sizeof(double **));
462
moinfo.evects_beta[0] = block_matrix(moinfo.nao, moinfo.nmo);
464
chkpt_init(PSIO_OPEN_OLD);
465
so2ao = chkpt_rd_usotao();
468
mmult(so2ao,1,moinfo.scf_vector_alpha,0,moinfo.evects_alpha[0],0,moinfo.nao,
469
moinfo.nso,moinfo.nmo,0);
470
mmult(so2ao,1,moinfo.scf_vector_beta,0,moinfo.evects_beta[0],0,moinfo.nao,
471
moinfo.nso,moinfo.nmo,0);
476
tmat = block_matrix(src_orbs, src_orbs);
477
psio_open(PSIF_MO_LAG, PSIO_OPEN_OLD);
479
psio_read_entry(PSIF_MO_LAG, "MO-basis Alpha Lagrangian", (char *) tmat[0],
480
sizeof(double)*src_orbs*src_orbs);
482
lag_a = init_array(src_ntri);
483
for (p=0; p<src_orbs; p++) {
484
for (q=0; q<=p; q++) {
485
P = moinfo.corr2pitz_nofzv_a[p];
486
Q = moinfo.corr2pitz_nofzv_a[q];
489
tval = 0.5 * (tmat[p][q] + tmat[q][p]);
491
if (params.lagran_double) tval *= 2.0;
492
if (params.lagran_halve) tval *= 0.5;
498
new_lag_a = init_array(dst_ntri);
499
tran_one(nirreps, moinfo.evects_alpha, src_orbs, src_first, src_last, src_orbspi, dst_orbs,
500
dst_first, dst_last, dst_orbspi, lag_a, new_lag_a, moinfo.order,
501
"\n\n\tAlpha Contribution to the Lagrangian (AO Basis):\n",
502
params.backtr, nfzc, print_integrals, outfile);
505
psio_read_entry(PSIF_MO_LAG, "MO-basis Beta Lagrangian", (char *) tmat[0],
506
sizeof(double)*src_orbs*src_orbs);
508
lag_b = init_array(src_ntri);
509
for (p=0; p<src_orbs; p++) {
510
for (q=0; q<=p; q++) {
511
P = moinfo.corr2pitz_nofzv_b[p];
512
Q = moinfo.corr2pitz_nofzv_b[q];
515
tval = 0.5 * (tmat[p][q] + tmat[q][p]);
517
if (params.lagran_double) tval *= 2.0;
518
if (params.lagran_halve) tval *= 0.5;
524
new_lag_b = init_array(dst_ntri);
525
tran_one(nirreps, moinfo.evects_beta, src_orbs, src_first, src_last, src_orbspi, dst_orbs,
526
dst_first, dst_last, dst_orbspi, lag_b, new_lag_b, moinfo.order,
527
"\n\n\tBeta Contribution to the Lagrangian (AO Basis):\n",
528
params.backtr, nfzc, print_integrals, outfile);
531
psio_close(PSIF_MO_LAG, 1);
534
/* combine the alpha and beta contributions */
535
for(p=0; p < dst_ntri; p++)
536
new_lag_a[p] += new_lag_b[p];
538
if(print_integrals) {
539
fprintf(outfile, "\n\tTotal AO-basis Lagrangian:\n");
540
print_array(new_lag_a, dst_orbs, outfile);
543
/* append the AO Lagrangian to the AO one-particle density matrix */
544
psio_open(itap, PSIO_OPEN_OLD);
545
psio_write_entry(itap, "AO-basis Lagrangian", (char *) new_lag_a, dst_ntri*sizeof(double));
553
/* we need to re-construct the C matrix because we now need the */
554
/* columns corresponding to frozen virtual orbitals */
555
destruct_evects(moinfo.backtr_nirreps, moinfo.evects);
556
moinfo.evects = (double ***) malloc (1 * sizeof(double **));
557
moinfo.evects[0] = block_matrix(moinfo.nao, moinfo.nmo);
558
chkpt_init(PSIO_OPEN_OLD);
559
so2ao = chkpt_rd_usotao();
561
mmult(so2ao,1,moinfo.scf_vector,0,moinfo.evects[0],0,moinfo.nao,moinfo.nso,
563
if (params.print_mos) {
564
fprintf(outfile, "C matrix (including AO to SO)\n");
565
print_mat(moinfo.evects[0], moinfo.nao, moinfo.nmo, outfile);
568
tmat = block_matrix(src_orbs, src_orbs);
571
/* read in the MO Lagrangian */
572
opdm = init_array(src_ntri);
573
psio_open(params.lag_in_file, PSIO_OPEN_OLD);
574
psio_read_entry(params.lag_in_file, "MO-basis Lagrangian", (char *) tmat[0],
575
src_orbs*src_orbs*sizeof(double));
576
psio_close(params.lag_in_file, 1);
577
if (params.print_lvl > 3) {
578
fprintf(outfile, "Lagrangian (MO Basis):\n");
579
print_mat(tmat, src_orbs, src_orbs, outfile);
582
/* reorder the Lagrangian to the back-transform order
583
* and enforce permutational symmetry.
585
for (p=0; p<src_orbs; p++) {
586
for (q=0; q<=p; q++) {
587
P = moinfo.corr2pitz[p];
588
Q = moinfo.corr2pitz[q];
589
tval = 0.5 * (tmat[p][q] + tmat[q][p]);
590
if (params.lagran_double) tval *= 2.0;
591
if (params.lagran_halve) tval *= 0.5;
598
if (params.print_lvl > 3) {
599
fprintf(outfile, "Reordered, Symmetrized Lagrangian in MO basis\n");
600
print_array(opdm, src_orbs, outfile);
603
/* now do the backtransformation */
604
tran_one(nirreps, moinfo.evects, src_orbs, src_first, src_last, src_orbspi, dst_orbs,
605
dst_first, dst_last, dst_orbspi, opdm, new_oe_ints, moinfo.order,
606
"\n\n\tLagrangian (AO Basis):\n",
607
params.backtr, nfzc, print_integrals, outfile);
609
psio_open(itap, PSIO_OPEN_OLD);
610
psio_write_entry(itap, "AO-basis Lagrangian", (char *) new_oe_ints, dst_ntri*sizeof(double));
618
if (params.print_lvl) {
619
fprintf(outfile, "\tOne-pdm and lagrangian written to file%d.\n", itap);
627
void tran_one(int nirreps, double ***C, int src_orbs, int *src_first, int *src_last,
628
int *src_orbspi, int dst_orbs, int *dst_first, int *dst_last,
629
int *dst_orbspi, double *src_ints, double *dst_ints, int *order,
630
char *label, int backtran, int nfzc, int printflg, FILE *outfile)
633
int psym, p, q, P, Q, pfirst, plast, pq;
634
int i, j, ifirst, ilast, I, J, i2, j2, ij;
636
int A_cols, B_cols, *C_cols;
638
A_cols = MAX0(src_orbs,dst_orbs);
639
A = block_matrix(A_cols, A_cols);
640
B = block_matrix(src_orbs,dst_orbs);
642
C_cols = backtran ? src_orbspi : dst_orbspi;
645
fprintf(outfile, "%s\n", label);
648
for (psym=0; psym < nirreps; psym++) {
649
pfirst = src_first[psym];
650
plast = src_last[psym];
651
for (p=pfirst,P=0; p <= plast; p++,P++) {
652
for (q=pfirst,Q=0; q <= plast; q++,Q++) {
654
A[P][Q] = src_ints[pq];
658
if (src_orbspi[psym] && dst_orbspi[psym]) {
661
C_DGEMM('n', backtran ? 't' : 'n', src_orbspi[psym], dst_orbspi[psym],
662
src_orbspi[psym], 1.0, A[0], A_cols, C[psym][0], C_cols[psym],
665
C_DGEMM(backtran ? 'n' : 't', 'n', dst_orbspi[psym], dst_orbspi[psym],
666
src_orbspi[psym], 1.0, C[psym][0], C_cols[psym], B[0], B_cols,
669
mmult(A,0,C[psym],(backtran ? 1 : 0),B,0,
670
src_orbspi[psym],src_orbspi[psym],dst_orbspi[psym],0);
671
mmult(C[psym],(backtran ? 0 : 1),B,0,A,0,
672
dst_orbspi[psym],src_orbspi[psym],dst_orbspi[psym],0);
678
if (dst_orbspi[psym]) {
679
fprintf(outfile, " Irrep %s\n", moinfo.labels[psym]);
680
print_mat(A,dst_orbspi[psym],dst_orbspi[psym],outfile);
684
ifirst = dst_first[psym];
685
ilast = dst_last[psym];
686
for (i=ifirst,I=0; i <= ilast; i++,I++) {
687
for (j=ifirst,J=0; (j <= ilast) && (j <= i); j++, J++) {
697
dst_ints[ij] = A[I][J];
700
} /* end loop over irreps */