← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/mp2/mp2.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/mp2/mp2.cc
 
1
/*! \defgroup MP2 mp2: Canonical Evaluation of MP2 Energy and Gradients */
 
2
 
 
3
/*! 
 
4
** \file
 
5
** \ingroup MP2
 
6
** \brief Canonical evaluation of MP2 energy and gradients
 
7
*/
 
8
 
 
9
#include <cstdio>
 
10
#include <cstring>
 
11
#include <cstdlib>
 
12
#include <cmath>
 
13
#include <libipv1/ip_lib.h>
 
14
#include <libciomr/libciomr.h>
 
15
#include <libdpd/dpd.h>
 
16
#include <libchkpt/chkpt.h>
 
17
#include <libqt/qt.h>
 
18
#include <libiwl/iwl.h>
 
19
#include <physconst.h>
 
20
#include <psifiles.h>
 
21
#include "globals.h"
 
22
 
 
23
namespace psi{ namespace mp2{
 
24
 
 
25
void init_io(int argc, char *argv[]);
 
26
void title(void);
 
27
void get_moinfo(void);
 
28
void get_params(void);
 
29
void init_ioff(void);
 
30
int **cacheprep_rhf(int level, int *cachefiles);
 
31
int **cacheprep_uhf(int level, int *cachefiles);
 
32
void cachedone_rhf(int **cachelist);
 
33
double energy(void);
 
34
void amps(void);
 
35
void sort_amps(void);
 
36
void opdm(void);
 
37
void twopdm(void);
 
38
void lag(void);
 
39
void build_X(void);
 
40
void build_A(void);
 
41
void Zvector(void);
 
42
void relax_I(void);
 
43
void relax_opdm(void);
 
44
void sort_opdm(void);
 
45
void sort_I(void);
 
46
void fold(void);
 
47
void deanti(void);
 
48
void write_data(void);
 
49
void check_energy(int);
 
50
void sort_twopdm(void);
 
51
void cleanup(void);
 
52
void exit_io(void);
 
53
 
 
54
}} /* End namespaces */
 
55
 
 
56
int main(int argc, char *argv[])
 
57
{
 
58
  using namespace psi::mp2;
 
59
 
 
60
  int *cachefiles;
 
61
  int **cachelist;
 
62
 
 
63
  init_io(argc,argv);
 
64
  title();
 
65
 
 
66
  get_moinfo();
 
67
  get_params();
 
68
  init_ioff();
 
69
  
 
70
  cachefiles = init_int_array(PSIO_MAXUNIT);
 
71
 
 
72
  if(params.ref == 2) { /** UHF **/
 
73
    cachelist = cacheprep_uhf(params.cachelev,cachefiles);
 
74
    dpd_init(0,mo.nirreps,params.memory,params.cachetype,cachefiles,cachelist,
 
75
             NULL,4,mo.aoccpi,mo.aocc_sym,mo.avirpi,mo.avir_sym,mo.boccpi,
 
76
             mo.bocc_sym,mo.bvirpi,mo.bvir_sym);
 
77
  }
 
78
  else { /** RHF or ROHF **/
 
79
    cachelist = cacheprep_rhf(params.cachelev,cachefiles);
 
80
    dpd_init(0,mo.nirreps,params.memory,params.cachetype,cachefiles,cachelist,
 
81
             NULL,2,mo.occpi,mo.occ_sym,mo.virpi,mo.vir_sym);
 
82
  }
 
83
  
 
84
  amps();
 
85
 
 
86
  mo.Emp2 = energy();
 
87
  
 
88
  fprintf(outfile,"\n");
 
89
  fprintf(outfile,"\tScaled_OS correlation energy      = %20.15f\n",mo.escsmp2_os);
 
90
  fprintf(outfile,"\tScaled_SS correlation energy      = %20.15f\n",mo.escsmp2_ss);
 
91
  fprintf(outfile,"\tSCS-MP2 correlation energy        = %20.15f\n",mo.escsmp2_os+mo.escsmp2_ss);
 
92
  fprintf(outfile,"      * SCS-MP2 total energy              = %20.15f\n",mo.Escf+mo.escsmp2_os+mo.escsmp2_ss);
 
93
 
 
94
  fprintf(outfile,"\n\tOpposite-spin correlation energy  = %20.15f\n",mo.emp2_os);
 
95
  fprintf(outfile,"\tSame-spin correlation energy      = %20.15f\n",mo.emp2_ss);
 
96
  fprintf(outfile,"\tMP2 correlation energy            = %20.15f\n",mo.Emp2);
 
97
  fprintf(outfile,"      * MP2 total energy                  = %20.15f\n\n",mo.Escf + mo.Emp2);
 
98
  fflush(outfile);
 
99
 
 
100
  chkpt_init(PSIO_OPEN_OLD);
 
101
  // Save MP2 contribution to Chkpt
 
102
  chkpt_wt_emp2(mo.Emp2);
 
103
  chkpt_wt_etot(mo.Escf+mo.Emp2);
 
104
  chkpt_close();
 
105
  
 
106
  if(params.opdm) {
 
107
    sort_amps();
 
108
    opdm();
 
109
    if(params.relax_opdm) {
 
110
      lag();
 
111
      build_A();
 
112
      Zvector();
 
113
    }
 
114
    sort_opdm();
 
115
    //dipole();
 
116
  }
 
117
 
 
118
  if(params.gradient) {
 
119
    sort_amps();
 
120
    opdm();
 
121
    twopdm();
 
122
    lag();
 
123
    build_X();
 
124
    build_A();
 
125
    Zvector();
 
126
    relax_I();
 
127
    relax_opdm();
 
128
    sort_I(); 
 
129
    sort_opdm();
 
130
    fold();
 
131
    deanti();
 
132
    write_data();
 
133
  }
 
134
 
 
135
  dpd_close(0);
 
136
  
 
137
  cleanup();
 
138
 
 
139
  exit_io();
 
140
  
 
141
  exit(0);
 
142
}
 
143
 
 
144
extern "C"{ const char *gprgid() { const char *prgid = "MP2"; return(prgid); } }
 
145
 
 
146
namespace psi{ namespace mp2{
 
147
 
 
148
void init_io(int argc, char *argv[])
 
149
{
 
150
  int i=0;
 
151
  int num_extra_args=0;
 
152
  char **extra_args;
 
153
  char *progid;
 
154
 
 
155
  extra_args = (char **) malloc(argc*sizeof(char *));
 
156
  progid = (char *) malloc(strlen(gprgid())+2);
 
157
  sprintf(progid, ":%s", gprgid());
 
158
  params.opdm = 0;
 
159
 
 
160
  for(i=1; i<argc; i++) {
 
161
    if(strcmp(argv[i], "--opdm") == 0) {
 
162
      params.opdm = 1;
 
163
    }
 
164
    else {
 
165
      extra_args[num_extra_args++] = argv[i];
 
166
    }
 
167
  }
 
168
  
 
169
  psi_start(&infile,&outfile,&psi_file_prefix,num_extra_args,extra_args,0);
 
170
  ip_cwk_add(progid);
 
171
  free(progid);
 
172
  tstart(outfile);
 
173
  
 
174
  psio_init(); psio_ipv1_config();
 
175
  for(i=CC_MIN; i <= CC_MAX; i++) 
 
176
    psio_open(i,1);
 
177
 
 
178
  free(extra_args);
 
179
}
 
180
 
 
181
void title(void)
 
182
{
 
183
  fprintf(outfile, "\t\t\t*************************\n");
 
184
  fprintf(outfile, "\t\t\t*                       *\n");
 
185
  fprintf(outfile, "\t\t\t*          MP2          *\n");
 
186
  fprintf(outfile, "\t\t\t*                       *\n");
 
187
  fprintf(outfile, "\t\t\t*************************\n");
 
188
  fflush(outfile);
 
189
}
 
190
 
 
191
void init_ioff(void)
 
192
{
 
193
  int i;
 
194
  
 
195
  ioff = init_int_array(MAXIOFF);
 
196
  ioff[0] = 0;
 
197
  for(i=1; i < MAXIOFF; i++) {
 
198
    ioff[i] = ioff[i-1] + i;
 
199
  }
 
200
  
 
201
}
 
202
 
 
203
void cleanup(void)
 
204
{
 
205
  int i;
 
206
  
 
207
  free(params.wfn);
 
208
  
 
209
  free(mo.doccpi);
 
210
  free(mo.soccpi);
 
211
  free(mo.mopi);
 
212
  free(mo.fzdoccpi);
 
213
  free(mo.fzvirtpi);
 
214
  for(i=0; i < mo.nirreps; i++)
 
215
    free(mo.irreplabels[i]);
 
216
  free(mo.irreplabels);
 
217
  free(ioff);
 
218
  
 
219
  if(params.ref == 2) {
 
220
    free(mo.aoccpi);
 
221
    free(mo.boccpi);
 
222
    free(mo.avirpi);
 
223
    free(mo.bvirpi);
 
224
    free(mo.aocc_sym);
 
225
    free(mo.bocc_sym);
 
226
    free(mo.avir_sym);
 
227
    free(mo.bvir_sym);
 
228
    free(mo.aocc_off);
 
229
    free(mo.bocc_off);
 
230
    free(mo.avir_off);
 
231
    free(mo.bvir_off);
 
232
    free(mo.qt_aocc);
 
233
    free(mo.qt_bocc);
 
234
    free(mo.qt_avir);
 
235
    free(mo.qt_bvir);
 
236
  }
 
237
  else {
 
238
    free(mo.occpi);
 
239
    free(mo.virpi);
 
240
    free(mo.occ_sym);
 
241
    free(mo.vir_sym);
 
242
    free(mo.occ_off);
 
243
    free(mo.vir_off);
 
244
    free(mo.qt_occ);
 
245
    free(mo.qt_vir);
 
246
  }
 
247
}
 
248
 
 
249
void exit_io(void)
 
250
{
 
251
  int i;
 
252
 
 
253
  psio_done();
 
254
  tstop(outfile);
 
255
  psi_stop(infile,outfile,psi_file_prefix);
 
256
}
 
257
 
 
258
}} /* End namespaces */