← Back to branch summary

~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/cphf/cphf.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/cphf/cphf.cc
 
1
/*! \file
 
2
    \ingroup CPHF
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
/*
 
6
** CPHF: Program to solve the Coupled Perturbed Hartree-Fock equations
 
7
** for nuclear and electric field perturbations and to compute
 
8
** electric polarizabilities, harmonic vibrational frequencies, and IR
 
9
** intensities.
 
10
**
 
11
** Limitations of and future plans for this code:
 
12
**
 
13
** (1) Spin-restricted closed-shell Hartree-Fock (RHF) wave functions
 
14
** only.  Extension to ROHF and UHF cases is needed.
 
15
**
 
16
** (2) All two-electron integrals are held in core and used in the MO
 
17
** basis.  Out-of-core and AO-based algorithms are needed in mohess.c,
 
18
** cphf_X.c, and build_hessian.c in order to handle larger basis sets
 
19
** and to avoid the two-electron integral transformation.
 
20
**
 
21
** (3) Symmetry-blocking is used in most of the loops, but is not
 
22
** actually used to improve storage or computational order.  I've put
 
23
** this off because the nuclear perturbations (x, y, and z on each
 
24
** nucleus) are not yet symmetry-adapted.  Some effort in this area is
 
25
** needed.
 
26
**
 
27
** (4) Thermodynamic functions, including enthalpies, entropies, heat
 
28
** capacities, free energies, and partition functions can be computed,
 
29
** given the vibrational data computed in vibration.c.
 
30
**
 
31
** TDC, December 2001 and October 2002
 
32
*/
 
33
 
 
34
#include <cstdio>
 
35
#include <cstdlib>
 
36
#include <cstring>
 
37
#include <libipv1/ip_lib.h>
 
38
#include <libciomr/libciomr.h>
 
39
#include <libpsio/psio.h>
 
40
#include <libiwl/iwl.h>
 
41
#include <libchkpt/chkpt.h>
 
42
#include <libqt/qt.h>
 
43
#include <psifiles.h>
 
44
#include "globals.h"
 
45
 
 
46
namespace psi { namespace cphf {
 
47
 
 
48
void init_io(int argc, char *argv[]);
 
49
void exit_io(void);
 
50
void title(void);
 
51
void init_ioff(void);
 
52
 
 
53
void setup(void);
 
54
void cleanup(void);
 
55
 
 
56
void out_of_core(double ***, double ***, double ***, double **);
 
57
void sort_B(double ***, double **);
 
58
void sort_A(double **, double **);
 
59
void mohess(double **);
 
60
void cphf_X(double ***, double **, double **, double ***);
 
61
void cphf_F(double **, double ***);
 
62
void polarize(double ***);
 
63
void build_hessian(double ***, double ***, double **, double ***, double **);
 
64
void build_dipder(double ***);
 
65
void vibration(double **, double **);
 
66
void cphf_B(double ***, double **);
 
67
 
 
68
}} // namespace psi::cphf
 
69
 
 
70
int main(int argc, char *argv[])
 
71
{
 
72
  using namespace psi::cphf;
 
73
  int errcod = 0;
 
74
  int coord = 0;
 
75
  double ***F;
 
76
  double ***S;
 
77
  double ***B;
 
78
  double **A;
 
79
  double **Baijk;
 
80
  double **Aaibj;
 
81
  double ***UX;
 
82
  double ***UF;
 
83
  double **hessian;
 
84
  double **L;
 
85
  
 
86
  init_io(argc,argv);
 
87
  title();
 
88
  init_ioff();
 
89
 
 
90
  timer_init();
 
91
  timer_on("CPHF Main");
 
92
 
 
93
  print_lvl = 0;
 
94
  errcod = ip_data("PRINT", "%d", &(print_lvl), 0);
 
95
 
 
96
  setup();
 
97
 
 
98
  F = (double ***) malloc(natom*3 * sizeof(double **));
 
99
  S = (double ***) malloc(natom*3 * sizeof(double **));
 
100
  B = (double ***) malloc(natom*3 * sizeof(double **));
 
101
  for(coord=0; coord < natom*3; coord++) {
 
102
    F[coord] = block_matrix(nmo, nmo);
 
103
    S[coord] = block_matrix(nmo, nmo);
 
104
    B[coord] = block_matrix(nmo, nmo);
 
105
  }
 
106
 
 
107
  A = block_matrix(num_pq,num_pq);
 
108
 
 
109
  out_of_core(F, S, B, A);
 
110
 
 
111
  Baijk = block_matrix(natom*3, num_ai);
 
112
  
 
113
  sort_B(B, Baijk);
 
114
 
 
115
  for(coord=0; coord < natom*3; coord++) {
 
116
    free_block(B[coord]);
 
117
  }
 
118
  free(B);
 
119
 
 
120
  Aaibj = block_matrix(num_ai,num_ai);
 
121
  
 
122
  sort_A(A, Aaibj);
 
123
  
 
124
  mohess(Aaibj);
 
125
  
 
126
  UX = (double ***) malloc(natom*3 * sizeof(double **));
 
127
  for(coord=0; coord < natom*3; coord++) {
 
128
    UX[coord] = block_matrix(nmo,nmo);
 
129
  }
 
130
  
 
131
  cphf_X(S, Baijk, Aaibj, UX);
 
132
 
 
133
  if (X_only) { /* only compute cphf coefficients, then quit */
 
134
    free_block(A);
 
135
    free_block(Baijk);
 
136
    free_block(Aaibj);
 
137
 
 
138
    for(coord=0; coord < natom*3; coord++) {
 
139
      free_block(UX[coord]);
 
140
      free_block(F[coord]);
 
141
      free_block(S[coord]);
 
142
    }
 
143
    free(UX); free(F); free(S);
 
144
    timer_off("CPHF Main");
 
145
    timer_done();
 
146
    exit_io();
 
147
    exit(PSI_RETURN_SUCCESS);
 
148
  }
 
149
  
 
150
  UF = (double ***) malloc(3 * sizeof(double **));
 
151
  for(coord=0; coord < 3; coord++)
 
152
    UF[coord] = block_matrix(nmo,nmo);
 
153
 
 
154
  cphf_F(Aaibj, UF);
 
155
  
 
156
  polarize(UF);
 
157
  
 
158
  hessian = block_matrix(natom*3, natom*3);
 
159
  build_hessian(F, S, A, UX, hessian);
 
160
  
 
161
  dipder = block_matrix(3, natom*3);
 
162
  dipder_q = block_matrix(3, natom*3);
 
163
  build_dipder(UX);
 
164
  
 
165
  L = block_matrix(natom*3, natom*3);
 
166
  vibration(hessian,L);
 
167
 
 
168
  cphf_B(UX,L);
 
169
  
 
170
  cleanup(); 
 
171
 
 
172
  free_block(A);
 
173
  free_block(Baijk);
 
174
  free_block(Aaibj);
 
175
 
 
176
  for(coord=0; coord < natom*3; coord++) { 
 
177
    free_block(UX[coord]);
 
178
    free_block(F[coord]);
 
179
    free_block(S[coord]);
 
180
  }
 
181
  for(coord=0; coord < 3; coord++) { 
 
182
    free_block(UF[coord]);
 
183
  }
 
184
  free(UX); free(UF); free(F); free(S);
 
185
  free_block(hessian);
 
186
  free_block(dipder);
 
187
  free_block(L);
 
188
  
 
189
  timer_off("CPHF Main");
 
190
  timer_done();
 
191
 
 
192
  exit_io();
 
193
  exit(PSI_RETURN_SUCCESS);
 
194
}
 
195
 
 
196
namespace psi { namespace cphf {
 
197
 
 
198
extern "C" {
 
199
extern const char *gprgid(void);
 
200
};
 
201
 
 
202
void init_io(int argc, char *argv[])
 
203
{
 
204
  int i, num_unparsed;
 
205
  char *progid, *argv_unparsed[100];
 
206
 
 
207
  progid = (char *) malloc(strlen(gprgid())+2);
 
208
  sprintf(progid, ":%s",gprgid());
 
209
 
 
210
  X_only = 0;
 
211
  for (i=1, num_unparsed=0; i<argc; ++i) {
 
212
    if (!strcmp(argv[i],"--X_only")) /* only compute cphf coefficients, then quit */
 
213
      X_only = 1;
 
214
    else
 
215
      argv_unparsed[num_unparsed++] = argv[i];
 
216
  }
 
217
 
 
218
  psi_start(&infile,&outfile,&psi_file_prefix,num_unparsed,argv_unparsed,0);
 
219
  ip_cwk_add(progid);
 
220
  free(progid);
 
221
  tstart(outfile);
 
222
 
 
223
  psio_init(); psio_ipv1_config();
 
224
}
 
225
 
 
226
void title(void)
 
227
{
 
228
  fprintf(outfile, "\t\t\t**************************\n");
 
229
  fprintf(outfile, "\t\t\t*                        *\n");
 
230
  fprintf(outfile, "\t\t\t*          CPHF          *\n");
 
231
  fprintf(outfile, "\t\t\t*                        *\n");
 
232
  fprintf(outfile, "\t\t\t**************************\n");
 
233
}
 
234
 
 
235
void exit_io(void)
 
236
{
 
237
  int i;
 
238
  psio_done();
 
239
  tstop(outfile);
 
240
  psi_stop(infile,outfile,psi_file_prefix);
 
241
}
 
242
 
 
243
extern "C" const char *gprgid(void)
 
244
{
 
245
   const char *prgid = "CPHF";
 
246
 
 
247
   return(prgid);
 
248
}
 
249
 
 
250
void init_ioff(void)
 
251
{
 
252
  int i;
 
253
  ioff = (int *) malloc(IOFF_MAX * sizeof(int));
 
254
  ioff[0] = 0;
 
255
  for(i=1; i < IOFF_MAX; i++) {
 
256
      ioff[i] = ioff[i-1] + i;
 
257
  }
 
258
}
 
259
 
 
260
}} // namespace psi::cphf