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

« back to all changes in this revision

Viewing changes to src/bin/ccenergy/diagnostic.c

  • 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:
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <math.h>
4
 
#include <libciomr/libciomr.h>
5
 
#include <libdpd/dpd.h>
6
 
#define EXTERN
7
 
#include "globals.h"
8
 
 
9
 
double diagnostic(void)
10
 
{
11
 
  int h, nirreps, Gi, Ga;
12
 
  int i, a, I, A, row, col;
13
 
  int num_elec, num_elec_a, num_elec_b;
14
 
  int *occpi, *virtpi;
15
 
  int *occ_sym, *vir_sym;
16
 
  int *clsdpi, *uoccpi;
17
 
  int *openpi;
18
 
  double t1diag, t1diag_a, t1diag_b;
19
 
  dpdfile2 T1A, T1B;
20
 
 
21
 
  nirreps = moinfo.nirreps;
22
 
  clsdpi = moinfo.clsdpi; 
23
 
  uoccpi = moinfo.uoccpi;
24
 
  openpi = moinfo.openpi;
25
 
 
26
 
  if(params.ref != 2) {
27
 
    occpi = moinfo.occpi; virtpi = moinfo.virtpi;
28
 
    occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
29
 
  }
30
 
 
31
 
  /* Compute the number of electrons */
32
 
  for(h=0,num_elec_a=0,num_elec_b=0; h < nirreps; h++) {
33
 
    num_elec_a += clsdpi[h] + openpi[h];
34
 
    num_elec_b += clsdpi[h];
35
 
  }
36
 
  num_elec = num_elec_a + num_elec_b;
37
 
 
38
 
  if(params.ref == 0) { /** RHF **/
39
 
 
40
 
    dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
41
 
    t1diag = dpd_file2_dot_self(&T1A);
42
 
    dpd_file2_close(&T1A);
43
 
 
44
 
    t1diag /= num_elec;
45
 
    t1diag = sqrt(t1diag);
46
 
 
47
 
  }
48
 
  else if(params.ref == 1) { /** ROHF **/
49
 
 
50
 
    dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
51
 
    dpd_file2_mat_init(&T1A);
52
 
    dpd_file2_mat_rd(&T1A);
53
 
    dpd_file2_init(&T1B, CC_OEI, 0, 0, 1, "tia");
54
 
    dpd_file2_mat_init(&T1B);
55
 
    dpd_file2_mat_rd(&T1B);
56
 
 
57
 
    t1diag = 0.0;
58
 
    for(h=0; h < nirreps; h++) {
59
 
 
60
 
      for(i=0; i < (occpi[h] - openpi[h]); i++) {
61
 
        for(a=0; a < (virtpi[h] - openpi[h]); a++) {
62
 
 
63
 
          t1diag += (T1A.matrix[h][i][a] + T1B.matrix[h][i][a]) *
64
 
            (T1A.matrix[h][i][a] + T1B.matrix[h][i][a]);
65
 
        }
66
 
      }
67
 
 
68
 
      for(i=0; i < (occpi[h] - openpi[h]); i++) {
69
 
        for(a=0; a < openpi[h]; a++) {
70
 
          A = a + uoccpi[h];
71
 
 
72
 
          t1diag += 2 * T1B.matrix[h][i][A] * T1B.matrix[h][i][A];
73
 
        }
74
 
      }
75
 
 
76
 
      for(i=0; i < openpi[h]; i++) {
77
 
        I = i + clsdpi[h];
78
 
        for(a=0; a < (virtpi[h] - openpi[h]); a++) {
79
 
                  
80
 
          t1diag += 2 * T1A.matrix[h][I][a] * T1A.matrix[h][I][a];
81
 
        }
82
 
      }
83
 
    }
84
 
 
85
 
    t1diag /= num_elec;
86
 
    t1diag = sqrt(t1diag);
87
 
    t1diag *= 0.5;
88
 
 
89
 
    dpd_file2_mat_close(&T1A);
90
 
    dpd_file2_close(&T1A);
91
 
    dpd_file2_mat_close(&T1B);
92
 
    dpd_file2_close(&T1B);
93
 
 
94
 
  }
95
 
  else if(params.ref == 2) { /** UHF **/
96
 
 
97
 
    dpd_file2_init(&T1A, CC_OEI, 0, 0, 1, "tIA");
98
 
    dpd_file2_mat_init(&T1A);
99
 
    dpd_file2_mat_rd(&T1A);
100
 
    dpd_file2_init(&T1B, CC_OEI, 0, 2, 3, "tia");
101
 
    dpd_file2_mat_init(&T1B);
102
 
    dpd_file2_mat_rd(&T1B);
103
 
 
104
 
    t1diag_a = 0.0;
105
 
    t1diag_b = 0.0;
106
 
    for(h=0; h < nirreps; h++) {
107
 
 
108
 
      for(row=0; row < T1A.params->rowtot[h]; row++) 
109
 
        for(col=0; col < T1A.params->coltot[h]; col++) 
110
 
          t1diag_a += (T1A.matrix[h][row][col] * T1A.matrix[h][row][col]);
111
 
 
112
 
      for(row=0; row < T1B.params->rowtot[h]; row++) 
113
 
        for(col=0; col < T1B.params->coltot[h]; col++) 
114
 
          t1diag_b += (T1B.matrix[h][row][col] * T1B.matrix[h][row][col]);
115
 
    }
116
 
 
117
 
 
118
 
    t1diag = sqrt((t1diag_a + t1diag_b)/(num_elec_a + num_elec_b));
119
 
 
120
 
    dpd_file2_mat_close(&T1A);
121
 
    dpd_file2_mat_close(&T1B);
122
 
    dpd_file2_close(&T1A);
123
 
    dpd_file2_close(&T1B);
124
 
              
125
 
  }
126
 
 
127
 
  return t1diag; 
128
 
}