← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/oscillator_strength.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/ccdensity/oscillator_strength.c
 
1
#include <stdio.h>
 
2
#include <math.h>
 
3
#include <libciomr/libciomr.h>
 
4
#include <libiwl/iwl.h>
 
5
#include <libchkpt/chkpt.h>
 
6
#include <libdpd/dpd.h>
 
7
#include <libqt/qt.h>
 
8
#include <psifiles.h>
 
9
#define EXTERN
 
10
#include "globals.h"
 
11
#include <physconst.h>
 
12
 
 
13
#define IOFF_MAX 32641
 
14
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
 
15
 
 
16
void oscillator_strength(struct TD_Params *S)
 
17
{
 
18
  int nmo, nso, nao, noei, stat, i, I, h, j, nirreps, *ioff;
 
19
  int *order, *order_A, *order_B, *doccpi, natom, *clsdpi, *openpi, *orbspi;
 
20
  double **scf_pitzer, **scf_pitzer_A, **scf_pitzer_B;
 
21
  double **scf_qt, **scf_qt_A, **scf_qt_B, **X, **usotao;
 
22
  double *zvals, **geom;
 
23
  double *mu_x_ints, *mu_y_ints, *mu_z_ints;
 
24
  double **MUX_AO, **MUY_AO, **MUZ_AO;
 
25
  double **MUX_MO, **MUY_MO, **MUZ_MO;
 
26
  double **MUX_MO_A, **MUY_MO_A, **MUZ_MO_A;
 
27
  double **MUX_MO_B, **MUY_MO_B, **MUZ_MO_B;
 
28
  double lt_x, lt_y, lt_z;
 
29
  double rt_x, rt_y, rt_z;
 
30
  double ds_x, ds_y, ds_z;
 
31
  double f_x, f_y, f_z;
 
32
  double f;
 
33
 
 
34
  chkpt_init(PSIO_OPEN_OLD);
 
35
  if ((params.ref == 0) || (params.ref == 1))
 
36
    scf_pitzer = chkpt_rd_scf();
 
37
  else if(params.ref == 2) {
 
38
    scf_pitzer_A = chkpt_rd_alpha_scf();
 
39
    scf_pitzer_B = chkpt_rd_beta_scf();
 
40
  }
 
41
 
 
42
  nso = chkpt_rd_nso();
 
43
  nao = chkpt_rd_nao();
 
44
  nmo = chkpt_rd_nmo();
 
45
  usotao = chkpt_rd_usotao();
 
46
  clsdpi = chkpt_rd_clsdpi();
 
47
  openpi = chkpt_rd_openpi();
 
48
  orbspi = chkpt_rd_orbspi();
 
49
  nirreps = chkpt_rd_nirreps();
 
50
  natom = chkpt_rd_natom();
 
51
  zvals = chkpt_rd_zvals();
 
52
  geom = chkpt_rd_geom();
 
53
  chkpt_close();
 
54
 
 
55
  lt_x = lt_y = lt_z = 0.0;
 
56
  rt_x = rt_y = rt_z = 0.0;
 
57
  ds_x = ds_y = ds_z = 0.0;
 
58
  f_x = f_y = f_z = 0.0;
 
59
  f = 0;
 
60
 
 
61
  ioff = init_int_array(IOFF_MAX);
 
62
  ioff[0] = 0;
 
63
  for(i=1; i < IOFF_MAX; i++) ioff[i] = ioff[i-1] + i;
 
64
 
 
65
  doccpi = init_int_array(moinfo.nirreps);
 
66
  for(h=0; h < moinfo.nirreps; h++) 
 
67
    doccpi[h] = moinfo.frdocc[h] + moinfo.clsdpi[h];
 
68
 
 
69
  if((params.ref == 0) || (params.ref == 1)) {
 
70
    order = init_int_array(nmo);
 
71
 
 
72
    reorder_qt(doccpi, moinfo.openpi, moinfo.frdocc, moinfo.fruocc, 
 
73
               order, moinfo.orbspi, moinfo.nirreps);
 
74
 
 
75
    scf_qt = block_matrix(nmo, nmo);
 
76
    for(i=0; i < nmo; i++) {
 
77
      I = order[i];  /* Pitzer --> QT */
 
78
      for(j=0; j < nmo; j++) scf_qt[j][I] = scf_pitzer[j][i];
 
79
    }
 
80
 
 
81
    free(order);
 
82
    free_block(scf_pitzer);
 
83
  }
 
84
  else if(params.ref == 2) {
 
85
    order_A = init_int_array(nmo); 
 
86
    order_B = init_int_array(nmo); 
 
87
 
 
88
    reorder_qt_uhf(doccpi, moinfo.openpi, moinfo.frdocc, moinfo.fruocc, 
 
89
                  order_A,order_B, moinfo.orbspi, moinfo.nirreps);
 
90
 
 
91
    scf_qt_A = block_matrix(nmo, nmo);
 
92
    for(i=0; i < nmo; i++) {
 
93
      I = order_A[i];  /* Pitzer --> QT */
 
94
      for(j=0; j < nmo; j++) scf_qt_A[j][I] = scf_pitzer_A[j][i];
 
95
    }
 
96
 
 
97
    scf_qt_B = block_matrix(nmo, nmo);
 
98
    for(i=0; i < nmo; i++) {
 
99
      I = order_B[i];  /* Pitzer --> QT */
 
100
      for(j=0; j < nmo; j++) scf_qt_B[j][I] = scf_pitzer_B[j][i];
 
101
    }
 
102
 
 
103
    free(order_A);
 
104
    free(order_B);
 
105
    free_block(scf_pitzer_A);
 
106
    free_block(scf_pitzer_B);
 
107
  }
 
108
  free(doccpi);
 
109
 
 
110
  /*** Read in dipole moment integrals in the AO basis ***/
 
111
  noei = nao * (nao + 1)/2;
 
112
 
 
113
  mu_x_ints = init_array(noei);
 
114
  stat = iwl_rdone(PSIF_OEI,PSIF_AO_MX,mu_x_ints,noei,0,0,outfile);
 
115
  mu_y_ints = init_array(noei);
 
116
  stat = iwl_rdone(PSIF_OEI,PSIF_AO_MY,mu_y_ints,noei,0,0,outfile);
 
117
  mu_z_ints = init_array(noei);
 
118
  stat = iwl_rdone(PSIF_OEI,PSIF_AO_MZ,mu_z_ints,noei,0,0,outfile);
 
119
 
 
120
  MUX_AO = block_matrix(nao,nao);
 
121
  MUY_AO = block_matrix(nao,nao);
 
122
  MUZ_AO = block_matrix(nao,nao);
 
123
 
 
124
  for(i=0; i < nao; i++)
 
125
    for(j=0; j < nao; j++) {
 
126
      MUX_AO[i][j] = mu_x_ints[INDEX(i,j)];
 
127
      MUY_AO[i][j] = mu_y_ints[INDEX(i,j)];
 
128
      MUZ_AO[i][j] = mu_z_ints[INDEX(i,j)];
 
129
    }
 
130
 
 
131
/*   fprintf(outfile, "MUX_AOs\n"); */
 
132
/*   print_mat(MUX_AO, nao, nao, outfile); */
 
133
/*   fprintf(outfile, "MUY_AOs\n"); */
 
134
/*   print_mat(MUY_AO, nao, nao, outfile); */
 
135
/*   fprintf(outfile, "MUZ_AOs\n"); */
 
136
/*   print_mat(MUZ_AO, nao, nao, outfile); */
 
137
 
 
138
  MUX_MO = block_matrix(nso,nso);
 
139
  MUY_MO = block_matrix(nso,nso);
 
140
  MUZ_MO = block_matrix(nso,nso);
 
141
 
 
142
  /*** Transform the AO dipole integrals to the SO basis ***/
 
143
  X = block_matrix(nso,nao); /* just a temporary matrix */
 
144
 
 
145
  C_DGEMM('n','n',nso,nao,nao,1,&(usotao[0][0]),nao,&(MUX_AO[0][0]),nao,
 
146
          0,&(X[0][0]),nao);
 
147
  C_DGEMM('n','t',nso,nso,nao,1,&(X[0][0]),nao,&(usotao[0][0]),nao,
 
148
          0,&(MUX_MO[0][0]),nso);
 
149
 
 
150
  C_DGEMM('n','n',nso,nao,nao,1,&(usotao[0][0]),nao,&(MUY_AO[0][0]),nao,
 
151
          0,&(X[0][0]),nao);
 
152
  C_DGEMM('n','t',nso,nso,nao,1,&(X[0][0]),nao,&(usotao[0][0]),nao,
 
153
          0,&(MUY_MO[0][0]),nso);
 
154
 
 
155
  C_DGEMM('n','n',nso,nao,nao,1,&(usotao[0][0]),nao,&(MUZ_AO[0][0]),nao,
 
156
          0,&(X[0][0]),nao);
 
157
  C_DGEMM('n','t',nso,nso,nao,1,&(X[0][0]),nao,&(usotao[0][0]),nao,
 
158
          0,&(MUZ_MO[0][0]),nso);
 
159
 
 
160
  free(mu_x_ints); free(mu_y_ints); free(mu_z_ints);
 
161
  free_block(X);
 
162
  free_block(usotao);
 
163
  free_block(MUX_AO);
 
164
  free_block(MUY_AO);
 
165
  free_block(MUZ_AO);
 
166
 
 
167
  /*** Transform the SO dipole integrals to the MO basis ***/
 
168
 
 
169
  X = block_matrix(nmo,nmo); /* just a temporary matrix */
 
170
 
 
171
  if((params.ref == 0) || (params.ref == 1)) {
 
172
 
 
173
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt[0][0]),nmo,&(MUX_MO[0][0]),nmo,
 
174
            0,&(X[0][0]),nmo);
 
175
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt[0][0]),nmo,
 
176
            0,&(MUX_MO[0][0]),nmo);
 
177
 
 
178
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt[0][0]),nmo,&(MUY_MO[0][0]),nmo,
 
179
            0,&(X[0][0]),nmo);
 
180
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt[0][0]),nmo,
 
181
            0,&(MUY_MO[0][0]),nmo);
 
182
 
 
183
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt[0][0]),nmo,&(MUZ_MO[0][0]),nmo,
 
184
            0,&(X[0][0]),nmo);
 
185
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt[0][0]),nmo,
 
186
            0,&(MUZ_MO[0][0]),nmo);
 
187
 
 
188
    free_block(scf_qt);
 
189
  }
 
190
  else if((params.ref == 2)) {
 
191
 
 
192
    MUX_MO_A = block_matrix(nso,nso);
 
193
    MUY_MO_A = block_matrix(nso,nso);
 
194
    MUZ_MO_A = block_matrix(nso,nso);
 
195
    MUX_MO_B = block_matrix(nso,nso);
 
196
    MUY_MO_B = block_matrix(nso,nso);
 
197
    MUZ_MO_B = block_matrix(nso,nso);
 
198
 
 
199
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt_A[0][0]),nmo,&(MUX_MO[0][0]),nmo,
 
200
            0,&(X[0][0]),nmo);
 
201
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt_A[0][0]),nmo,
 
202
            0,&(MUX_MO_A[0][0]),nmo);
 
203
 
 
204
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt_B[0][0]),nmo,&(MUX_MO[0][0]),nmo,
 
205
            0,&(X[0][0]),nmo);
 
206
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt_B[0][0]),nmo,
 
207
            0,&(MUX_MO_B[0][0]),nmo);
 
208
 
 
209
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt_A[0][0]),nmo,&(MUY_MO[0][0]),nmo,
 
210
            0,&(X[0][0]),nmo);
 
211
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt_A[0][0]),nmo,
 
212
            0,&(MUY_MO_A[0][0]),nmo);
 
213
 
 
214
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt_B[0][0]),nmo,&(MUY_MO[0][0]),nmo,
 
215
            0,&(X[0][0]),nmo);
 
216
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt_B[0][0]),nmo,
 
217
            0,&(MUY_MO_B[0][0]),nmo);
 
218
 
 
219
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt_A[0][0]),nmo,&(MUZ_MO[0][0]),nmo,
 
220
            0,&(X[0][0]),nmo);
 
221
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt_A[0][0]),nmo,
 
222
            0,&(MUZ_MO_A[0][0]),nmo);
 
223
 
 
224
    C_DGEMM('t','n',nmo,nmo,nmo,1,&(scf_qt_B[0][0]),nmo,&(MUZ_MO[0][0]),nmo,
 
225
            0,&(X[0][0]),nmo);
 
226
    C_DGEMM('n','n',nmo,nmo,nmo,1,&(X[0][0]),nmo,&(scf_qt_B[0][0]),nmo,
 
227
            0,&(MUZ_MO_B[0][0]),nmo);
 
228
 
 
229
    free_block(MUX_MO);
 
230
    free_block(MUY_MO);
 
231
    free_block(MUZ_MO);
 
232
    free_block(scf_qt_A);
 
233
    free_block(scf_qt_B);
 
234
  }
 
235
 
 
236
  free_block(X);
 
237
 
 
238
/*   fprintf(outfile, "MUX_MOs\n"); */
 
239
/*   print_mat(MUX_MO, nmo, nmo, outfile); */
 
240
/*   fprintf(outfile, "MUY_MOs\n"); */
 
241
/*   print_mat(MUY_MO, nmo, nmo, outfile); */
 
242
/*   fprintf(outfile, "MUZ_MOs\n"); */
 
243
/*   print_mat(MUZ_MO, nmo, nmo, outfile); */
 
244
 
 
245
  fprintf(outfile,"\n\tOscillator Strength for %d%3s\n",S->root+1,
 
246
          moinfo.labels[S->irrep]);
 
247
  fprintf(outfile,"\t                              X    \t       Y    \t       Z\n");
 
248
 
 
249
  if((params.ref == 0) || (params.ref == 1)) {
 
250
 
 
251
    for(i=0; i < nmo; i++)
 
252
      for(j=0; j < nmo; j++) {
 
253
        lt_x += MUX_MO[i][j] * moinfo.ltd[i][j];
 
254
        lt_y += MUY_MO[i][j] * moinfo.ltd[i][j];
 
255
        lt_z += MUZ_MO[i][j] * moinfo.ltd[i][j];
 
256
      }
 
257
 
 
258
    for(i=0; i < nmo; i++)
 
259
      for(j=0; j < nmo; j++) {
 
260
        rt_x += MUX_MO[i][j] * moinfo.rtd[i][j];
 
261
        rt_y += MUY_MO[i][j] * moinfo.rtd[i][j];
 
262
        rt_z += MUZ_MO[i][j] * moinfo.rtd[i][j];
 
263
      }
 
264
  }
 
265
  else if(params.ref == 2) {
 
266
 
 
267
    for(i=0; i < nmo; i++)
 
268
      for(j=0; j < nmo; j++) {
 
269
        lt_x += MUX_MO_A[i][j] * moinfo.ltd_a[i][j];
 
270
        lt_y += MUY_MO_A[i][j] * moinfo.ltd_a[i][j];
 
271
        lt_z += MUZ_MO_A[i][j] * moinfo.ltd_a[i][j];
 
272
      }
 
273
 
 
274
    for(i=0; i < nmo; i++)
 
275
      for(j=0; j < nmo; j++) {
 
276
        rt_x += MUX_MO_A[i][j] * moinfo.rtd_a[i][j];
 
277
        rt_y += MUY_MO_A[i][j] * moinfo.rtd_a[i][j];
 
278
        rt_z += MUZ_MO_A[i][j] * moinfo.rtd_a[i][j];
 
279
      }
 
280
 
 
281
    for(i=0; i < nmo; i++)
 
282
      for(j=0; j < nmo; j++) {
 
283
        lt_x += MUX_MO_B[i][j] * moinfo.ltd_b[i][j];
 
284
        lt_y += MUY_MO_B[i][j] * moinfo.ltd_b[i][j];
 
285
        lt_z += MUZ_MO_B[i][j] * moinfo.ltd_b[i][j];
 
286
      }
 
287
 
 
288
    for(i=0; i < nmo; i++)
 
289
      for(j=0; j < nmo; j++) {
 
290
        rt_x += MUX_MO_B[i][j] * moinfo.rtd_b[i][j];
 
291
        rt_y += MUY_MO_B[i][j] * moinfo.rtd_b[i][j];
 
292
        rt_z += MUZ_MO_B[i][j] * moinfo.rtd_b[i][j];
 
293
      }
 
294
  }
 
295
 
 
296
  ds_x = lt_x * rt_x;
 
297
  ds_y = lt_y * rt_y;
 
298
  ds_z = lt_z * rt_z;
 
299
  
 
300
  f_x = (2*S->cceom_energy*ds_x)/3;
 
301
  f_y = (2*S->cceom_energy*ds_y)/3;
 
302
  f_z = (2*S->cceom_energy*ds_z)/3;
 
303
  
 
304
  f = f_x + f_y + f_z;
 
305
  S->OS = f;
 
306
  
 
307
  fprintf(outfile,"\t<0|mu_e|n>              %11.8lf \t %11.8lf \t %11.8lf\n",
 
308
          lt_x,lt_y,lt_z);
 
309
  fprintf(outfile,"\t<n|mu_e|0>              %11.8lf \t %11.8lf \t %11.8lf\n",
 
310
          rt_x,rt_y,rt_z);
 
311
  fprintf(outfile,"\tDipole Strength         %11.8lf \n",ds_x+ds_y+ds_z);
 
312
  fprintf(outfile,"\tOscillator Strength     %11.8lf \n",f_x+f_y+f_z);
 
313
  fflush(outfile);
 
314
 
 
315
  if((params.ref == 0) || (params.ref == 1)) {
 
316
    free_block(MUX_MO); 
 
317
    free_block(MUY_MO); 
 
318
    free_block(MUZ_MO);
 
319
  }
 
320
  else if(params.ref == 2) {
 
321
    free_block(MUX_MO_A); 
 
322
    free_block(MUY_MO_A); 
 
323
    free_block(MUZ_MO_A);
 
324
    free_block(MUX_MO_B); 
 
325
    free_block(MUY_MO_B); 
 
326
    free_block(MUZ_MO_B);
 
327
  }
 
328
 
 
329
  return;
 
330
}