← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/psimrcc/mrcc_Heff.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_Heff.cc
 
1
#include <libmoinfo/libmoinfo.h>
 
2
#include <liboptions/liboptions.h>
 
3
#include "mrcc.h"
 
4
#include "matrix.h"
 
5
#include "blas.h"
 
6
#include "debugging.h"
 
7
#include <libutil/libutil.h>
 
8
 
 
9
#include <libchkpt/chkpt.h>
 
10
 
 
11
extern FILE* outfile;
 
12
 
 
13
namespace psi{ namespace psimrcc{
 
14
 
 
15
using namespace std;
 
16
 
 
17
bool CCMRCC::build_diagonalize_Heff(int cycle, double time)
 
18
{
 
19
  total_time     = time;
 
20
  bool converged = false;
 
21
 
 
22
  build_Heff_diagonal();
 
23
  build_Heff_offdiagonal();
 
24
 
 
25
  if(cycle==0){
 
26
    current_energy=diagonalize_Heff(moinfo->get_root(),moinfo->get_nrefs(),Heff,eigenvector,true);
 
27
    old_energy=current_energy;
 
28
    print_eigensystem(moinfo->get_nrefs(),Heff,eigenvector);
 
29
  }
 
30
  if((cycle>0) || (cycle==-1)){
 
31
    // Compute the energy difference
 
32
    old_energy=current_energy;
 
33
    current_energy=diagonalize_Heff(moinfo->get_root(),moinfo->get_nrefs(),Heff,eigenvector,false);
 
34
 
 
35
    if(options_get_bool("PRINT_HEFF"))
 
36
      print_eigensystem(moinfo->get_nrefs(),Heff,eigenvector);
 
37
    DEBUGGING(3,
 
38
      print_eigensystem(moinfo->get_nrefs(),Heff,eigenvector);
 
39
    )
 
40
    double delta_energy = current_energy-old_energy;
 
41
    if(fabs(log10(fabs(delta_energy))) > options_get_int("CONVERGENCE"))
 
42
      converged = true;
 
43
 
 
44
 
 
45
///    TODO fix this code which is temporarly not working
 
46
//     if(options_get_int("DAMPING_FACTOR")>0){
 
47
//       if(fabs(delta_energy) < moinfo->get_no_damp_convergence()){
 
48
//         double damping_factor = moinfo->get_damping_factor();
 
49
//         damping_factor *= 0.95;
 
50
//         fprintf(outfile,"\n\t# Scaling damp factor to zero, damping_factor = %lf",moinfo->get_damping_factor());
 
51
//         moinfo->set_damping_factor(damping_factor);
 
52
//       }
 
53
//     }
 
54
  }
 
55
  print_mrccsd_energy(cycle);
 
56
  if(converged){
 
57
    print_eigensystem(moinfo->get_nrefs(),Heff,eigenvector);
 
58
    chkpt_wt_etot(current_energy);
 
59
  }
 
60
  return(converged);
 
61
}
 
62
 
 
63
void CCMRCC::build_Heff_diagonal()
 
64
{
 
65
  // Compute the diagonal elements of the effective Hamiltonian
 
66
  // using a simple UCCSD energy expression
 
67
  blas->solve("Eaa{u}   = t1[o][v]{u} . fock[o][v]{u}");
 
68
  blas->solve("Ebb{u}   = t1[O][V]{u} . fock[O][V]{u}");
 
69
 
 
70
  blas->solve("Eaaaa{u} = 1/4 tau[oo][vv]{u} . <[oo]:[vv]>");
 
71
  blas->solve("Eabab{u} =     tau[oO][vV]{u} . <[oo]|[vv]>");
 
72
  blas->solve("Ebbbb{u} = 1/4 tau[OO][VV]{u} . <[oo]:[vv]>");
 
73
 
 
74
  blas->solve("ECCSD{u}  = Eaa{u} + Ebb{u} + Eaaaa{u} + Eabab{u} + Ebbbb{u} + ERef{u}");
 
75
 
 
76
  for(int n=0;n<moinfo->get_nrefs();n++)
 
77
    Heff[n][n]=blas->get_scalar("ECCSD",moinfo->get_ref_number("a",n));
 
78
}
 
79
 
 
80
void CCMRCC::build_Heff_offdiagonal()
 
81
{
 
82
  for(int i=0;i<moinfo->get_ref_size("a");i++){
 
83
    int i_unique = moinfo->get_ref_number("a",i);
 
84
    // Find the off_diagonal elements for reference i
 
85
    // Loop over reference j (in a safe way)
 
86
    for(int j=0;j<moinfo->get_ref_size("a");j++){
 
87
      if(i!=j){
 
88
        vector<pair<int,int> >  alpha_internal_excitation = moinfo->get_alpha_internal_excitation(i,j);
 
89
        vector<pair<int,int> >   beta_internal_excitation = moinfo->get_beta_internal_excitation(i,j);
 
90
        double                   sign_internal_excitation = moinfo->get_sign_internal_excitation(i,j);
 
91
 
 
92
        double element = 0.0;
 
93
        if(i==i_unique){
 
94
          // Set alpha-alpha single excitations
 
95
          if((alpha_internal_excitation.size()==1)&&(beta_internal_excitation.size()==0))
 
96
            element=sign_internal_excitation * blas->get_MatTmp("t1_eqns[o][v]",i_unique,none)->get_two_address_element(
 
97
                                               alpha_internal_excitation[0].first,
 
98
                                               alpha_internal_excitation[0].second);
 
99
 
 
100
          // Set beta-beta single excitations
 
101
          if((alpha_internal_excitation.size()==0)&&(beta_internal_excitation.size()==1))
 
102
            element=sign_internal_excitation * blas->get_MatTmp("t1_eqns[O][V]",i_unique,none)->get_two_address_element(
 
103
                                               beta_internal_excitation[0].first,
 
104
                                               beta_internal_excitation[0].second);
 
105
 
 
106
          // Set (alpha,alpha)->(alpha,alpha) double excitations
 
107
          if((alpha_internal_excitation.size()==2)&&(beta_internal_excitation.size()==0))
 
108
            element=sign_internal_excitation * blas->get_MatTmp("t2_eqns[oo][vv]",i_unique,none)->get_four_address_element(
 
109
                                               alpha_internal_excitation[0].first,
 
110
                                               alpha_internal_excitation[1].first,
 
111
                                               alpha_internal_excitation[0].second,
 
112
                                               alpha_internal_excitation[1].second);
 
113
 
 
114
          // Set (alpha,beta)->(alpha,beta) double excitations
 
115
          if((alpha_internal_excitation.size()==1)&&(beta_internal_excitation.size()==1))
 
116
            element=sign_internal_excitation * blas->get_MatTmp("t2_eqns[oO][vV]",i_unique,none)->get_four_address_element(
 
117
                                               alpha_internal_excitation[0].first,
 
118
                                                beta_internal_excitation[0].first,
 
119
                                               alpha_internal_excitation[0].second,
 
120
                                                beta_internal_excitation[0].second);
 
121
 
 
122
          // Set (beta,beta)->(beta,beta) double excitations
 
123
          if((alpha_internal_excitation.size()==0)&&(beta_internal_excitation.size()==2))
 
124
            element=sign_internal_excitation * blas->get_MatTmp("t2_eqns[OO][VV]",i_unique,none)->get_four_address_element(
 
125
                                               beta_internal_excitation[0].first,
 
126
                                               beta_internal_excitation[1].first,
 
127
                                               beta_internal_excitation[0].second,
 
128
                                               beta_internal_excitation[1].second);
 
129
        }else{
 
130
          // Set alpha-alpha single excitations
 
131
          if((alpha_internal_excitation.size()==1)&&(beta_internal_excitation.size()==0))
 
132
            element=sign_internal_excitation * blas->get_MatTmp("t1_eqns[O][V]",i_unique,none)->get_two_address_element(
 
133
                                               alpha_internal_excitation[0].first,
 
134
                                               alpha_internal_excitation[0].second);
 
135
 
 
136
          // Set beta-beta single excitations
 
137
          if((alpha_internal_excitation.size()==0)&&(beta_internal_excitation.size()==1))
 
138
            element=sign_internal_excitation * blas->get_MatTmp("t1_eqns[o][v]",i_unique,none)->get_two_address_element(
 
139
                                               beta_internal_excitation[0].first,
 
140
                                               beta_internal_excitation[0].second);
 
141
 
 
142
          // Set (alpha,alpha)->(alpha,alpha) double excitations
 
143
          if((alpha_internal_excitation.size()==2)&&(beta_internal_excitation.size()==0))
 
144
            element=sign_internal_excitation * blas->get_MatTmp("t2_eqns[OO][VV]",i_unique,none)->get_four_address_element(
 
145
                                               alpha_internal_excitation[0].first,
 
146
                                               alpha_internal_excitation[1].first,
 
147
                                               alpha_internal_excitation[0].second,
 
148
                                               alpha_internal_excitation[1].second);
 
149
 
 
150
          // Set (alpha,beta)->(alpha,beta) double excitations
 
151
          if((alpha_internal_excitation.size()==1)&&(beta_internal_excitation.size()==1))
 
152
            element=sign_internal_excitation * blas->get_MatTmp("t2_eqns[oO][vV]",i_unique,none)->get_four_address_element(
 
153
                                                beta_internal_excitation[0].first,
 
154
                                               alpha_internal_excitation[0].first,
 
155
                                                beta_internal_excitation[0].second,
 
156
                                               alpha_internal_excitation[0].second);
 
157
 
 
158
          // Set (beta,beta)->(beta,beta) double excitations
 
159
          if((alpha_internal_excitation.size()==0)&&(beta_internal_excitation.size()==2))
 
160
            element=sign_internal_excitation * blas->get_MatTmp("t2_eqns[oo][vv]",i_unique,none)->get_four_address_element(
 
161
                                               beta_internal_excitation[0].first,
 
162
                                               beta_internal_excitation[1].first,
 
163
                                               beta_internal_excitation[0].second,
 
164
                                               beta_internal_excitation[1].second);
 
165
        }
 
166
        Heff[j][i]=element;
 
167
      }
 
168
    }
 
169
  }
 
170
}
 
171
 
 
172
}} /* End Namespaces */