← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccresponse/analyze.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/ccresponse/analyze.cc
 
1
/*! \file
 
2
    \ingroup CCRESPONSE
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cmath>
 
7
#include <libdpd/dpd.h>
 
8
#include <libchkpt/chkpt.h>
 
9
#include <libqt/qt.h>
 
10
#include "globals.h"
 
11
 
 
12
double **Build_R(void);
 
13
double **Build_U(void);
 
14
 
 
15
void analyze(char *pert, char *cart, int irrep, double omega)
 
16
{
 
17
  int nirreps, h, i, j, a, b, ij, ab, u, v;
 
18
  int position, num_div, tot1, tot2, nvir, nso, nocc;
 
19
  double width, max, min, value, value2;
 
20
  double *amp_array;
 
21
  double **tmp, **T2trans, **T1trans;
 
22
  FILE *efile;
 
23
  dpdbuf4 I, T2, D;
 
24
  dpdfile2 T1;
 
25
  char lbl[32];
 
26
 
 
27
  nirreps = moinfo.nirreps;
 
28
  num_div = 500;
 
29
  max = 9;
 
30
  min = 0;
 
31
  width = (max-min) / (num_div);
 
32
 
 
33
 
 
34
  sprintf(lbl, "X_%s_%1s_%5.3f", pert, cart, omega);
 
35
  ffile(&efile, lbl, 1);
 
36
  amp_array = init_array(num_div);
 
37
 
 
38
  nvir = moinfo.virtpi[0];
 
39
  nocc = moinfo.occpi[0];
 
40
  nso = moinfo.nso;
 
41
 
 
42
  sprintf(lbl, "X_%s_%1s_IjAb (%5.3f)", pert, cart, omega);
 
43
  dpd_buf4_init(&T2, CC_LR, 0, 0, 5, 0, 5, 0, lbl);
 
44
  dpd_buf4_mat_irrep_init(&T2, 0);
 
45
  dpd_buf4_mat_irrep_rd(&T2, 0);
 
46
  T2trans = block_matrix(nocc*nocc, nso*nso);
 
47
  tmp = block_matrix(nvir, nso);
 
48
  tot1 = 0;
 
49
  tot2 = 0;
 
50
  for(ij=0; ij<T2.params->rowtot[0]; ij++) {
 
51
 
 
52
    C_DGEMM('n', 't', nvir, nso, nvir, 1.0, &(T2.matrix[0][ij][0]), nvir, 
 
53
            &(moinfo.C[0][0][0]), nvir, 0.0, &(tmp[0][0]), nso);
 
54
    C_DGEMM('n', 'n', nso, nso, nvir, 1.0, &(moinfo.C[0][0][0]), nvir,
 
55
            tmp[0], nso, 0.0, T2trans[ij], nso);
 
56
 
 
57
    for(ab=0; ab<nso*nso; ab++) {
 
58
      value = fabs(log10(fabs(T2trans[ij][ab])));
 
59
      tot2++;
 
60
      if ((value >= max) && (value <= (max+width))) {
 
61
        amp_array[num_div-1]++;
 
62
        tot1++;
 
63
      }
 
64
      else if ((value <= min) && (value >= (min-width))) {
 
65
        amp_array[0]++;
 
66
        tot1++;
 
67
      }
 
68
      else if ((value < max) && (value > min)) {
 
69
        position = floor((value-min)/width);
 
70
        amp_array[position]++;
 
71
        tot1++;
 
72
      }
 
73
    }
 
74
  }
 
75
  dpd_buf4_mat_irrep_close(&T2, 0);
 
76
  dpd_buf4_close(&T2);
 
77
  free_block(tmp);
 
78
  free_block(T2trans);
 
79
 
 
80
  value2 = 0;
 
81
  for (i = num_div-1; i >= 0; i--) {
 
82
    value = amp_array[i] / tot1;
 
83
    value2 += value;
 
84
    fprintf(efile, "%10.5lf %lf\n", -((i)*width)-min, value);
 
85
  }
 
86
  free(amp_array);
 
87
  fprintf(outfile, "Total number of converged T2 amplitudes = %d\n", tot2);
 
88
  fprintf(outfile, "Number of T2 amplitudes in analysis= %d\n", tot1);
 
89
  fclose(efile);
 
90
 
 
91
  num_div = 40;
 
92
  max = 2;
 
93
  min = -5;
 
94
  width = (max-min) / (num_div);
 
95
 
 
96
  sprintf(lbl, "X1_%s_%1s_%5.3f", pert, cart, omega);
 
97
  ffile(&efile, lbl, 1);
 
98
  amp_array = init_array(num_div);
 
99
 
 
100
  sprintf(lbl, "X_%s_%1s_IA (%5.3f)", pert, cart, omega);
 
101
  dpd_file2_init(&T1, CC_OEI, 0, 0, 1, lbl);
 
102
  dpd_file2_print(&T1, outfile);
 
103
  dpd_file2_mat_init(&T1);
 
104
  dpd_file2_mat_rd(&T1);
 
105
 
 
106
  /*
 
107
  T1trans = block_matrix(nocc, nso);
 
108
 
 
109
  C_DGEMM('n','t', nocc, nso, nvir, 1.0, &(T1.matrix[0][0][0]), nvir,
 
110
          &(moinfo.C[0][0][0]), nvir, 0.0, &(T1trans[0][0]), nso);
 
111
  */
 
112
 
 
113
  tot1 = tot2 = 0;
 
114
  for(i=0; i < nocc; i++) {
 
115
    for(a=0; a < nso; a++) {
 
116
      /*      value = fabs(log10(fabs(T1trans[i][a]))); */
 
117
      value = log10(fabs(T1.matrix[0][i][a]));
 
118
      tot2++;
 
119
      if ((value >= max) && (value <= (max+width))) {
 
120
        amp_array[num_div-1]++;
 
121
        tot1++;
 
122
      }
 
123
      else if ((value <= min) && (value >= (min-width))) {
 
124
        amp_array[0]++;
 
125
        tot1++;
 
126
      }
 
127
      else if ((value < max) && (value > min)) {
 
128
        position = floor((value-min)/width);
 
129
        amp_array[position]++;
 
130
        tot1++;
 
131
      }
 
132
    }
 
133
  }
 
134
  /*  free_block(T1trans); */
 
135
 
 
136
  dpd_file2_mat_close(&T1);
 
137
  dpd_file2_close(&T1);
 
138
 
 
139
  value2 = 0;
 
140
  for (i = num_div-1; i >= 0; i--) {
 
141
    value = amp_array[i] / tot1;
 
142
    value2 += value;
 
143
    fprintf(efile, "%10.5lf %lf\n", ((i)*width)-min, value);
 
144
  }
 
145
 
 
146
  free(amp_array);
 
147
  fclose(efile);
 
148
 
 
149
}
 
150