Skip to content

Commit

Permalink
Merge pull request #278 from m-a-d-n-e-s-s/pullreq-local-cis
Browse files Browse the repository at this point in the history 
Pullreq local cis
  • Loading branch information
@fbischoff
fbischoff committed Apr 20, 2018
2 parents 13839b4 + 1ca88b8 commit bb443fa
Show file tree
Hide file tree
Showing 4 changed files with 91 additions and 43 deletions.
14 changes: 7 additions & 7 deletions src/apps/chem/SCFOperators.cc
View file
Expand Up @@ -334,7 +334,7 @@ void Exchange::set_parameters(const vecfuncT& bra, const vecfuncT& ket,
CoulombOperatorPtr(world, lo, econv));
}

vecfuncT Exchange::operator()(const vecfuncT& vket) const {
vecfuncT Exchange::operator()(const vecfuncT& vket, const double& mul_tol) const {
const bool same = this->same();
int nocc = mo_bra.size();
int nf = vket.size();
Expand All @@ -352,11 +352,11 @@ vecfuncT Exchange::operator()(const vecfuncT& vket) const {
if (small_memory_) { // Smaller memory algorithm ... possible 2x saving using i-j sym
for(int i=0; i<nocc; ++i){
if(occ[i] > 0.0){
vecfuncT psif = mul_sparse(world, mo_bra[i], vket, tol); /// was vtol
vecfuncT psif = mul_sparse(world, mo_bra[i], vket, mul_tol); /// was vtol
truncate(world, psif);
psif = apply(world, *poisson.get(), psif);
truncate(world, psif);
psif = mul_sparse(world, mo_ket[i], psif, tol); /// was vtol
psif = mul_sparse(world, mo_ket[i], psif, mul_tol); /// was vtol
gaxpy(world, 1.0, Kf, occ[i], psif);
}
}
Expand All @@ -367,12 +367,12 @@ vecfuncT Exchange::operator()(const vecfuncT& vket) const {
if (same)
jtop = i + 1;
for (int j = 0; j < jtop; ++j) {
psif.push_back(mul_sparse(mo_bra[i], vket[j], tol, false));
psif.push_back(mul_sparse(mo_bra[i], vket[j], mul_tol, false));
}
}

world.gop.fence();
truncate(world, psif);
truncate(world, psif,tol);
psif = apply(world, *poisson.get(), psif);
truncate(world, psif, tol);
reconstruct(world, psif);
Expand All @@ -384,9 +384,9 @@ vecfuncT Exchange::operator()(const vecfuncT& vket) const {
if (same)
jtop = i + 1;
for (int j = 0; j < jtop; ++j, ++ij) {
psipsif[i * nf + j] = mul_sparse(psif[ij], mo_ket[i], false);
psipsif[i * nf + j] = mul_sparse(psif[ij], mo_ket[i],mul_tol ,false);
if (same && i != j) {
psipsif[j * nf + i] = mul_sparse(psif[ij], mo_ket[j], false);
psipsif[j * nf + i] = mul_sparse(psif[ij], mo_ket[j],mul_tol , false);
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion src/apps/chem/SCFOperators.h
View file
Expand Up @@ -359,7 +359,7 @@ class Exchange {
/// note that only one spin is used (either alpha or beta orbitals)
/// @param[in] vket the orbitals |i> that the operator is applied on
/// @return a vector of orbitals K| i>
vecfuncT operator()(const vecfuncT& vket) const;
vecfuncT operator()(const vecfuncT& vket,const double& mul_tol=0.0) const;

/// compute the matrix element <bra | K | ket>

Expand Down
115 changes: 81 additions & 34 deletions src/apps/chem/TDHF.cc
View file
Expand Up @@ -94,22 +94,22 @@ struct gauss_functor : public FunctionFunctorInterface<double,3> {

/// helper struct for computing the moments
struct xyz {
int direction;
xyz(int direction) : direction(direction) {}
double operator()(const coord_3d& r) const {
return r[direction];
}
int direction;
xyz(int direction) : direction(direction) {}
double operator()(const coord_3d& r) const {
return r[direction];
}
};

TDHF::TDHF(World &world, const CCParameters & param, const Nemo & nemo_):
world(world),
parameters(param),
nemo(nemo_),
g12(world,OT_G12,param),
mo_ket_(make_mo_ket(nemo_)),
mo_bra_(make_mo_bra(nemo_)),
Q(world,mo_bra_.get_vecfunction(),mo_ket_.get_vecfunction()),
msg(world) {
world(world),
parameters(param),
nemo(nemo_),
g12(world,OT_G12,param),
mo_ket_(make_mo_ket(nemo_)),
mo_bra_(make_mo_bra(nemo_)),
Q(world,mo_bra_.get_vecfunction(),mo_ket_.get_vecfunction()),
msg(world) {
msg.section("Initialize TDHF Class");
msg.debug = parameters.debug;
if (not param.no_compute_response) {
Expand All @@ -121,6 +121,23 @@ TDHF::TDHF(World &world, const CCParameters & param, const Nemo & nemo_):
if(world.rank()==0) std::cout << eps << "\n";
time.info();
}
if(nemo.get_calc()->param.localize){
Fock F(world, nemo.get_calc().get(), nemo.nuclear_correlation);
F_occ = F(get_active_mo_bra(),get_active_mo_ket());
for(size_t i=0;i<get_active_mo_ket().size();++i){
std::cout << std::scientific << std::setprecision(10);
if(world.rank()==0) std::cout << "F(" << i << "," << i << ")=" << F_occ(i,i) << "\n";
if(std::fabs(get_orbital_energy(i+param.freeze)-F_occ(i,i))>1.e-5){
if(world.rank()==0) std::cout << "eps(" << i << ")=" << get_orbital_energy(i) << " | diff=" << get_orbital_energy(i+param.freeze)-F_occ(i,i) << "\n";
}
}
}else{
F_occ = Tensor<double>(get_active_mo_bra().size(),get_active_mo_ket().size());
F_occ*=0.0;
for(size_t i=0;i<get_active_mo_ket().size();++i){
F_occ(i,i)=get_orbital_energy(i+param.freeze);
}
}
}

std::vector<CC_vecfunction> TDHF::read_vectors() const {
Expand Down Expand Up @@ -469,8 +486,9 @@ vecfuncT TDHF::get_tda_potential(const CC_vecfunction &x)const{
// Real Alpha density (with nuclear cusps)
const real_function_3d alpha_density=0.5*nemo.R_square*nemo_density;

// XC Potential
const XCOperator xc(world,xc_data, not nemo.get_calc()->param.spin_restricted,alpha_density,alpha_density);




// Apply Ground State Potential to x-states
vecfuncT Vpsi1;
Expand All @@ -492,20 +510,26 @@ vecfuncT TDHF::get_tda_potential(const CC_vecfunction &x)const{
CCTimer timeJ(world,"Jx");
const vecfuncT Jx=J(x.get_vecfunction());
timeJ.info(parameters.debug);
// Applied XC Potential
CCTimer timeXCx(world,"XCx");
real_function_3d xc_pot = xc.make_xc_potential();

// compute the asymptotic correction of exchange-correlation potential
if(nemo.do_ac()) {
double charge = double(nemo.molecule().total_nuclear_charge());
real_function_3d scaledJ = -1.0/charge*J.potential()*(1.0-hf_coeff);
xc_pot = nemo.get_ac().apply(xc_pot, scaledJ);
}

const vecfuncT XCx=mul(world, xc_pot, x.get_vecfunction());
// Ground State Potential applied to x, without exchange
Vpsi1 = Jx+XCx; // Nx removed
if(nemo.get_calc()->xc.is_dft()){
// XC Potential
const XCOperator xc(world,xc_data, not nemo.get_calc()->param.spin_restricted,alpha_density,alpha_density);

// Applied XC Potential
CCTimer timeXCx(world,"XCx");
real_function_3d xc_pot = xc.make_xc_potential();

// compute the asymptotic correction of exchange-correlation potential
if(nemo.do_ac()) {
double charge = double(nemo.molecule().total_nuclear_charge());
real_function_3d scaledJ = -1.0/charge*J.potential()*(1.0-hf_coeff);
xc_pot = nemo.get_ac().apply(xc_pot, scaledJ);
}

const vecfuncT XCx=mul(world, xc_pot, x.get_vecfunction());
// Ground State Potential applied to x, without exchange
Vpsi1 = Jx+XCx; // Nx removed
}else Vpsi1=Jx;
// add exchange if demanded
if(hf_coeff!=0.0){
CCTimer timeKx(world,"Kx");
Expand Down Expand Up @@ -538,13 +562,19 @@ vecfuncT TDHF::get_tda_potential(const CC_vecfunction &x)const{
Coulomb Jp(world);
real_function_3d density_pert=2.0*nemo.make_density(occ,active_bra,x.get_vecfunction());
Jp.potential()=Jp.compute_potential(density_pert);
// reconstruct the full perturbed density: do not truncate!
real_function_3d gamma=xc.apply_xc_kernel(density_pert);
vecfuncT XCp=mul(world,gamma,active_mo);
truncate(world,XCp);

vecfuncT XCp=zero_functions<double,3>(world,get_active_mo_ket().size());
if(nemo.get_calc()->xc.is_dft()){
// XC Potential
const XCOperator xc(world,xc_data, not nemo.get_calc()->param.spin_restricted,alpha_density,alpha_density);
// reconstruct the full perturbed density: do not truncate!
real_function_3d gamma=xc.apply_xc_kernel(density_pert);
vecfuncT XCp=mul(world,gamma,active_mo);
truncate(world,XCp);
}

if(parameters.tda_triplet){
if(gamma.norm2()!=0.0) MADNESS_EXCEPTION("Triplets only for CIS",1);
if(norm2(world,XCp)!=0.0) MADNESS_EXCEPTION("Triplets only for CIS",1);
Vpsi2 = XCp;
}
else Vpsi2 = Jp(active_mo)+XCp;
Expand Down Expand Up @@ -583,8 +613,25 @@ vecfuncT TDHF::get_tda_potential(const CC_vecfunction &x)const{
}
// whole tda potential
vecfuncT Vpsi = Vpsi1 + Q(Vpsi2);
// if the ground state is localized add the coupling terms
// canonical: -ei|xi> (part of greens function)
// local: -fik|xk> (fii|xi> part of greens functions, rest needs to be added)

if(nemo.get_calc()->param.localize){
const vecfuncT vx=x.get_vecfunction();
vecfuncT fock_coupling=madness::transform(world,vx,F_occ);
// subtract the diagonal terms
for(size_t i=0;i<fock_coupling.size();++i){
fock_coupling[i]=(fock_coupling[i]-F_occ(i,i)*vx[i]);
}
Vpsi-=fock_coupling;
}


truncate(world,Vpsi);



// debug output
if(parameters.debug or parameters.plot){
plot_plane(world,Vpsi,"Vpsi");
Expand Down Expand Up @@ -1135,7 +1182,7 @@ void TDHF::analyze(const std::vector<CC_vecfunction> &x) const {
<< std::fixed << std::setprecision(1) << root.excitation <<": "
<< std::fixed << std::setprecision(10) << root.omega << " Eh "
<< root.omega*constants::hartree_electron_volt_relationship << " eV\n";
std::cout << std::scientific;
std::cout << std::scientific;
print(" oscillator strength (length) ", osl);
print(" oscillator strength (velocity) ", osv);
// print out the most important amplitudes
Expand Down
3 changes: 2 additions & 1 deletion src/apps/chem/TDHF.h
View file
Expand Up @@ -209,7 +209,8 @@ class TDHF{

/// analyze the root: oscillator strength and contributions from occupied orbitals
void analyze(const std::vector<CC_vecfunction> &x) const;

/// Fock matrix for occupied orbitals
Tensor<double> F_occ;
/// The MPI Communicator
World& world;
/// The Parameters for the Calculations
Expand Down

0 comments on commit bb443fa

Please sign in to comment.