← Back to branch summary

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

« back to all changes in this revision

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