Skip to content

Commit

Permalink
fixing expectation value of using the singles
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Aug 13, 2024
1 parent 598acd6 commit 072796c
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 38 deletions.
9 changes: 7 additions & 2 deletions src/madness/chem/BSHApply.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ class BSHApply {
double levelshift=0.0;
double lo=1.e-6;
double bshtol=1.e-5;
bool printme=false;
bool printme=true;
bool destroy_Vpsi=false;
Function<double,NDIM> metric;
return_value ret_value=residual; // return the new orbitals/functions or the residuals
Expand Down Expand Up @@ -64,6 +64,7 @@ class BSHApply {
std::vector < std::shared_ptr<SeparatedConvolution<double,NDIM> > > ops(psi.size());
for (int i=0; i<eps.dim(0); ++i) {
T e_i= (eps.ndim()==2) ? eps(i,i) : eps(i);
if (printme) print("orbital energy for the BSH operator",e_i);
ops[i]=std::shared_ptr<SeparatedConvolution<double,NDIM> >(
BSHOperatorPtr<NDIM>(world, sqrt(-2.0*eps_in_green(e_i)), lo, bshtol));
ops[i]->destructive()=true;
Expand Down Expand Up @@ -130,7 +131,7 @@ class BSHApply {
const Tensor<T> fock1) const {

// check dimensions
bool consistent=(psi.size()==size_t(fock1.dim(0)));
bool consistent=(psi.size()==size_t(fock1.dim(0)));
if ((fock1.ndim()==2) and not (psi.size()==size_t(fock1.dim(1)))) consistent=false;

if (not consistent) {
Expand All @@ -150,6 +151,10 @@ class BSHApply {
for (int i=0; i<fock.dim(0); ++i) {
fock(i,i)-=eps_in_green(fock(i,i));
}
if (printme) {
print("coupling fock matrix");
print(fock);
}
return transform(world, psi, fock);

} else {
Expand Down
57 changes: 30 additions & 27 deletions src/madness/chem/CC2.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,6 +239,7 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
if (ctype == CT_LRCC2) omega = singles.omega;
else if (ctype == CT_LRCCS) omega = singles.omega;
else if (ctype == CT_ADC2) omega = singles.omega;
print("omega 1" ,omega);

// consistency check
switch (ctype) {
Expand Down Expand Up @@ -276,47 +277,35 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
else if (ctype == CT_LRCCS) V = CCPotentials::get_CCS_potential_ex(world,singles,false, info);
// else if (ctype == CT_ADC2) V = CCOPS.get_ADC2_singles_potential(world, gs_doubles, singles, ex_doubles, info);
else MADNESS_EXCEPTION("iterate singles: unknown type", 1);

// add local coupling
V-=compute_local_coupling(singles.get_vecfunction(),info);
truncate(world, V);
time_V.info(true, norm2(world, V));

if (ctype == CT_LRCCS or ctype == CT_LRCC2 or ctype == CT_ADC2) {
// update excitation energy
if (ctype==CT_LRCC2 or ctype==CT_LRCCS or ctype==CT_ADC2) {
old_omega=omega;
omega = singles.omega; // computed with the potential
omega = CCPotentials::compute_cis_expectation_value(world, singles, V, true, info);
singles.omega = omega;
}

scale(world, V, -2.0); // moved to BSHApply
truncate(world, V);

// BSHApply<double,3> bsh_apply(world);
// bsh_apply.ret_value=BSHApply<double,3>::update; // return the new singles functions, not the residual
// bsh_apply.metric=info.R_square;
//
// // coupling between singles involves the active fock matrix shifted by the excitation energy
// auto nfreeze=info.parameters.freeze();
// Tensor<double> fock=info.fock(Slice(nfreeze,-1),Slice(nfreeze,-1));
// for (int i=0; i<fock.dim(0); ++i) fock(i,i)+=omega;
// if (world.rank()==0 and info.parameters.debug()) {
// print("active Fock matrix shifted by omega=",omega);
// print(fock);
// }

if (world.rank()==0 and info.parameters.debug())
print("omega entering the update in the singles" ,omega);

// make bsh operators
CCTimer time_makebsh(world, "Make G-Operators");
scale(world, V, -2.0); // moved to BSHApply
std::vector<std::shared_ptr<SeparatedConvolution<double, 3> > > G(singles.size());
for (size_t i = 0; i < G.size(); i++) {
const double bsh_eps = info.orbital_energies[i + info.parameters.freeze()] + omega;
G[i] = std::shared_ptr<SeparatedConvolution<double, 3> >(
BSHOperatorPtr3D(world, sqrt(-2.0 * bsh_eps), info.parameters.lo(), info.parameters.thresh_bsh_3D()));
}
world.gop.fence();
time_makebsh.info();
//
// // apply bsh operators

// apply bsh operators
CCTimer time_applyG(world, "Apply G-Operators");
// auto [GV, energy_update] = bsh_apply(singles.get_vecfunction(), fock, V);
vector_real_function_3d GV = apply<SeparatedConvolution<double, 3>, double, 3>(world, G, V);
// world.gop.fence();
// auto GV=res-singles.get_vecfunction();
world.gop.fence();
time_applyG.info();

// apply Q-Projector to result
Expand Down Expand Up @@ -513,9 +502,23 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface {
for (auto& tmp_pair : vpairs) pairs.insert(tmp_pair.i, tmp_pair.j, tmp_pair);
auto ccpair2function = [](const CCPair& a) {return a.function();};
return compute_local_coupling(pairs.convert<real_function_6d>(pairs,ccpair2function), info);

};

/// compute the coupling of singles function if orbitals are localized

/// @return the coupling terms c_i = -\sum_(j\neq i) f_ij |\phi_j> (for whatever phi is)
static std::vector<real_function_3d> compute_local_coupling(const std::vector<real_function_3d>& singles,
const Info& info) {

MADNESS_CHECK_THROW(singles.size()>0,"compute_local_coupling: singles vector is empty");
World& world=singles.front().world();
auto active=Slice(info.parameters.freeze(),-1);
Tensor<double> Fact=info.fock(active,active);
for (int i=0; i<Fact.dim(0); ++i) Fact(i,i)=0.0;
vector_real_function_3d fock_coupling = madness::transform(world, singles, Fact);
return fock_coupling;
}

/// add the coupling terms for local MP2

/// \sum_{k\neq i} f_ki |u_kj> + \sum_{l\neq j} f_lj |u_il>
Expand Down
14 changes: 8 additions & 6 deletions src/madness/chem/CCPotentials.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -486,10 +486,16 @@ CCPotentials::compute_kinetic_energy(World& world, const vector_real_function_3d
return kinetic;
}


double
CCPotentials::compute_cis_expectation_value(World& world, const CC_vecfunction& x,
const vector_real_function_3d& V, const bool print, const Info& info)
{
// following eq. (34) of the CIS paper Kottmann et al, PCCP, 17, 31453, (2015)
// doi: https://doi.org/10.1039/C5CP00345H
// the expectation value of the CIS wave function is computed by projecting the
// CIS wave function onto eq. (22)
// the potential V must contain the coupling term when using localized orbitals
const vector_real_function_3d xbra = info.R_square*(x.get_vecfunction());
const vector_real_function_3d xket = x.get_vecfunction();
const double kinetic = compute_kinetic_energy(world, xbra, xket);
Expand Down Expand Up @@ -2879,7 +2885,7 @@ CCPotentials::get_CC2_singles_potential_gs(World& world, const CC_vecfunction& s
}

madness::vector_real_function_3d
CCPotentials::get_CCS_potential_ex(World& world, CC_vecfunction& x, const bool print, Info& info) {
CCPotentials::get_CCS_potential_ex(World& world, const CC_vecfunction& x, const bool print, Info& info) {
if (x.type != RESPONSE) error("get_CCS_response_potential: Wrong type of input singles");

Pairs<CCPair> empty_doubles;
Expand All @@ -2895,14 +2901,12 @@ CCPotentials::get_CCS_potential_ex(World& world, CC_vecfunction& x, const bool p
info.intermediate_potentials.insert(copy(world, potential), x, POT_singles_);
vector_real_function_3d result = add(world, fock_residue, potential);
truncate(world, result);
const double omega = compute_cis_expectation_value(world, x, result, print, info);
x.omega = omega;
return result;
}

madness::vector_real_function_3d
CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& gs_singles,
const Pairs<CCPair>& gs_doubles, CC_vecfunction& ex_singles,
const Pairs<CCPair>& gs_doubles, const CC_vecfunction& ex_singles,
const Pairs<CCPair>& response_doubles, Info& info)
{
MADNESS_ASSERT(gs_singles.type == PARTICLE);
Expand Down Expand Up @@ -2949,8 +2953,6 @@ CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& g
info.intermediate_potentials.insert(copy(world, potential), ex_singles, POT_singles_);
vector_real_function_3d result = add(world, fock_residue, potential);
truncate(world, result);
const double omega = compute_cis_expectation_value(world, ex_singles, result, true, info);
ex_singles.omega = omega;
return result;
}

Expand Down
6 changes: 3 additions & 3 deletions src/madness/chem/CCPotentials.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ class CCPotentials {
static double
compute_kinetic_energy(World& world, const vector_real_function_3d& xbra, const vector_real_function_3d& xket);

/// returns \f$ <x|T|x> + <x|V|x> \f$
/// compute the expectation value excitation energy using the CIS/CCS/CC2 singles
static double
compute_cis_expectation_value(World& world, const CC_vecfunction& x,
const vector_real_function_3d& V, const bool print, const Info& info);
Expand Down Expand Up @@ -681,13 +681,13 @@ class CCPotentials {
/// the V part is stored in the intermediate_potentials structure
/// the expectation value is calculated and updated
static vector_real_function_3d
get_CCS_potential_ex(World& world, CC_vecfunction& x, const bool print, Info& info);
get_CCS_potential_ex(World& world, const CC_vecfunction& x, const bool print, Info& info);

/// Calculates the CC2 singles potential for the Excited state: result = Fock_residue + V
/// the V part is stored in the intermediate_potentials structure
static vector_real_function_3d
get_CC2_singles_potential_ex(World& world, const CC_vecfunction& gs_singles,
const Pairs<CCPair>& gs_doubles, CC_vecfunction& ex_singles,
const Pairs<CCPair>& gs_doubles, const CC_vecfunction& ex_singles,
const Pairs<CCPair>& response_doubles, Info& info);

/// Calculates the CC2 singles potential for the Excited state: result = Fock_residue + V
Expand Down

0 comments on commit 072796c

Please sign in to comment.