4
#include <libmoinfo/libmoinfo.h>
5
#include <libciomr/libciomr.h>
8
#include "memory_manager.h"
9
#include "sblock_matrix.h"
13
namespace psi{ namespace mcscf{
15
double SCF::energy(int cycle,double old_energy)
17
double electronic_energy = 0.0;
21
electronic_energy += dot(Dc,T);
23
if(reference == rohf){
27
electronic_energy += dot(Do,T);
30
total_energy = electronic_energy + moinfo_scf->get_nuclear_energy();
32
if(reference == tcscf){
33
// SBlockMatrix Dtc_sum("Dtc sum",nirreps,sopi,sopi);
35
// // Compute diagonal elements of H
36
// for(int I = 0 ; I < nci; ++I){
39
// construct_G(Dtc_sum,G,"PK");
43
// H_tcscf[I][I] = dot(Dtc_sum,T) + moinfo_scf->get_nuclear_energy();
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);
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]);
62
// Compute the CI gradient
64
for(int I = 0 ; I < nci; ++I){
66
for(int J = 0; J < nci; ++J){
67
ci_grad[I] += H_tcscf[I][J] * ci[J];
69
ci_grad[I] -= old_energy * ci[I];
70
norm_ci_grad += fabs(ci_grad[I]);
74
double** eigenvectors;
75
allocate1(double,eigenvalues,nci);
76
allocate2(double,eigenvectors,nci,nci);
78
sq_rsp(nci,nci,H_tcscf,eigenvalues,1,eigenvectors,1.0e-14);
80
total_energy = eigenvalues[root];
82
for(int I = 0 ; I < nci; ++I)
83
ci[I] = eigenvectors[I][root];
85
release1(eigenvalues);
86
release2(eigenvectors);
89
fprintf(outfile,"\n @SCF %4d %20.12f %20.12f",cycle,total_energy,total_energy-old_energy);
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 ? "," : "");
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,"]");
105
return(total_energy);
108
}} /* End Namespaces */