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

« back to all changes in this revision

Viewing changes to src/bin/mcscf/scf_energy.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:
 
1
#include <iostream>
 
2
#include <cmath>
 
3
 
 
4
#include <libmoinfo/libmoinfo.h>
 
5
#include <libciomr/libciomr.h>
 
6
 
 
7
#include "scf.h"
 
8
#include "memory_manager.h"
 
9
#include "sblock_matrix.h"
 
10
 
 
11
extern FILE* outfile;
 
12
 
 
13
namespace psi{ namespace mcscf{
 
14
 
 
15
double SCF::energy(int cycle,double old_energy)
 
16
{
 
17
  double electronic_energy = 0.0;
 
18
 
 
19
  T  = H;
 
20
  T += Fc;
 
21
  electronic_energy += dot(Dc,T);
 
22
 
 
23
  if(reference == rohf){
 
24
    T  = H;
 
25
    T.scale(0.5);
 
26
    T += Fo;
 
27
    electronic_energy += dot(Do,T);
 
28
  }
 
29
 
 
30
  total_energy = electronic_energy + moinfo_scf->get_nuclear_energy();
 
31
 
 
32
  if(reference == tcscf){
 
33
//     SBlockMatrix Dtc_sum("Dtc sum",nirreps,sopi,sopi);
 
34
// 
 
35
//     // Compute diagonal elements of H
 
36
//     for(int I = 0 ; I < nci; ++I){
 
37
//       Dtc_sum  = Dc;
 
38
//       Dtc_sum += Dtc[I];
 
39
//       construct_G(Dtc_sum,G,"PK");
 
40
//       T  = H;
 
41
//       T.scale(2.0);
 
42
//       T += G;
 
43
//       H_tcscf[I][I] = dot(Dtc_sum,T) + moinfo_scf->get_nuclear_energy();
 
44
//     }
 
45
// 
 
46
//     // Compute off-diagonal elements of H
 
47
//     for(int I = 0 ; I < nci; ++I){
 
48
//       for(int J = I + 1; J < nci; ++J){
 
49
//         construct_G(Dtc[I],G,"K");
 
50
//         H_tcscf[I][J] = H_tcscf[J][I] = - dot(Dtc[J],G);
 
51
//       }
 
52
//     }
 
53
 
 
54
//     fprintf(outfile,"\n  Hamiltonian");
 
55
//     for(int I = 0 ; I < nci; ++I){
 
56
//       fprintf(outfile,"\n    ");
 
57
//       for(int J = 0 ; J < nci; ++J)
 
58
//         fprintf(outfile," %11.8f ",H_tcscf[I][J]);
 
59
//     }
 
60
    
 
61
 
 
62
    // Compute the CI gradient
 
63
    norm_ci_grad = 0.0;
 
64
    for(int I = 0 ; I < nci; ++I){
 
65
      ci_grad[I] = 0.0;
 
66
      for(int J = 0; J < nci; ++J){
 
67
        ci_grad[I] += H_tcscf[I][J] * ci[J];
 
68
      }
 
69
      ci_grad[I] -= old_energy * ci[I];
 
70
      norm_ci_grad += fabs(ci_grad[I]);
 
71
    }
 
72
    
 
73
    double*  eigenvalues;
 
74
    double** eigenvectors;
 
75
    allocate1(double,eigenvalues,nci);
 
76
    allocate2(double,eigenvectors,nci,nci);
 
77
 
 
78
    sq_rsp(nci,nci,H_tcscf,eigenvalues,1,eigenvectors,1.0e-14);
 
79
  
 
80
    total_energy = eigenvalues[root];
 
81
 
 
82
    for(int I = 0 ; I < nci; ++I)
 
83
      ci[I] = eigenvectors[I][root]; 
 
84
 
 
85
    release1(eigenvalues);
 
86
    release2(eigenvectors);
 
87
  }
 
88
 
 
89
  fprintf(outfile,"\n  @SCF %4d  %20.12f  %20.12f",cycle,total_energy,total_energy-old_energy);
 
90
 
 
91
  if(reference == tcscf){
 
92
    fprintf(outfile,"\n    ci      = [");
 
93
    for(int I = 0 ; I < nci; ++I)
 
94
      fprintf(outfile,"%11.8f%s",ci[I], I != nci -1 ? "," : "");
 
95
    fprintf(outfile,"]");
 
96
 
 
97
    fprintf(outfile,"\n    ci_grad = [");
 
98
    for(int I = 0 ; I < nci; ++I)
 
99
      fprintf(outfile,"%11.8f%s",ci_grad[I], I != nci -1 ? "," : "");
 
100
    fprintf(outfile,"]");
 
101
  }
 
102
 
 
103
  fflush(outfile);
 
104
 
 
105
  return(total_energy); 
 
106
}
 
107
 
 
108
}} /* End Namespaces */
 
109