← Back to branch summary

~ubuntu-branches/ubuntu/trusty/psicode/trusty

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/ccdensity.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/ccdensity.cc
 
1
/*! \file
 
2
    \ingroup CCDENSITY
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
/*
 
6
**  CCDENSITY: Program to calculate the coupled-cluster one- and 
 
7
**             two-particle densities.
 
8
*/
 
9
 
 
10
#include <cstdio>
 
11
#include <cstdlib>
 
12
#include <cstring>
 
13
#include <libipv1/ip_lib.h>
 
14
#include <libpsio/psio.h>
 
15
#include <libciomr/libciomr.h>
 
16
#include <libdpd/dpd.h>
 
17
#include <libiwl/iwl.h>
 
18
#include <cmath>
 
19
#include <psifiles.h>
 
20
#include "MOInfo.h"
 
21
#include "Params.h"
 
22
#include "Frozen.h"
 
23
#include "globals.h"
 
24
 
 
25
namespace psi { namespace ccdensity {
 
26
 
 
27
void init_io(int argc, char *argv[]);
 
28
void title(void);
 
29
void get_moinfo(void);
 
30
void get_frozen(void);
 
31
void get_params(void);
 
32
void exit_io(void);
 
33
void onepdm(struct RHO_Params);
 
34
void sortone(struct RHO_Params);
 
35
void twopdm(void);
 
36
void energy(struct RHO_Params);
 
37
void resort_tei(void);
 
38
void resort_gamma(void);
 
39
void lag(struct RHO_Params rho_params);
 
40
void build_X(void);
 
41
void build_A(void);
 
42
void build_Z(void);
 
43
void relax_I(void);
 
44
void relax_D(struct RHO_Params rho_params);
 
45
void sortI(void);
 
46
void fold(struct RHO_Params rho_params);
 
47
void deanti(struct RHO_Params rho_params);
 
48
void add_ref_RHF(struct iwlbuf *);
 
49
void add_ref_ROHF(struct iwlbuf *);
 
50
void add_ref_UHF(struct iwlbuf *, struct iwlbuf *, struct iwlbuf *);
 
51
void add_core_ROHF(struct iwlbuf *);
 
52
void add_core_UHF(struct iwlbuf *, struct iwlbuf *, struct iwlbuf *);
 
53
void dump_RHF(struct iwlbuf *, struct RHO_Params rho_params);
 
54
void dump_ROHF(struct iwlbuf *, struct RHO_Params rho_params);
 
55
void dump_UHF(struct iwlbuf *, struct iwlbuf *, struct iwlbuf *, struct RHO_Params rho_params);
 
56
void kinetic(void);
 
57
void dipole(void);
 
58
void probable(void);
 
59
int **cacheprep_rhf(int level, int *cachefiles);
 
60
int **cacheprep_uhf(int level, int *cachefiles);
 
61
void cachedone_rhf(int **cachelist);
 
62
void cachedone_uhf(int **cachelist);
 
63
void setup_LR(struct RHO_Params);
 
64
void G_build(void);
 
65
void x_oe_intermediates(struct RHO_Params);
 
66
void x_onepdm(struct RHO_Params);
 
67
void x_te_intermediates(void);
 
68
void x_Gijkl(void);
 
69
void x_Gabcd(void);
 
70
void x_Gibja(void);
 
71
void x_Gijka(void);
 
72
void x_Gijab(void);
 
73
void x_Gciab(void);
 
74
void V_build_x(void);
 
75
void x_xi1(void);
 
76
void x_xi_zero(void);
 
77
void x_xi2(void);
 
78
void x_xi_oe_intermediates(void);
 
79
void G_norm(void);
 
80
void zero_onepdm(struct RHO_Params rho_params);
 
81
void zero_twopdm(void);
 
82
void get_rho_params(void);
 
83
void get_td_params(void);
 
84
void td_setup(struct TD_Params S);
 
85
void tdensity(struct TD_Params S);
 
86
void td_print(void);
 
87
void oscillator_strength(struct TD_Params *S);
 
88
void rotational_strength(struct TD_Params *S);
 
89
void ael(struct RHO_Params *rho_params);
 
90
void cleanup(void);
 
91
void td_cleanup(void);
 
92
void x_oe_intermediates_rhf(struct RHO_Params rho_params);
 
93
void x_te_intermediates_rhf(void);
 
94
void x_xi_intermediates(void);
 
95
void V_build(void);
 
96
 
 
97
}} // namespace psi::cdensity 
 
98
 
 
99
using namespace psi::ccdensity;
 
100
 
 
101
int main(int argc, char *argv[])
 
102
{
 
103
  int i;
 
104
  int **cachelist, *cachefiles;
 
105
  struct iwlbuf OutBuf;
 
106
  struct iwlbuf OutBuf_AA, OutBuf_BB, OutBuf_AB;
 
107
  dpdfile2 D;
 
108
  double tval;
 
109
  
 
110
  init_io(argc,argv);
 
111
  title();
 
112
  /*  get_frozen(); */
 
113
  get_params();
 
114
  get_moinfo();
 
115
  get_rho_params();
 
116
 
 
117
  if ((moinfo.nfzc || moinfo.nfzv) && params.relax_opdm) {
 
118
    fprintf(outfile, "\n\tGradients/orbital relaxation involving frozen orbitals not yet available.\n");
 
119
    exit(PSI_RETURN_FAILURE);
 
120
  }
 
121
 
 
122
  cachefiles = init_int_array(PSIO_MAXUNIT);
 
123
 
 
124
  if(params.ref == 0 || params.ref == 1) { /** RHF or ROHF **/
 
125
    cachelist = cacheprep_rhf(params.cachelev, cachefiles);
 
126
 
 
127
    dpd_init(0, moinfo.nirreps, params.memory, 0, cachefiles, cachelist, NULL, 
 
128
             2, moinfo.occpi, moinfo.occ_sym, moinfo.virtpi, moinfo.vir_sym);
 
129
 
 
130
  }
 
131
  else if(params.ref == 2) { /** UHF **/
 
132
    cachelist = cacheprep_uhf(params.cachelev, cachefiles);
 
133
 
 
134
    dpd_init(0, moinfo.nirreps, params.memory, 0, cachefiles, 
 
135
             cachelist, NULL, 4, moinfo.aoccpi, moinfo.aocc_sym, moinfo.avirtpi,
 
136
             moinfo.avir_sym, moinfo.boccpi, moinfo.bocc_sym, moinfo.bvirtpi, moinfo.bvir_sym);
 
137
  }
 
138
 
 
139
  for (i=0; i<params.nstates; ++i) {
 
140
 
 
141
    /* CC_GLG will contain L, or R0*L + Zeta, if relaxed and zeta is available */
 
142
    /* CC_GL will contain L */
 
143
    setup_LR(rho_params[i]);
 
144
 
 
145
    /* Calculate Xi, put Xi in EOM_XI, and quit */
 
146
    if ( params.calc_xi ) {
 
147
      /* these intermediates go into EOM_TMP and are used to compute Xi;
 
148
         they may be reused to compute the excited-state density matrix */
 
149
      if (params.ref == 0) {
 
150
        x_oe_intermediates_rhf(rho_params[i]);
 
151
        x_te_intermediates_rhf();
 
152
      }
 
153
      else {
 
154
        x_oe_intermediates(rho_params[i]);
 
155
        x_te_intermediates();
 
156
      }
 
157
      x_xi_intermediates(); /*Xi intermediates put in EOM_TMP_XI */
 
158
      x_xi_zero(); /* make blank Xi */
 
159
      x_xi1();
 
160
      x_xi2();
 
161
      dpd_close(0);
 
162
      if(params.ref == 2) cachedone_uhf(cachelist);
 
163
      else cachedone_rhf(cachelist);
 
164
      free(cachefiles);
 
165
      cleanup();
 
166
      psio_close(EOM_TMP_XI,0); /* delete EOM_TMP_XI */
 
167
      psio_open(EOM_TMP_XI,PSIO_OPEN_NEW);
 
168
      exit_io();
 
169
      exit(PSI_RETURN_SUCCESS);
 
170
    }
 
171
 
 
172
    /* compute ground state parts of onepdm or put zeroes there */
 
173
    if ( ((rho_params[i].L_irr == rho_params[i].G_irr) || (params.use_zeta)) ) {
 
174
      zero_onepdm(rho_params[i]);
 
175
      onepdm(rho_params[i]);
 
176
    }
 
177
    else
 
178
      zero_onepdm(rho_params[i]);
 
179
 
 
180
    /* if the one-electron excited-state intermediates are not already on disk (from a Xi
 
181
       calculation, compute them.  They are nearly all necessary to compute the excited-state
 
182
       onepdm. Then complete excited-state onepdm.*/
 
183
    if (!rho_params[i].R_ground) {
 
184
      x_oe_intermediates(rho_params[i]); /* change to x_oe_intermediates_rhf() when rho gets spin-adapted */
 
185
      x_onepdm(rho_params[i]);
 
186
    }
 
187
 
 
188
    /* begin construction of twopdm */
 
189
    if (!params.onepdm) {
 
190
 
 
191
      /* Compute intermediates for construction of ground-state twopdm */
 
192
      if ( (params.L_irr == params.G_irr) || (params.use_zeta) ) {
 
193
        V_build(); /* uses CC_GLG, writes tau2*L2 to CC_MISC */
 
194
        G_build(); /* uses CC_GLG, writes t2*L2 to CC_GLG */
 
195
      }
 
196
 
 
197
      /* Compute ground-state twopdm or ground-state-like contributions to the excited twodpm */
 
198
      if ( (params.L_irr == params.G_irr) || (params.use_zeta) )
 
199
        twopdm();
 
200
      else
 
201
        zero_twopdm();
 
202
 
 
203
      /* Compute intermediates for construction of excited-state twopdm */
 
204
      if (!params.ground) {
 
205
        x_te_intermediates(); /* change to x_te_intermediates_rhf() when rho gets spin-adapted */
 
206
        V_build_x(); /* uses CC_GL, writes t2*L2 to EOM_TMP */
 
207
 
 
208
        /* add in non-R0 parts of onepdm and twopdm */
 
209
        x_Gijkl();
 
210
        x_Gabcd();
 
211
        x_Gibja();
 
212
        x_Gijka();
 
213
        x_Gciab();
 
214
        x_Gijab();
 
215
      }
 
216
    }
 
217
 
 
218
    sortone(rho_params[i]); /* puts full 1-pdm into moinfo.opdm */
 
219
    /* dipole(); */
 
220
 
 
221
    if (!params.onepdm) {
 
222
      if(!params.aobasis) energy(rho_params[i]);
 
223
 
 
224
      kinetic(); /* puts kinetic energy integrals into MO basis */
 
225
 
 
226
      lag(rho_params[i]); /* builds the orbital lagrangian pieces, I */
 
227
 
 
228
      /* dpd_init(1, moinfo.nirreps, params.memory, 2, frozen.occpi, frozen.occ_sym,
 
229
         frozen.virtpi, frozen.vir_sym); */
 
230
 
 
231
      /*  if(moinfo.nfzc || moinfo.nfzv) {
 
232
          resort_gamma();
 
233
          resort_tei();
 
234
          } */
 
235
 
 
236
      build_X(); /* builds orbital rotation gradient X */
 
237
      build_A(); /* construct MO Hessian A */
 
238
      build_Z(); /* solves the orbital Z-vector equations */
 
239
 
 
240
      relax_I(); /* adds orbital response contributions to Lagrangian */
 
241
 
 
242
      if (params.relax_opdm) {
 
243
        relax_D(rho_params[i]); /* adds orbital response contributions to onepdm */
 
244
      }
 
245
      sortone(rho_params[i]); /* builds large moinfo.opdm matrix */
 
246
      sortI(); /* builds large lagrangian matrix I */
 
247
      fold(rho_params[i]);
 
248
      deanti(rho_params[i]);
 
249
    }
 
250
 
 
251
    /*  dpd_close(0); dpd_close(1); */
 
252
 
 
253
    if(params.ref == 0) { /** RHF **/
 
254
 
 
255
      iwl_buf_init(&OutBuf, PSIF_MO_TPDM, params.tolerance, 0, 0);
 
256
 
 
257
      add_core_ROHF(&OutBuf);
 
258
      add_ref_RHF(&OutBuf);
 
259
      dump_RHF(&OutBuf, rho_params[i]);
 
260
 
 
261
      iwl_buf_flush(&OutBuf, 1);
 
262
      iwl_buf_close(&OutBuf, 1);
 
263
    }
 
264
    else if(params.ref == 1) { /** ROHF **/
 
265
 
 
266
      iwl_buf_init(&OutBuf, PSIF_MO_TPDM, params.tolerance, 0, 0);
 
267
 
 
268
      add_core_ROHF(&OutBuf);
 
269
      add_ref_ROHF(&OutBuf);
 
270
      dump_ROHF(&OutBuf, rho_params[i]);
 
271
 
 
272
      iwl_buf_flush(&OutBuf, 1);
 
273
      iwl_buf_close(&OutBuf, 1);
 
274
    }
 
275
    else if(params.ref == 2) { /** UHF **/
 
276
 
 
277
      iwl_buf_init(&OutBuf_AA, PSIF_MO_AA_TPDM, params.tolerance, 0, 0);
 
278
      iwl_buf_init(&OutBuf_BB, PSIF_MO_BB_TPDM, params.tolerance, 0, 0);
 
279
      iwl_buf_init(&OutBuf_AB, PSIF_MO_AB_TPDM, params.tolerance, 0, 0);
 
280
 
 
281
      /*    add_core_UHF(&OutBuf_AA, &OutBuf_BB, &OutBuf_AB); */
 
282
      add_ref_UHF(&OutBuf_AA, &OutBuf_BB, &OutBuf_AB);
 
283
      dump_UHF(&OutBuf_AA, &OutBuf_BB, &OutBuf_AB, rho_params[i]);
 
284
 
 
285
      iwl_buf_flush(&OutBuf_AA, 1);
 
286
      iwl_buf_flush(&OutBuf_BB, 1);
 
287
      iwl_buf_flush(&OutBuf_AB, 1);
 
288
      iwl_buf_close(&OutBuf_AA, 1);
 
289
      iwl_buf_close(&OutBuf_BB, 1);
 
290
      iwl_buf_close(&OutBuf_AB, 1);
 
291
    }
 
292
    free_block(moinfo.opdm);
 
293
 
 
294
    psio_close(CC_TMP,0);   psio_open(CC_TMP,PSIO_OPEN_NEW);
 
295
    psio_close(EOM_TMP0,0); psio_open(EOM_TMP0,PSIO_OPEN_NEW);
 
296
    psio_close(EOM_TMP1,0); psio_open(EOM_TMP1,PSIO_OPEN_NEW);
 
297
    psio_close(CC_GLG,0);   psio_open(CC_GLG,PSIO_OPEN_NEW);
 
298
    psio_close(CC_GL,0);    psio_open(CC_GL,PSIO_OPEN_NEW);
 
299
    psio_close(CC_GR,0);    psio_open(CC_GR,PSIO_OPEN_NEW);
 
300
    if (!params.calc_xi) {
 
301
      psio_close(EOM_TMP,0);
 
302
      psio_open(EOM_TMP,PSIO_OPEN_NEW);
 
303
    }
 
304
  }
 
305
 
 
306
  /* if ( params.ael && (params.nstates > 1) ) 
 
307
    ael(rho_params); */
 
308
 
 
309
  if(params.transition) {
 
310
 
 
311
    get_td_params();
 
312
    for(i=0; i < params.nstates; i++) {
 
313
      td_setup(td_params[i]);
 
314
      tdensity(td_params[i]);
 
315
      oscillator_strength(&(td_params[i]));
 
316
      if(params.ref == 0) {
 
317
        rotational_strength(&(td_params[i]));
 
318
      }
 
319
      td_cleanup();
 
320
    }
 
321
    td_print();
 
322
  }
 
323
 
 
324
  dpd_close(0);
 
325
 
 
326
  if(params.ref == 2) cachedone_uhf(cachelist);
 
327
  else cachedone_rhf(cachelist);
 
328
  free(cachefiles);
 
329
 
 
330
  cleanup(); 
 
331
  exit_io();
 
332
  exit(PSI_RETURN_SUCCESS);
 
333
}
 
334
 
 
335
extern "C" {const char *gprgid() { const char *prgid = "CCDENSITY"; return(prgid); }}
 
336
 
 
337
namespace psi { namespace ccdensity {
 
338
 
 
339
void init_io(int argc, char *argv[])
 
340
{
 
341
  int i, num_unparsed;
 
342
  char *progid, *argv_unparsed[100];;
 
343
 
 
344
  progid = (char *) malloc(strlen(gprgid())+2);
 
345
  sprintf(progid, ":%s",gprgid());
 
346
 
 
347
  params.onepdm = 0;
 
348
  params.prop_all = 0;
 
349
  params.calc_xi = 0;
 
350
  params.restart = 0;
 
351
  params.use_zeta = 0;
 
352
  params.transition = 0;
 
353
 
 
354
  for (i=1, num_unparsed=0; i<argc; ++i) {
 
355
    if(!strcmp(argv[i], "--onepdm")) {
 
356
      params.onepdm = 1; /* generate ONLY the onepdm (for one-electron properties) */
 
357
    }
 
358
    else if (!strcmp(argv[i],"--use_zeta")) {
 
359
      params.use_zeta = 1;
 
360
      params.ground = 0;
 
361
      params.restart = 1;
 
362
    }
 
363
    else if (!strcmp(argv[i],"--calc_xi")) {
 
364
      params.calc_xi = 1;
 
365
      params.ground = 0;
 
366
      params.restart = 0;
 
367
    }
 
368
    else if (!strcmp(argv[i],"--prop_all")) {
 
369
      params.prop_all = 1;
 
370
    }
 
371
    else if (!strcmp(argv[i],"--transition")) {
 
372
      params.transition = 1;
 
373
      params.relax_opdm = 0;
 
374
      params.ground = 0;
 
375
    }
 
376
    else {
 
377
      argv_unparsed[num_unparsed++] = argv[i];
 
378
    }
 
379
  }
 
380
 
 
381
  psi_start(&infile,&outfile,&psi_file_prefix,num_unparsed,argv_unparsed,0);
 
382
  ip_cwk_add(progid);
 
383
  free(progid);
 
384
  tstart(outfile);
 
385
  psio_init(); psio_ipv1_config();
 
386
 
 
387
  /* Open all dpd data files here */
 
388
  /* erase files for easy debugging */
 
389
  for(i=CC_MIN; i <= CC_MAX; i++) psio_open(i,PSIO_OPEN_OLD);
 
390
  /*
 
391
  psio_close(CC_GR,0);
 
392
  psio_close(CC_GL,0);
 
393
  psio_close(EOM_TMP0,0);
 
394
  psio_open(CC_GR,PSIO_OPEN_NEW);
 
395
  psio_open(CC_GL,PSIO_OPEN_NEW);
 
396
  psio_open(EOM_TMP0,PSIO_OPEN_NEW);
 
397
  */
 
398
}
 
399
 
 
400
void title(void)
 
401
{
 
402
  fprintf(outfile, "\n");
 
403
  fprintf(outfile, "\t\t\t**************************\n");
 
404
  fprintf(outfile, "\t\t\t*                        *\n");
 
405
  fprintf(outfile, "\t\t\t*        CCDENSITY       *\n");
 
406
  fprintf(outfile, "\t\t\t*                        *\n");
 
407
  fprintf(outfile, "\t\t\t**************************\n");
 
408
  fprintf(outfile, "\n");
 
409
}
 
410
 
 
411
void exit_io(void)
 
412
{
 
413
  int i;
 
414
 
 
415
  /* delete temporary EOM files */
 
416
  psio_close(EOM_TMP0,0);
 
417
  psio_close(EOM_TMP1,0);
 
418
  psio_close(CC_GLG,0);
 
419
  psio_open(EOM_TMP0,PSIO_OPEN_NEW);
 
420
  psio_open(EOM_TMP1,PSIO_OPEN_NEW);
 
421
  psio_open(CC_GLG,PSIO_OPEN_NEW);
 
422
  if (!params.calc_xi) {
 
423
    psio_close(EOM_TMP,0);
 
424
    psio_open(EOM_TMP,PSIO_OPEN_NEW);
 
425
  }
 
426
  if (params.use_zeta) { /* we're done with Xi amplitudes */
 
427
    psio_close(EOM_XI,0);
 
428
    psio_open(EOM_XI,PSIO_OPEN_NEW);
 
429
  }
 
430
 
 
431
  /* Close all dpd data files here */
 
432
  for(i=CC_MIN; i <= CC_MAX; i++) psio_close(i,1);
 
433
 
 
434
  psio_done();
 
435
  tstop(outfile);
 
436
  psi_stop(infile,outfile,psi_file_prefix);
 
437
}
 
438
 
 
439
}} // namespace psi::ccdensity