← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/get_rho_params.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/ccdensity/get_rho_params.cc
 
1
/*! \file
 
2
    \ingroup CCDENSITY
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cstring>
 
8
#include <libdpd/dpd.h>
 
9
#include <libqt/qt.h>
 
10
#include <libchkpt/chkpt.h>
 
11
#include <libipv1/ip_lib.h>
 
12
#include "MOInfo.h"
 
13
#include "Params.h"
 
14
#include "Frozen.h"
 
15
#define EXTERN
 
16
#include "globals.h"
 
17
 
 
18
namespace psi { namespace ccdensity {
 
19
 
 
20
void get_rho_params(void)
 
21
{
 
22
  int i,j,k,l,prop_sym,prop_root, lambda_and_Ls=0, errcod, prop_all,cnt;
 
23
  char lbl[32];
 
24
  int *states_per_irrep;
 
25
 
 
26
  /* check WFN keyword in input */
 
27
  errcod = ip_string("WFN", &(params.wfn), 0);
 
28
  if ( !cc_excited(params.wfn) ) params.ground = 1;
 
29
  else params.ground = 0;
 
30
 
 
31
  /* setup propery variables for excited states */
 
32
  if (!params.ground) {
 
33
    chkpt_init(PSIO_OPEN_OLD);
 
34
    if (chkpt_rd_override_occ()) {
 
35
      states_per_irrep = chkpt_rd_statespi();
 
36
    }
 
37
    else {
 
38
      ip_count("STATES_PER_IRREP", &i, 0);
 
39
      states_per_irrep = (int *) malloc(moinfo.nirreps * sizeof(int));
 
40
      for (i=0;i<moinfo.nirreps;++i)
 
41
        errcod = ip_data("STATES_PER_IRREP","%d",&(states_per_irrep[i]),1,i);
 
42
    }
 
43
    chkpt_close();
 
44
 
 
45
    prop_all = 0;
 
46
    if (ip_exist("PROP_ALL",0)) ip_boolean("PROP_ALL",&prop_all,0);
 
47
    /* can also be turned on by command-line (at least for now) */
 
48
    if (params.prop_all) prop_all = 1;
 
49
    if (params.dertype == 1) prop_all = 0; /* don't do relaxed, multiple excited states */
 
50
    params.prop_all = prop_all;
 
51
 
 
52
    if (ip_exist("PROP_SYM",0)) {  /* read symmetry of state for properties */
 
53
      ip_data("PROP_SYM","%d",&(prop_sym),0);
 
54
      prop_sym -= 1;
 
55
      prop_sym = moinfo.sym^prop_sym;
 
56
    }
 
57
    else { /* just use last irrep of states requested for symmetry of states */
 
58
      for (i=0;i<moinfo.nirreps;++i) {
 
59
        if (states_per_irrep[i] > 0)
 
60
          prop_sym = i^moinfo.sym;
 
61
      }
 
62
    }
 
63
    if (ip_exist("PROP_ROOT",0)) { /* read prop_root */
 
64
      ip_data("PROP_ROOT","%d",&(prop_root),0);
 
65
      prop_root -= 1;
 
66
    }
 
67
    else { /* just use highest root, if you need only one of them */
 
68
      prop_root = states_per_irrep[prop_sym^moinfo.sym];
 
69
      prop_root -= 1;
 
70
    }
 
71
  }
 
72
 
 
73
  /* setup density parameters */
 
74
  if (params.ground) { /* just compute ground state density */
 
75
    params.nstates = 1;
 
76
    rho_params = (struct RHO_Params *) malloc(params.nstates * sizeof(struct RHO_Params));
 
77
    rho_params[0].L_irr = 0;
 
78
    rho_params[0].R_irr = 0;
 
79
    rho_params[0].L_root = -1;
 
80
    rho_params[0].R_root = -1;
 
81
    rho_params[0].L_ground = 1;
 
82
    rho_params[0].R_ground = 1;
 
83
    rho_params[0].R0 = 1;
 
84
    rho_params[0].L0 = 0;
 
85
    rho_params[0].cceom_energy = 0;
 
86
  }
 
87
  else if (params.calc_xi) { /* just compute Xi for excited-state density */
 
88
    params.nstates = 1;
 
89
    rho_params = (struct RHO_Params *) malloc(params.nstates * sizeof(struct RHO_Params));
 
90
    rho_params[0].L_irr = prop_sym;
 
91
    rho_params[0].R_irr = prop_sym;
 
92
    rho_params[0].L_root = prop_root;
 
93
    rho_params[0].R_root = prop_root;
 
94
    rho_params[0].L_ground = 0;
 
95
    rho_params[0].R_ground = 0;
 
96
  }
 
97
  else { /* excited state density(ies) are involved */
 
98
    if (params.prop_all)  { /* do all roots */
 
99
      params.nstates = 1;
 
100
      for(i=0; i<moinfo.nirreps; i++) {
 
101
//        cnt = 0;
 
102
//        ip_data("STATES_PER_IRREP","%d",&cnt, 1, i);
 
103
        params.nstates += states_per_irrep[i];
 
104
      }
 
105
      rho_params = (struct RHO_Params *) malloc(params.nstates * sizeof(struct RHO_Params));
 
106
      rho_params[0].L_irr = 0;
 
107
      rho_params[0].R_irr = 0;
 
108
      rho_params[0].L_root = -1;
 
109
      rho_params[0].R_root = -1;
 
110
      rho_params[0].L_ground = 1;
 
111
      rho_params[0].R_ground = 1;
 
112
 
 
113
      cnt = 0;
 
114
      for(i=0; i<moinfo.nirreps; i++) { /* loop over irrep of R */
 
115
//        ip_data("STATES_PER_IRREP","%d",&j, 1, i);
 
116
//        for (k=0; k<j; ++k) {
 
117
        for (k=0; k<states_per_irrep[i]; ++k) {
 
118
          ++cnt;
 
119
          rho_params[cnt].L_irr = i^moinfo.sym;
 
120
          rho_params[cnt].R_irr = i^moinfo.sym;
 
121
          rho_params[cnt].L_root = k;
 
122
          rho_params[cnt].R_root = k;
 
123
          rho_params[cnt].L_ground = 0;
 
124
          rho_params[cnt].R_ground = 0;
 
125
        }
 
126
      }
 
127
    }
 
128
    else { /* do only one root */
 
129
      params.nstates = 1;
 
130
      rho_params = (struct RHO_Params *) malloc(params.nstates * sizeof(struct RHO_Params));
 
131
      rho_params[0].L_irr = prop_sym;
 
132
      rho_params[0].R_irr = prop_sym;
 
133
      rho_params[0].L_root = prop_root;
 
134
      rho_params[0].R_root = prop_root;
 
135
      rho_params[0].L_ground = 0;
 
136
      rho_params[0].R_ground = 0;
 
137
    }
 
138
    if(params.onepdm) params.transition = 1;
 
139
  }
 
140
 
 
141
  /* for each state, determine G_irr, cceom_energy, R0, and labels for files */
 
142
  for(i=0; i<params.nstates; i++) {
 
143
    rho_params[i].G_irr = (rho_params[i].L_irr) ^ (rho_params[i].R_irr);
 
144
 
 
145
    if (rho_params[i].R_ground) {
 
146
      rho_params[i].cceom_energy = 0.0;
 
147
      rho_params[i].R0 = 1.0;
 
148
    }
 
149
    else {
 
150
      if(!strcmp(params.wfn,"EOM_CC2")) {
 
151
        sprintf(lbl,"EOM CC2 Energy for root %d %d", rho_params[i].R_irr, rho_params[i].R_root);
 
152
        psio_read_entry(CC_INFO,lbl,(char*)&(rho_params[i].cceom_energy), sizeof(double));
 
153
        sprintf(lbl,"EOM CC2 R0 for root %d %d",rho_params[i].R_irr, rho_params[i].R_root);
 
154
        psio_read_entry(CC_INFO,lbl,(char*)&(rho_params[i].R0),sizeof(double));
 
155
      }
 
156
      else if(!strcmp(params.wfn,"EOM_CCSD")) {
 
157
        sprintf(lbl,"EOM CCSD Energy for root %d %d", rho_params[i].R_irr, rho_params[i].R_root);
 
158
        psio_read_entry(CC_INFO,lbl,(char*)&(rho_params[i].cceom_energy), sizeof(double));
 
159
        sprintf(lbl,"EOM CCSD R0 for root %d %d",rho_params[i].R_irr, rho_params[i].R_root);
 
160
        psio_read_entry(CC_INFO,lbl,(char*)&(rho_params[i].R0),sizeof(double));
 
161
      }
 
162
      else if(!strcmp(params.wfn,"EOM_CC3")) {
 
163
        sprintf(lbl,"EOM CC3 Energy for root %d %d", rho_params[i].R_irr, rho_params[i].R_root);
 
164
        psio_read_entry(CC_INFO,lbl,(char*)&(rho_params[i].cceom_energy), sizeof(double));
 
165
        sprintf(lbl,"EOM CC3 R0 for root %d %d",rho_params[i].R_irr, rho_params[i].R_root);
 
166
        psio_read_entry(CC_INFO,lbl,(char*)&(rho_params[i].R0),sizeof(double));
 
167
      }
 
168
    }
 
169
    if (rho_params[i].L_ground)
 
170
      rho_params[i].L0 = 1.0;
 
171
    else
 
172
      rho_params[i].L0 = 0.0;
 
173
 
 
174
    sprintf(rho_params[i].R1A_lbl, "RIA %d %d",  rho_params[i].R_irr, rho_params[i].R_root);
 
175
    sprintf(rho_params[i].R1B_lbl, "Ria %d %d",  rho_params[i].R_irr, rho_params[i].R_root);
 
176
    sprintf(rho_params[i].R2AA_lbl,"RIJAB %d %d",rho_params[i].R_irr, rho_params[i].R_root);
 
177
    sprintf(rho_params[i].R2BB_lbl,"Rijab %d %d",rho_params[i].R_irr, rho_params[i].R_root);
 
178
    sprintf(rho_params[i].R2AB_lbl,"RIjAb %d %d",rho_params[i].R_irr, rho_params[i].R_root);
 
179
    sprintf(rho_params[i].L1A_lbl, "LIA %d %d",  rho_params[i].L_irr, rho_params[i].L_root);
 
180
    sprintf(rho_params[i].L1B_lbl, "Lia %d %d",  rho_params[i].L_irr, rho_params[i].L_root);
 
181
    sprintf(rho_params[i].L2AA_lbl,"LIJAB %d %d",rho_params[i].L_irr, rho_params[i].L_root);
 
182
    sprintf(rho_params[i].L2BB_lbl,"Lijab %d %d",rho_params[i].L_irr, rho_params[i].L_root);
 
183
    sprintf(rho_params[i].L2AB_lbl,"LIjAb %d %d",rho_params[i].L_irr, rho_params[i].L_root);
 
184
 
 
185
    sprintf(rho_params[i].DIJ_lbl, "DIJ %d", i-1); /* change to a different label ? */
 
186
    sprintf(rho_params[i].Dij_lbl, "Dij %d", i-1);
 
187
    sprintf(rho_params[i].DAB_lbl, "DAB %d", i-1);
 
188
    sprintf(rho_params[i].Dab_lbl, "Dab %d", i-1);
 
189
    sprintf(rho_params[i].DAI_lbl, "DAI %d", i-1);
 
190
    sprintf(rho_params[i].Dai_lbl, "Dai %d", i-1);
 
191
    sprintf(rho_params[i].DIA_lbl, "DIA %d", i-1);
 
192
    sprintf(rho_params[i].Dia_lbl, "Dia %d", i-1);
 
193
 
 
194
    if (params.ref == 0) {
 
195
      if (i == 0) sprintf(rho_params[i].opdm_lbl, "MO-basis OPDM");
 
196
      else sprintf(rho_params[i].opdm_lbl, "MO-basis OPDM Root %d", i);
 
197
    } 
 
198
    else if (params.ref == 1) { /* ROHF */
 
199
      if (i == 0) sprintf(rho_params[i].opdm_lbl, "MO-basis OPDM");
 
200
      else sprintf(rho_params[i].opdm_lbl, "MO-basis OPDM Root %d", i);
 
201
    }
 
202
    else if (params.ref == 2) { /* UHF */
 
203
      if (i == 0) {
 
204
        sprintf(rho_params[i].opdm_a_lbl, "MO-basis Alpha OPDM");
 
205
        sprintf(rho_params[i].opdm_b_lbl, "MO-basis Beta OPDM");
 
206
      }
 
207
      else {
 
208
        sprintf(rho_params[i].opdm_a_lbl, "MO-basis Alpha OPDM Root %d", i);
 
209
        sprintf(rho_params[i].opdm_b_lbl, "MO-basis Beta OPDM Root %d", i);
 
210
      }
 
211
    }
 
212
  }
 
213
 
 
214
  psio_write_entry(CC_INFO, "Num. of CC densities", (char *) &(params.nstates), sizeof(int));
 
215
 
 
216
  fprintf(outfile,"\tNumber of States = %-d\n",params.nstates);
 
217
 
 
218
  /*
 
219
  fprintf(outfile,"\n\tGround? L_root L_irr R_root R_irr G_irr     EOM Energy        R0\n");
 
220
  for(i=0; i<params.nstates; i++) {
 
221
    fprintf(outfile,"\t%5s %6d %7s %4d %7s %4s %15.10lf %12.8lf\n",
 
222
            (rho_params[i].R_ground ? "Yes":"No"),
 
223
            rho_params[i].L_root+1, moinfo.labels[rho_params[i].L_irr],
 
224
            rho_params[i].R_root+1, moinfo.labels[rho_params[i].R_irr],
 
225
            moinfo.labels[rho_params[i].G_irr],
 
226
            rho_params[i].cceom_energy,
 
227
            rho_params[i].R0);
 
228
  }
 
229
  */
 
230
 
 
231
  fprintf(outfile,"\n\tGround?  State     EOM Energy       R0\n");
 
232
  for(i=0; i<params.nstates; i++) {
 
233
    fprintf(outfile,"\t%5s     %d%3s %15.10lf %12.8lf\n",
 
234
            (rho_params[i].R_ground ? "Yes":"No"),
 
235
            rho_params[i].L_root+1, moinfo.labels[rho_params[i].L_irr],
 
236
            rho_params[i].cceom_energy,
 
237
            rho_params[i].R0);
 
238
  }
 
239
 
 
240
  /* set variables for xi and twopdm code in the old non-state specific structure */
 
241
  params.G_irr = rho_params[0].G_irr;
 
242
  params.R_irr = rho_params[0].R_irr;
 
243
  params.L_irr = rho_params[0].L_irr;
 
244
  params.R0 = rho_params[0].R0;
 
245
  params.L0 = rho_params[0].L0;
 
246
  params.cceom_energy = rho_params[0].cceom_energy;
 
247
 
 
248
  return;
 
249
}
 
250
 
 
251
}} // namespace psi::ccdensity