~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/transqt/transform_one.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
ImportĀ upstreamĀ versionĀ 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
#include <stdio.h>
 
2
#include <stdlib.h>
 
3
#include <libipv1/ip_lib.h>
 
4
#include <libciomr/libciomr.h>
 
5
#include <libpsio/psio.h>
 
6
#include <libqt/qt.h>
 
7
#include <libiwl/iwl.h>
 
8
#include <libchkpt/chkpt.h>
 
9
#include <psifiles.h>
 
10
#include "MOInfo.h"
 
11
#include "Params.h"
 
12
#include "globals.h"
 
13
 
 
14
 
 
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)))
 
18
 
 
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);
 
28
 
 
29
/*
 
30
** TRANSFORM_ONE(): This function does the SO to MO transformation of
 
31
**    the one-electron integrals.
 
32
**
 
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.
 
36
**
 
37
** C. David Sherrill
 
38
** 29 April 1998
 
39
**
 
40
*/
 
41
void transform_one()
 
42
{
 
43
  if (params.backtr) { 
 
44
    trans_one_backwards();
 
45
  }
 
46
  else {
 
47
    trans_one_forwards();
 
48
  }
 
49
 
 
50
}
 
51
 
 
52
 
 
53
 
 
54
/*
 
55
** TRANS_ONE_FORWARDS()
 
56
**
 
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.
 
65
**
 
66
** C. David Sherrill
 
67
** 29 April 1998
 
68
*/
 
69
void trans_one_forwards(void)
 
70
{
 
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;
 
74
  int print_integrals;
 
75
  double *oe_ints, *new_oe_ints;
 
76
 
 
77
  if (params.print_lvl) {
 
78
    fprintf(outfile, "\n\tTransforming one-electron integrals...\n");
 
79
    fflush(outfile);
 
80
  }
 
81
 
 
82
  nirreps = moinfo.nirreps;
 
83
  nmo = moinfo.nmo;
 
84
  nso = moinfo.nso;
 
85
  print_integrals = params.print_oe_ints;
 
86
  nfzc = 0;  nfzv = 0;
 
87
 
 
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;
 
94
 
 
95
  src_orbs = nso;
 
96
  dst_orbs = nmo - nfzv - nfzc;
 
97
  dst_ntri = (dst_orbs * (dst_orbs + 1)) / 2;
 
98
  src_ntri = (src_orbs * (src_orbs + 1)) / 2;
 
99
 
 
100
  new_oe_ints = init_array(dst_ntri);
 
101
 
 
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
 
107
   */ 
 
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);
 
112
    }
 
113
    else destruct_evects(nirreps, moinfo.evects);
 
114
 
 
115
    if(!strcmp(params.ref,"UHF")) {
 
116
      moinfo.evects_alpha = construct_evects("alpha", moinfo.nirreps, 
 
117
                                             moinfo.orbspi,
 
118
                                             moinfo.sopi, moinfo.orbspi,
 
119
                                             moinfo.first_so, moinfo.last_so,
 
120
                                             moinfo.first, moinfo.last,
 
121
                                             moinfo.first, moinfo.last, 
 
122
                                             params.print_mos);
 
123
      moinfo.evects_beta = construct_evects("beta", moinfo.nirreps, 
 
124
                                            moinfo.orbspi,
 
125
                                            moinfo.sopi, moinfo.orbspi,
 
126
                                            moinfo.first_so, moinfo.last_so,
 
127
                                            moinfo.first, moinfo.last,
 
128
                                            moinfo.first, moinfo.last, 
 
129
                                            params.print_mos);
 
130
    }
 
131
    else {
 
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, 
 
137
                                       params.print_mos);
 
138
    }
 
139
  }
 
140
 
 
141
  if (params.do_h_bare) {
 
142
    
 
143
    oe_ints = moinfo.oe_ints;
 
144
 
 
145
    if(!strcmp(params.ref,"UHF")) {
 
146
    
 
147
      itap = params.h_bare_a_file;
 
148
 
 
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);
 
154
 
 
155
      iwl_wrtone(itap, PSIF_MO_A_OEI, dst_ntri, new_oe_ints);
 
156
 
 
157
      if (params.print_lvl) {
 
158
        fprintf(outfile, "\tOne-electron A integrals written to file %d.\n", 
 
159
                itap);
 
160
        fflush(outfile);
 
161
      }
 
162
 
 
163
      itap = params.h_bare_b_file;
 
164
 
 
165
      /* Clear the new_oe_ints */
 
166
      zero_arr(new_oe_ints,dst_ntri);
 
167
 
 
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);
 
173
 
 
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", 
 
177
                itap);
 
178
        fflush(outfile);
 
179
      }
 
180
    }
 
181
    else {
 
182
 
 
183
      itap = params.h_bare_file;
 
184
 
 
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);
 
190
 
 
191
 
 
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);
 
195
        fflush(outfile);
 
196
      }
 
197
    }
 
198
 
 
199
 
 
200
  }
 
201
 
 
202
 
 
203
  if (params.do_h_fzc) {
 
204
 
 
205
    if(!strcmp(params.ref,"UHF")) {
 
206
 
 
207
      itap = params.h_fzc_a_file;
 
208
      oe_ints = moinfo.fzc_operator_alpha;
 
209
 
 
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);
 
215
 
 
216
      iwl_wrtone(itap, PSIF_MO_A_FZC, dst_ntri, new_oe_ints);
 
217
 
 
218
      if (params.print_lvl) {
 
219
        fprintf(outfile, "\tAlpha frozen-core operator written to file %d.\n", 
 
220
                itap);
 
221
        fflush(outfile);
 
222
      }
 
223
 
 
224
      /* Clear the new_oe_ints */
 
225
      zero_arr(new_oe_ints,dst_ntri);
 
226
 
 
227
      itap = params.h_fzc_b_file;
 
228
      oe_ints = moinfo.fzc_operator_beta;
 
229
 
 
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);
 
235
 
 
236
      iwl_wrtone(itap, PSIF_MO_B_FZC, dst_ntri, new_oe_ints);
 
237
 
 
238
      if (params.print_lvl) {
 
239
        fprintf(outfile, "\tBeta frozen-core operator written to file %d.\n", 
 
240
                itap);
 
241
        fflush(outfile);
 
242
      }
 
243
    }
 
244
    else {
 
245
 
 
246
      itap = params.h_fzc_file;
 
247
      oe_ints = moinfo.fzc_operator;
 
248
 
 
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);
 
254
 
 
255
      iwl_wrtone(itap, PSIF_MO_FZC, dst_ntri, new_oe_ints);
 
256
 
 
257
      if (params.print_lvl) {
 
258
        fprintf(outfile, "\tFrozen-core operator written to file %d.\n", itap);
 
259
        fflush(outfile);
 
260
      }
 
261
 
 
262
    }
 
263
 
 
264
  }
 
265
 
 
266
  free(new_oe_ints);
 
267
  
 
268
}
 
269
 
 
270
 
 
271
/*
 
272
** TRANS_ONE_BACKWARDS()
 
273
**
 
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().
 
279
**
 
280
** C. David Sherrill
 
281
** 29 April 1998
 
282
*/
 
283
void trans_one_backwards(void)
 
284
{
 
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;
 
288
  int print_integrals;
 
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;
 
293
  int p,q,pq;
 
294
  int I,J,P,Q,PQ;
 
295
  PSI_FPTR pdm_fp_in=0, pdm_fp_out=0;
 
296
 
 
297
  if (params.print_lvl) {
 
298
    fprintf(outfile, "\n\tTransforming one-electron integrals...\n");
 
299
    fflush(outfile);
 
300
  }
 
301
 
 
302
  itap = params.opdm_out_file;
 
303
  nirreps = moinfo.backtr_nirreps;
 
304
  nmo = moinfo.nmo;
 
305
  print_integrals = params.print_oe_ints;
 
306
  nfzc = 0;
 
307
  nfzv = moinfo.nfzv;   
 
308
  
 
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;
 
315
  
 
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;
 
320
 
 
321
  if(!strcmp(params.ref, "UHF")) {
 
322
 
 
323
    tmat = block_matrix(src_orbs, src_orbs);
 
324
    psio_open(PSIF_MO_OPDM, PSIO_OPEN_OLD);
 
325
 
 
326
    psio_read_entry(PSIF_MO_OPDM, "MO-basis Alpha OPDM", (char *) tmat[0], 
 
327
                    sizeof(double)*src_orbs*src_orbs);
 
328
 
 
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];
 
334
        PQ = INDEX(P,Q);
 
335
        opdm_a[PQ] = 0.5 * (tmat[p][q] + tmat[q][p]);
 
336
      }
 
337
    } 
 
338
 
 
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);
 
344
    free(opdm_a);
 
345
 
 
346
    psio_read_entry(PSIF_MO_OPDM, "MO-basis Beta OPDM", (char *) tmat[0], 
 
347
                    sizeof(double)*src_orbs*src_orbs);
 
348
 
 
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];
 
354
        PQ = INDEX(P,Q);
 
355
        opdm_b[PQ] = 0.5 * (tmat[p][q] + tmat[q][p]);
 
356
      }
 
357
    } 
 
358
 
 
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);
 
364
    free(opdm_b);
 
365
 
 
366
    psio_close(PSIF_MO_OPDM, 1);
 
367
    free_block(tmat);
 
368
 
 
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);
 
374
    }
 
375
 
 
376
    /* combine the alpha and beta contributions */
 
377
    for(p=0; p < dst_ntri; p++)
 
378
      new_opdm_a[p] += new_opdm_b[p];
 
379
 
 
380
    if(print_integrals) {
 
381
      fprintf(outfile, "\n\tTotal AO-basis OPDM:\n");
 
382
      print_array(new_opdm_a, dst_orbs, outfile);
 
383
    }
 
384
 
 
385
    psio_open(itap, PSIO_OPEN_OLD);
 
386
    psio_write_entry(itap, "AO-basis OPDM", (char *) new_opdm_a, dst_ntri*sizeof(double));
 
387
    psio_close(itap, 1);
 
388
    
 
389
    free(new_opdm_a);
 
390
    free(new_opdm_b);
 
391
 
 
392
  }
 
393
  else {
 
394
    /* read in the MO one-particle density matrix */
 
395
    opdm = init_array(src_ntri);
 
396
 
 
397
    tmat = block_matrix(src_orbs, src_orbs);
 
398
 
 
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));
 
403
    else {
 
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];
 
410
      free_block(tmat2);
 
411
    }
 
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);
 
416
    }
 
417
 
 
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];
 
423
        PQ = INDEX(P,Q);
 
424
        opdm[PQ] = 0.5 * (tmat[p][q] + tmat[q][p]);
 
425
      }
 
426
    } 
 
427
  
 
428
    oe_ints = opdm;
 
429
    new_oe_ints = init_array(dst_ntri);
 
430
  
 
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);
 
435
 
 
436
    psio_open(itap, PSIO_OPEN_OLD);
 
437
    psio_write_entry(itap, "AO-basis OPDM", (char *) new_oe_ints, dst_ntri*sizeof(double));
 
438
    psio_close(itap, 1);
 
439
    
 
440
    free(opdm);
 
441
    free_block(tmat);
 
442
  }  
 
443
 
 
444
 
 
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;
 
449
  src_orbs = nmo; 
 
450
  src_ntri = (src_orbs * (src_orbs + 1)) / 2;
 
451
 
 
452
  if(!strcmp(params.ref,"UHF")) {
 
453
 
 
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);
 
459
 
 
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);
 
463
 
 
464
    chkpt_init(PSIO_OPEN_OLD);
 
465
    so2ao = chkpt_rd_usotao();
 
466
    chkpt_close();
 
467
 
 
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);
 
472
 
 
473
    free_block(so2ao);
 
474
 
 
475
 
 
476
    tmat = block_matrix(src_orbs, src_orbs);
 
477
    psio_open(PSIF_MO_LAG, PSIO_OPEN_OLD);
 
478
 
 
479
    psio_read_entry(PSIF_MO_LAG, "MO-basis Alpha Lagrangian", (char *) tmat[0], 
 
480
                    sizeof(double)*src_orbs*src_orbs);
 
481
 
 
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];
 
487
        PQ = INDEX(P,Q);
 
488
 
 
489
        tval = 0.5 * (tmat[p][q] + tmat[q][p]);
 
490
 
 
491
        if (params.lagran_double) tval *= 2.0;
 
492
        if (params.lagran_halve) tval *= 0.5;
 
493
 
 
494
        lag_a[PQ] = tval;
 
495
      }
 
496
    } 
 
497
 
 
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);
 
503
    free(lag_a);
 
504
 
 
505
    psio_read_entry(PSIF_MO_LAG, "MO-basis Beta Lagrangian", (char *) tmat[0], 
 
506
                    sizeof(double)*src_orbs*src_orbs);
 
507
 
 
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];
 
513
        PQ = INDEX(P,Q);        
 
514
 
 
515
        tval = 0.5 * (tmat[p][q] + tmat[q][p]);
 
516
 
 
517
        if (params.lagran_double) tval *= 2.0;
 
518
        if (params.lagran_halve) tval *= 0.5;
 
519
 
 
520
        lag_b[PQ] = tval;
 
521
      }
 
522
    } 
 
523
 
 
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);
 
529
    free(lag_b);
 
530
 
 
531
    psio_close(PSIF_MO_LAG, 1);
 
532
    free_block(tmat);
 
533
 
 
534
    /* combine the alpha and beta contributions */
 
535
    for(p=0; p < dst_ntri; p++)
 
536
      new_lag_a[p] += new_lag_b[p];
 
537
 
 
538
    if(print_integrals) {
 
539
      fprintf(outfile, "\n\tTotal AO-basis Lagrangian:\n");
 
540
      print_array(new_lag_a, dst_orbs, outfile);
 
541
    }
 
542
 
 
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));
 
546
    psio_close(itap, 1);
 
547
    
 
548
    free(new_lag_a);
 
549
    free(new_lag_b);
 
550
 
 
551
  }
 
552
  else {  
 
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();
 
560
    chkpt_close();
 
561
    mmult(so2ao,1,moinfo.scf_vector,0,moinfo.evects[0],0,moinfo.nao,moinfo.nso,
 
562
          moinfo.nmo,0);
 
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);
 
566
    }
 
567
    free_block(so2ao);
 
568
    tmat = block_matrix(src_orbs, src_orbs);
 
569
 
 
570
 
 
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);
 
580
    }
 
581
  
 
582
    /* reorder the Lagrangian to the back-transform order 
 
583
     * and enforce permutational symmetry.  
 
584
     */
 
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;
 
592
        PQ = INDEX(P,Q);
 
593
        opdm[PQ] += tval;
 
594
      }
 
595
    } 
 
596
    free_block(tmat);
 
597
  
 
598
    if (params.print_lvl > 3) {
 
599
      fprintf(outfile, "Reordered, Symmetrized Lagrangian in MO basis\n");
 
600
      print_array(opdm, src_orbs, outfile);
 
601
    }
 
602
  
 
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);
 
608
 
 
609
    psio_open(itap, PSIO_OPEN_OLD);
 
610
    psio_write_entry(itap, "AO-basis Lagrangian", (char *) new_oe_ints, dst_ntri*sizeof(double));
 
611
    psio_close(itap, 1);
 
612
    free(opdm);
 
613
 
 
614
    free(new_oe_ints);
 
615
  
 
616
  }
 
617
 
 
618
  if (params.print_lvl) {
 
619
    fprintf(outfile, "\tOne-pdm and lagrangian written to file%d.\n", itap);
 
620
  }
 
621
 
 
622
}
 
623
 
 
624
 
 
625
 
 
626
 
 
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)
 
631
{
 
632
 
 
633
  int psym, p, q, P, Q, pfirst, plast, pq;
 
634
  int i, j, ifirst, ilast, I, J, i2, j2, ij;
 
635
  double **A, **B;
 
636
  int A_cols, B_cols, *C_cols;
 
637
  
 
638
  A_cols = MAX0(src_orbs,dst_orbs);
 
639
  A = block_matrix(A_cols, A_cols);
 
640
  B = block_matrix(src_orbs,dst_orbs);
 
641
  B_cols = dst_orbs;
 
642
  C_cols = backtran ? src_orbspi : dst_orbspi;
 
643
 
 
644
  if (printflg) {
 
645
    fprintf(outfile, "%s\n", label);
 
646
  }
 
647
 
 
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++) {
 
653
        pq = INDEX(p,q);
 
654
        A[P][Q] = src_ints[pq];
 
655
      }
 
656
    }
 
657
    
 
658
    if (src_orbspi[psym] && dst_orbspi[psym]) {
 
659
 
 
660
#ifdef USE_BLAS
 
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],
 
663
              0.0, B[0], B_cols);
 
664
 
 
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,
 
667
              0.0, A[0], A_cols);
 
668
#else
 
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);
 
673
#endif
 
674
 
 
675
    }    
 
676
 
 
677
    if (printflg) {
 
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);
 
681
      }
 
682
    }
 
683
 
 
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++) {
 
688
 
 
689
        if (!backtran) {
 
690
          i2 = order[i]-nfzc;
 
691
          j2 = order[j]-nfzc;
 
692
          ij = INDEX(i2,j2);
 
693
        }
 
694
        else
 
695
          ij = INDEX(i,j);
 
696
 
 
697
        dst_ints[ij] = A[I][J];
 
698
      }
 
699
    }
 
700
  } /* end loop over irreps */
 
701
 
 
702
  free_block(A);
 
703
  free_block(B);
 
704
 
 
705
}
 
706
 
 
707