← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/psimrcc/mrcc_bwccsd.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/psimrcc/mrcc_bwccsd.cc
 
1
/***************************************************************************
 
2
 *  PSIMRCC : Copyright (C) 2007 by Francesco Evangelista and Andrew Simmonett
 
3
 *  frank@ccc.uga.edu   andysim@ccc.uga.edu
 
4
 *  A multireference coupled cluster code
 
5
 ***************************************************************************/
 
6
#include <libmoinfo/libmoinfo.h>
 
7
#include "mrcc.h"
 
8
#include "matrix.h"
 
9
#include "blas.h"
 
10
#include "debugging.h"
 
11
#include <libutil/libutil.h>
 
12
 
 
13
extern FILE* outfile;
 
14
 
 
15
namespace psi{ namespace psimrcc{
 
16
 
 
17
using namespace std;
 
18
 
 
19
void CCMRCC::update_amps_bwccsd_wrapper()
 
20
{
 
21
  ptr->update_amps_bwccsd();
 
22
}
 
23
 
 
24
void CCMRCC::update_amps_bwccsd()
 
25
{
 
26
  update_t1_amps_bwccsd();
 
27
  update_t2_amps_bwccsd();
 
28
 
 
29
  zero_internal_delta_amps();
 
30
 
 
31
  blas->solve("||Delta_t1||{u}  = t1_delta[o][v]{u} . t1_delta[o][v]{u}");
 
32
  blas->solve("||Delta_t1||{u} += t1_delta[O][V]{u} . t1_delta[O][V]{u}");
 
33
 
 
34
  blas->solve("||Delta_t2||{u}  = t2_delta[oo][vv]{u} . t2_delta[oo][vv]{u}");
 
35
  blas->solve("||Delta_t2||{u} += t2_delta[oO][vV]{u} . t2_delta[oO][vV]{u}");
 
36
  blas->solve("||Delta_t2||{u} += t2_delta[OO][VV]{u} . t2_delta[OO][VV]{u}");
 
37
 
 
38
  // Compute the T-AMPS difference
 
39
  delta_t1_amps=0.0;
 
40
  delta_t2_amps=0.0;
 
41
  for(int n=0;n<moinfo->get_ref_size("a");n++){
 
42
    delta_t1_amps+=blas->get_scalar("||Delta_t1||",moinfo->get_ref_number("a",n));
 
43
    delta_t2_amps+=blas->get_scalar("||Delta_t2||",moinfo->get_ref_number("a",n));
 
44
  }
 
45
  delta_t1_amps=pow(delta_t1_amps,0.5)/((double)moinfo->get_nrefs());
 
46
  delta_t2_amps=pow(delta_t2_amps,0.5)/((double)moinfo->get_nrefs());
 
47
}
 
48
 
 
49
void CCMRCC::update_t1_amps_bwccsd()
 
50
{
 
51
  blas->solve("d'1[o][v]{u}  = d1[o][v]{u}");
 
52
  blas->solve("d'1[O][V]{u}  = d1[O][V]{u}");
 
53
 
 
54
  for(int n=0;n<moinfo->get_ref_size("u");n++){
 
55
    double shift = current_energy-Heff[moinfo->get_ref_number("u",n)][moinfo->get_ref_number("u",n)];
 
56
    string n_str = to_string(moinfo->get_ref_number("u",n));
 
57
    blas->solve("d'1[o][v]{" + n_str + "} += " + to_string(shift));
 
58
    blas->solve("d'1[O][V]{" + n_str + "} += " + to_string(shift));
 
59
  }
 
60
 
 
61
  blas->solve("t1_delta[o][v]{u}  =   t1_eqns[o][v]{u} / d'1[o][v]{u} - t1[o][v]{u}");
 
62
  blas->solve("t1_delta[O][V]{u}  =   t1_eqns[O][V]{u} / d'1[O][V]{u} - t1[O][V]{u}");
 
63
 
 
64
  blas->solve("t1[o][v]{u} = t1_eqns[o][v]{u} / d'1[o][v]{u}");
 
65
  blas->solve("t1[O][V]{u} = t1_eqns[O][V]{u} / d'1[O][V]{u}");
 
66
 
 
67
  blas->solve("t1_norm{u}  = t1[o][v]{u} . t1[o][v]{u}");
 
68
  blas->solve("t1_norm{u} += t1[O][V]{u} . t1[O][V]{u}");
 
69
 
 
70
  zero_t1_internal_amps();
 
71
 
 
72
  DEBUGGING(3,
 
73
    blas->print("t1[o][v]{u}");
 
74
    blas->print("t1[O][V]{u}");
 
75
  );
 
76
}
 
77
 
 
78
void CCMRCC::update_t2_amps_bwccsd()
 
79
{
 
80
  blas->solve("d'2[oo][vv]{u}  = d2[oo][vv]{u}");
 
81
  blas->solve("d'2[oO][vV]{u}  = d2[oO][vV]{u}");
 
82
  blas->solve("d'2[OO][VV]{u}  = d2[OO][VV]{u}");
 
83
 
 
84
  for(int n=0;n<moinfo->get_ref_size("u");n++){
 
85
    string shift = to_string(current_energy-Heff[moinfo->get_ref_number("u",n)][moinfo->get_ref_number("u",n)]);
 
86
    string n_str = to_string(moinfo->get_ref_number("u",n));
 
87
    blas->solve("d'2[oo][vv]{" + n_str + "} += " + shift);
 
88
    blas->solve("d'2[oO][vV]{" + n_str + "} += " + shift);
 
89
    blas->solve("d'2[OO][VV]{" + n_str + "} += " + shift);
 
90
  }
 
91
 
 
92
 
 
93
 
 
94
  // (a) Compute eq. (20) of J. Chem. Phys. 110, 10275 (1999)
 
95
  // Comment : Look at eq. (21) of J. Chem. Phys. 110, 10275 (1999)
 
96
  blas->solve("t1_eqns[o][v]{u} += - d1[o][v]{u} * t1[o][v]{u}");
 
97
  blas->solve("t1_eqns[O][V]{u} += - d1[O][V]{u} * t1[O][V]{u}");
 
98
 
 
99
  // aaaa case
 
100
  // (b) Add PijPab (term from a) to the T2 equations 
 
101
  blas->solve("t2_eqns[oo][vv]{u} += #1324#   t1[o][v]{u} X t1_eqns[o][v]{u}");
 
102
  blas->solve("t2_eqns[oo][vv]{u} += #2314# - t1[o][v]{u} X t1_eqns[o][v]{u}");
 
103
  blas->solve("t2_eqns[oo][vv]{u} += #1423# - t1[o][v]{u} X t1_eqns[o][v]{u}");
 
104
  blas->solve("t2_eqns[oo][vv]{u} += #2413#   t1[o][v]{u} X t1_eqns[o][v]{u}");
 
105
  // (c) Subtract (term from c) from the T2 equations
 
106
  for(int n=0;n<moinfo->get_ref_size("u");n++){
 
107
    int ref = moinfo->get_ref_number("u",n);
 
108
    string neg_shift = to_string(-current_energy+Heff[ref][ref]);
 
109
    string shift     = to_string(current_energy-Heff[ref][ref]);
 
110
    string n_str     = to_string(ref);
 
111
    blas->solve("t2_eqns[oo][vv]{" + n_str + "} += #1324# " + neg_shift + "  t1[o][v]{" + n_str + "} X t1[o][v]{" + n_str + "}");
 
112
    blas->solve("t2_eqns[oo][vv]{" + n_str + "} += #2314# " + shift  + "  t1[o][v]{" + n_str + "} X t1[o][v]{" + n_str + "}");
 
113
  }
 
114
 
 
115
  // abab case
 
116
  // (b) Add PijPab (term from a) to the T2 equations 
 
117
  blas->solve("t2_eqns[oO][vV]{u} += #1324# t1[o][v]{u} X t1_eqns[O][V]{u}");
 
118
  blas->solve("t2_eqns[oO][vV]{u} += #2413# t1[O][V]{u} X t1_eqns[o][v]{u}");
 
119
  // (c) Subtract (term from c) from the T2 equations
 
120
  for(int n=0;n<moinfo->get_ref_size("u");n++){
 
121
    int ref = moinfo->get_ref_number("u",n);
 
122
    string n_str = to_string(ref);
 
123
    string neg_shift = to_string(-current_energy+Heff[ref][ref]);
 
124
    blas->solve("t2_eqns[oO][vV]{" + n_str + "} += #1324# " + neg_shift + "  t1[o][v]{" + n_str + "} X t1[O][V]{" + n_str + "}");
 
125
  }
 
126
 
 
127
  // bbbb case
 
128
  // (b) Add PijPab (term from a) to the T2 equations 
 
129
  blas->solve("t2_eqns[OO][VV]{u} += #1324#   t1[O][V]{u} X t1_eqns[O][V]{u}");
 
130
  blas->solve("t2_eqns[OO][VV]{u} += #2314# - t1[O][V]{u} X t1_eqns[O][V]{u}");
 
131
  blas->solve("t2_eqns[OO][VV]{u} += #1423# - t1[O][V]{u} X t1_eqns[O][V]{u}");
 
132
  blas->solve("t2_eqns[OO][VV]{u} += #2413#   t1[O][V]{u} X t1_eqns[O][V]{u}");
 
133
  // (c) Subtract (term from c) from the T2 equations
 
134
  for(int n=0;n<moinfo->get_nunique();n++){
 
135
    int ref = moinfo->get_ref_number("u",n);
 
136
    string neg_shift = to_string(-current_energy+Heff[ref][ref]);
 
137
    string shift     = to_string(current_energy-Heff[ref][ref]);
 
138
    string n_str     = to_string(ref);
 
139
    blas->solve("t2_eqns[OO][VV]{" + n_str + "} += #1324# " + neg_shift + "  t1[O][V]{" + n_str + "} X t1[O][V]{" + n_str + "}");
 
140
    blas->solve("t2_eqns[OO][VV]{" + n_str + "} += #2314# " + shift     + "  t1[O][V]{" + n_str + "} X t1[O][V]{" + n_str + "}");
 
141
  }
 
142
 
 
143
  blas->solve("t2_delta[oo][vv]{u} = t2_eqns[oo][vv]{u} / d'2[oo][vv]{u} - t2[oo][vv]{u}");
 
144
  blas->solve("t2_delta[oO][vV]{u} = t2_eqns[oO][vV]{u} / d'2[oO][vV]{u} - t2[oO][vV]{u}");
 
145
  blas->solve("t2_delta[OO][VV]{u} = t2_eqns[OO][VV]{u} / d'2[OO][VV]{u} - t2[OO][VV]{u}");
 
146
 
 
147
  blas->solve("t2[oo][vv]{u} = t2_eqns[oo][vv]{u} / d'2[oo][vv]{u}");
 
148
  blas->solve("t2[oO][vV]{u} = t2_eqns[oO][vV]{u} / d'2[oO][vV]{u}");
 
149
  blas->solve("t2[OO][VV]{u} = t2_eqns[OO][VV]{u} / d'2[OO][VV]{u}");
 
150
 
 
151
  DEBUGGING(3,
 
152
    blas->print("t2_eqns[oo][vv]{u}");
 
153
    blas->print("t2[oo][vv]{u}");
 
154
    blas->print("t2_eqns[oO][vV]{u}");
 
155
    blas->print("t2[oO][vV]{u}");
 
156
    blas->print("t2_eqns[OO][VV]{u}");
 
157
    blas->print("t2[OO][VV]{u}");
 
158
  );
 
159
}
 
160
 
 
161
}} /* End Namespaces */