TDA.h

Go to the documentation of this file.

/*
 * TDA.h
 *
 *  Created on: Jul 14, 2014
 *      Author: kottmanj
 */

#ifndef TDA_H_
#define TDA_H_

//#include<examples/dft_solver.h>
#include <chem/projector.h>
//#include<examples/nonlinsol.h> not used anymore
#include<chem/SCF.h>
#include <madness/mra/operator.h>
#include <madness/mra/mra.h>
#include <madness/mra/vmra.h>
#include <madness/mra/lbdeux.h>
#include <madness/misc/ran.h>
#include <chem/TDA_XC.h>

#include <chem/TDA_exops.h>

// Kain solver
#include<examples/nonlinsol.h>

// std::sort
#include <algorithm>

namespace madness {

typedef std::vector<Function<double,3> > vecfuncT;
typedef SeparatedConvolution<double,3> operatorT;
typedef std::shared_ptr<operatorT> poperatorT;
typedef std::shared_ptr< FunctionFunctorInterface<double,3> > functorT;

struct TDA_TIMER{
        TDA_TIMER(World &world,std::string msg) : world(world),start_wall(wall_time()),start_cpu(cpu_time()),operation(msg){}
        World & world;
        const double start_wall;
        const double start_cpu;
        std::string operation;
        double end_wall;
        double end_cpu;
        void update_time(){
                end_wall = wall_time()-start_wall;
                end_cpu = cpu_time()-start_cpu;
        }
public:
        void info(bool debug = true){
                if(debug==true){
                        update_time();
                        if(world.rank()==0) std::cout<< std::setw(40) << operation << " : " << std::scientific << std::setprecision(1)
                        << end_wall << " (wall) "<< end_cpu << " (cpu)" << std::endl;
                }
        }

};

struct xfunction{
        xfunction(World &world) :world(world),omega(0.00001),converged(false),number(100),iterations(0),kain(false),f_length(999),f_velocity(999) {}
        xfunction(World& world, const vecfuncT& x1) : world(world), x(x1),omega(0.00001),converged(false),number(100),iterations(0),kain(true),f_length(999),f_velocity(999) {}
        xfunction(const xfunction &other) : world(other.world),x(other.x),Vx(other.Vx),omega(other.omega),expectation_value(other.expectation_value),error(other.error),
                        delta(other.delta),converged(other.converged),number(other.number),iterations(other.iterations),kain(other.kain),f_length(other.f_length),f_velocity(other.f_velocity){}

        World & world;
        vecfuncT x;
        vecfuncT Vx;
        double omega;
        std::vector<double> expectation_value;
        std::vector<double> error;
        std::vector<double> delta;
        bool converged;
        size_t number;
        size_t iterations;
        bool kain;
        vecfuncT current_residuals;

        double f_length;
        double f_velocity;

        xfunction& operator=(const xfunction &other){
                x=other.x;
                Vx=other.Vx;
                omega = other.omega;
                expectation_value = other.expectation_value;
                error=other.error;
                delta=other.delta;
                converged=other.converged;
                number=other.number;
                iterations=other.iterations;
                kain=other.kain;
                f_length = other.f_length;
                f_velocity = other.f_velocity;

                return *this;
        }

        // Operators needed by the KAIN solver
        xfunction operator-=(const xfunction& b) {
                x=sub(world,x,b.x);
                return *this;
        }

        xfunction operator-(const xfunction &b)const {
                return xfunction(world,sub(world,x,b.x));
        }

        xfunction operator+=(const xfunction& b) { // Operator+= necessary
                x=add(world,x,b.x);
                return *this;
        }

        xfunction operator*(double a) const { // Scale by a constant necessary

                PROFILE_BLOCK(Vscale);
                xfunction result(*this);
                for (unsigned int i=0; i<x.size(); ++i) {
                        result.x[i]=mul(a,x[i],false);
                }
                world.gop.fence();

                //              scale(world, x, a);
                return result;
        }

        // finally an operator that the result can be sorted after the energy
        // sorting of xfunctions should not happen during the iterations
        // therefore a warning is installed
        bool operator<=(const xfunction &b)const{return omega<=b.omega;std::cout << "WARNING XFUNCTIONS ARE SORTED" << std::endl;}
        bool operator< (const xfunction &b)const{return omega<b.omega; std::cout << "WARNING XFUNCTIONS ARE SORTED" << std::endl;}



};

struct allocator{
        World& world;
        const int noct;

        allocator(World& world, const int nnoct) : world(world), noct(nnoct) {}

        xfunction operator()(){
                return xfunction(world,zero_functions<double,3>(world,noct));
        }
        allocator operator=(const allocator &other){
                allocator tmp(world,other.noct);
                return tmp;
        }
};

static double inner(const xfunction &a, const xfunction &b) {
        if (a.x.size()!=b.x.size()) MADNESS_EXCEPTION("ERROR :Inner product of two xfunction structures: Different sizes in x-vectors",1);
        if (a.x.size()==0) return 0.0;
        return madness::inner(a.x[0].world(),a.x,b.x).sum();
}

// TYPEDEFS
typedef std::vector<xfunction> xfunctionsT;
typedef XNonlinearSolver<xfunction,double,allocator> solverT;

struct kain_solver_helper_struct{
        kain_solver_helper_struct(){}
private:
        size_t noct;
        bool is_used;
        std::vector<solverT> solver;
public:

        void sanity_check(const xfunctionsT &xfunctions)const{
                if(xfunctions.size()!=solver.size()){
                        std::cout << "KAIN SOLVER SANITY CHECK FAILED: Unequal sizes " << xfunctions.size() << " xfunctions and " << solver.size() << " solvers" << std::endl;
                        MADNESS_EXCEPTION("Error in kain solver helper struct",1);
                }
                if(xfunctions[0].x.size()!=noct){
                        MADNESS_EXCEPTION("Error in kain solver helper struct, unequal sizes in occupied orbitals (check freeze_ keyword)",1);
                }
        }

        void initialize(World &world,const size_t excitations,const size_t nnoct,const bool kain){
                noct=nnoct;
                is_used = kain;
                if(kain){
                        for(size_t i=0;i<excitations; i++){
                                solverT onesolver(allocator(world,noct));
                                onesolver.set_maxsub(3);
                                solver.push_back(onesolver);
                        }
                }
        }
        void update(xfunctionsT &xfunctions){
                if(not is_used) MADNESS_EXCEPTION("Kain update was requested, but kain should not be used",1);
                if(is_used){
                        for(size_t k=0;k<xfunctions.size();k++){
                                // we need ot push this back if we update or not, because of the subspace transformation
                                xfunction tmp = solver[k].update(xfunction(xfunctions[k].world,xfunctions[k].x),xfunction(xfunctions[k].world,xfunctions[k].current_residuals));
                                if(xfunctions[k].kain){
                                        xfunctions[k].x = tmp.x;
                                }else{
                                        std::cout << "(convergence too low) forcing full step for  " << k << std::endl;
                                        xfunctions[k].x = sub(xfunctions[k].world,xfunctions[k].x,xfunctions[k].current_residuals);
                                }
                        }
                }
        }
        void transform_subspace(World &world,const madness::Tensor<double> U){
                if(solver[0].get_ulist().empty()){ std::cout<<"kain: empty ulist, no transformation in "; return;}
                if(solver[0].get_rlist().empty()){ std::cout<<"kain: empty rlist, no transformation in "; return;}
                for(int k=0;k<solver[0].get_ulist().size();k++){
                        std::vector<vecfuncT> all_x;
                        std::vector<vecfuncT> all_r;
                        for(size_t i=0;i<solver.size();i++){
                                all_x.push_back(solver[i].get_ulist()[k].x);
                                all_r.push_back(solver[i].get_rlist()[k].x);
                        }
                        std::vector<vecfuncT> new_x = transform_vecfunctions(world,all_x,U);
                        std::vector<vecfuncT> new_r = transform_vecfunctions(world,all_r,U);
                        for(size_t i=0;i<solver.size();i++){
                                solver[i].get_ulist()[k].x = new_x[i];
                                solver[i].get_rlist()[k].x = new_r[i];
                        }
                }
        }
        void reduce_subspace(const size_t i){
                solver.erase(solver.begin()+i);
        }
        void erase_subspace(const size_t i){
                if(solver.size()<i) MADNESS_EXCEPTION("Tried to delete subspace of nonexisting solver",1);
                solver[i].get_ulist().clear();
                solver[i].get_rlist().clear();
        }
        void increase_subspace(World &world,const xfunctionsT &xfunctions){
                // if the subspace is increased (means more solvers are added) the subspaces of the existing solvers must be erased
                // else there will be problems with the next subspace transformation
                // the simplest solution is to add cimpletely new solvers for all xfunctions

                // clean up: erase all old solvers
                solver.clear();
                // add new solvers
                solverT onesolver(allocator(world,noct));
                onesolver.set_maxsub(3);
                for(size_t j=0;j<xfunctions.size();j++){
                        solver.push_back(onesolver);
                }
                sanity_check(xfunctions);
        }
        std::vector<vecfuncT> transform_vecfunctions(World &world,const std::vector<vecfuncT> &xfunctions,const madness::Tensor<double> U)const{

                std::vector<vecfuncT> new_xfunctions(xfunctions.size());
                for (std::size_t i = 0; i < xfunctions.size(); i++) {
                        new_xfunctions[i] = zero_functions<double, 3>(world,
                                        xfunctions[i].size());
                        compress(world, new_xfunctions[i]);
                        compress(world, xfunctions[i]);
                }

                for (size_t i = 0; i < xfunctions.size(); i++) {
                        for (size_t j = 0; j < xfunctions.size(); ++j) {
                                gaxpy(world, 1.0, new_xfunctions[i], U(j, i), xfunctions[j]);
                        }
                }

                // Return the transformed vector of vecfunctions
                return new_xfunctions;

        }
};

struct diffuse_functions : public FunctionFunctorInterface<double,3> {
public:
        diffuse_functions(const double L,const std::vector<coord_3d> coord,const std::vector<int> signs, const size_t natoms, const size_t mode):
                L(L), coord(coord), signs(signs),natoms(natoms),mode(mode)
{if(mode>1) MADNESS_EXCEPTION("mode in diffuse_function struct is not 0,1 or 2",1);}

        double operator()(const coord_3d &r)const{
                std::vector<double> diffuse_1s_functions;
                double exponent =1.0;
                if(mode==0) exponent = 15.0;
                if(mode==1) exponent = 20.0;
                for(size_t i=0;i<natoms;i++){
                        if (signs[i]!=0){
                                // make diffuse 1s functions localized at the atoms
                                double x = r[0]-coord[i][0];
                                double y = r[1]-coord[i][0];
                                double z = r[2]-coord[i][0];
                                double rad = sqrt(x*x+y*y+z*z);
                                double diffuse_tmp = signs[i]*exp(-exponent/L*rad);
                                diffuse_1s_functions.push_back(diffuse_tmp);
                        }else diffuse_1s_functions.push_back(0.0);
                }
                double result=0;
                for(size_t i=0;i<diffuse_1s_functions.size();i++) result+=diffuse_1s_functions[i];
                return result;
        }

private:
        const double L;
        const std::vector<coord_3d> coord;
        const std::vector<int> signs;
        const size_t natoms;
        const size_t mode;
        double diffuse_2s (const coord_3d &r)const{return r_function(r)*exp(-10.0/L*r_function(r));}
        double diffuse_2px(const coord_3d &r)const{return          r[0]*exp(-10.0/L*r_function(r));}
        double diffuse_2py(const coord_3d &r)const{return          r[1]*exp(-10.0/L*r_function(r));}
        double diffuse_2pz(const coord_3d &r)const{return          r[2]*exp(-10.0/L*r_function(r));}
        double diffuse_1s (const coord_3d &r)const{return               exp(-10.0/L*r_function(r));}
        double r_function(const coord_3d &r)const{return sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]);}

};



class TDA {
public:
        TDA(World &world,const SCF &calc,const vecfuncT &mos,const std::string input):
                world(world),
                dft_(false),
                calc_(calc),
                mos_(mos),
                print_grid_(false),
                guess_("physical"),
                guess_iter_(15),
                guess_mode_("physical"),
                guess_exop_("quadrupole"),
                guess_excitations_(6),
                excitations_(4),
                bsh_eps_(1.e-5),
                iter_max_(100),
                econv_(1.e-4),
                dconv_(1.e-3),
                hard_dconv_(1.e-3),
                hard_econv_(5.e-5),
                nfreeze_(0),
                plot_(false),
                debug_(false),
                only_fock_(false),
                only_GS_(false),
                on_the_fly_(true),
                read_(false),
                only_sequential_(false),
                xclib_interface_(world,calc),
                ipot_(0.0),
                rydberg_(false),
                kain_(false),
                kain_subspace_(3),
                shift_(0.0),
                safety_(1.0)
{
                setup(mos,input);
}
        void setup(const vecfuncT &mos,const std::string input){


                // so that the thresh can be changed from the outside
                mos_ = mos;

                size_t noct = calc_.aeps.size();
                // The highest possible excitation (-homo_energy)
                double highest_excitation_default = -calc_.aeps(noct-1);
                highest_excitation_ = highest_excitation_default;
                ipot_ = -calc_.aeps(noct-1)*2.0;


                // The guessed lowest excitation (if no guess_omega_ is in the input)
                double guess_omega_default = -0.1*calc_.aeps[noct-1];
                guess_omega_ = guess_omega_default;

                std::ifstream f(input.c_str());
                position_stream(f, "TDA");
                std::string s, tag;
                while (std::getline(f,s)) {
                        std::istringstream ss(s);
                        ss >> tag;
                        if (tag == "end") break;
                        else if (tag == "dft") dft_=true;
                        else if (tag == "excitations") ss >> excitations_;
                        else if (tag == "guess") ss >> guess_;
                        else if (tag == "hard_dconv") ss >> hard_dconv_;
                        else if (tag == "hard_econv") ss >> hard_econv_;
                        else if (tag == "guess_iter") ss >> guess_iter_;
                        else if (tag == "guess_omega") ss >> guess_omega_;
                        else if (tag == "guess_mode") ss >> guess_mode_;
                        else if (tag == "guess_exop") ss >> guess_exop_;
                        else if (tag == "guess_excitations") ss >> guess_excitations_;
                        else if (tag == "bsh_eps") ss >> bsh_eps_;
                        else if (tag == "iter_max") ss >> iter_max_;
                        else if (tag == "econv") ss >> econv_;
                        else if (tag == "dconv") ss >> dconv_;
                        else if (tag == "freeze") ss >> nfreeze_;
                        else if (tag == "print_grid") print_grid_=true;
                        else if (tag == "plot") plot_=true;
                        else if (tag == "debug") debug_=true;
                        else if (tag == "only_fock") only_fock_=true;
                        else if (tag == "only_GS") only_GS_=true;
                        else if (tag == "highest_excitation") ss >> highest_excitation_;
                        else if (tag == "no_otf") on_the_fly_=false;
                        else if (tag == "read") read_ = true;
                        else if (tag == "only_sequential") only_sequential_=true;
                        else if (tag == "ipot") ss >> ipot_;
                        else if (tag == "rydberg") rydberg_=true;
                        else if (tag == "kain") kain_=true;
                        else if (tag == "kain_subspace") ss>> kain_subspace_;
                        else if (tag == "exop1") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop2") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop3") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop4") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop5") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop6") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop7") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop8") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop9") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "exop10") {std::string tmp; ss >> tmp; custom_exops_.push_back(tmp);}
                        else if (tag == "truncate_safety") ss>>safety_;
                        else continue;
                }

                // make potential shift = -ipot - homo
                if(dft_) shift_= -ipot_ - get_calc().aeps[noct-1];
                highest_excitation_=highest_excitation_-shift_;


                if (world.rank() == 0) {
                        std::cout<< std::setw(60) <<"\n\n\n\n ======= TDA info =======\n\n\n" << std::endl;
                        if (nfreeze_==0) std::cout<< std::setw(40) <<"# frozen orbitals : "<<"none" << std::endl;
                        if (nfreeze_>0) std::cout<< std::setw(40) <<"# frozen orbitals : " <<  "0 to " << nfreeze_-1 << std::endl;
                        std::cout<< std::setw(40) <<"active orbitals : " << nfreeze_ << " to " << calc_.param.nalpha-1 << std::endl;
                        std::cout<< std::setw(40) << "guess from : " << guess_ << std::endl;
                        std::cout<< std::setw(40) << "Gram-Schmidt is used : " << !only_fock_ << std::endl;
                        std::cout<< std::setw(40) << "threshold 3D : " << FunctionDefaults<3>::get_thresh() << std::endl;
                        std::cout<< std::setw(40) << "energy convergence : " << econv_ << std::endl;
                        std::cout<< std::setw(40) << "max residual (dconv) : " << dconv_ << std::endl;
                        std::cout<< std::setw(40) << "number of excitations : " << excitations_ << std::endl;
                        std::cout<< std::setw(40) << "number of guess excitations : " << guess_excitations_ << std::endl;
                        std::cout<< std::setw(40) << "guessed lowest extitation energy : " << guess_omega_ << std::endl;
                        std::cout<< std::setw(40) << "guessed excitation operators : " << guess_exop_ << std::endl;
                        std::cout<< std::setw(40) << "highest possible excitation : " << highest_excitation_default << std::endl;
                        std::cout<< std::setw(40) << "used highest possible excitation : " << highest_excitation_ << std::endl;
                        std::cout<< std::setw(40) << "guessed ionization potential is : " << ipot_ << std::endl;
                        std::cout<< std::setw(40) << "potential shift is : " << shift_ << std::endl;
                        std::cout<< std::setw(40) << "guessed lowest excitation : " << guess_omega_default << std::endl;
                        std::cout<< std::setw(40) << "chosen lowest excitation : " << guess_omega_ << std::endl;
                        std::cout<< std::setw(40) << "orthonormalization : ";
                        if(only_fock_) std::cout << "only perturbed fock matrix"<< std::endl;
                        else if(only_GS_) std::cout << "only Gram-Schmidt"<< std::endl;
                        else std::cout << "use both"<< std::endl;
                        std::cout<< std::setw(40) << "potential calculation : " << "on_the_fly is " << on_the_fly_ << std::endl;
                        std::cout<< std::setw(40) << "use KAIN : " << kain_ << std::endl;
                }



                lo=get_calc().param.lo;
                bsh_eps_ = FunctionDefaults<3>::get_thresh()*0.1;
                for(size_t i=nfreeze_;i<mos_.size();i++){active_mo_.push_back(mos_[i]);}

                // Initialize the projector on the occupied space
                Projector<double,3> projector(mos_);
                rho0 = projector;

                // make the unperturbed density (closed shell)
                real_function_3d active_density = real_factory_3d(world);
                for(size_t i=0;i<active_mo_.size();i++){active_density += 2.0*active_mo_[i]*active_mo_[i];}
                active_density_ = active_density;
                real_function_3d density = real_factory_3d(world);
                for(size_t i=0;i<mos_.size();i++){density+=2.0*mos_[i]*mos_[i];}
                density_=density;
                std::cout <<std::setw(40) <<"Norm of unperturbed density is : " << density_.norm2() << std::endl;

                // Initialize the exchange intermediate
                if(not dft_) {
                        exchange_intermediate_ = make_exchange_intermediate();
                        std::cout << std::setw(40) << "CIS is used" << " : LIBXC Interface is not initialized" << std::endl;
                }if(dft_){
                        lda_intermediate_ = make_lda_intermediate();
                }

                // Initialize the KAIN solvers
                kain_solvers.initialize(world,excitations_,noct,kain_);


                // Prevent misstakes:
                if(shift_>0){MADNESS_EXCEPTION("Potential shift is positive",1);}
                if(not dft_ and shift_ !=0.0){MADNESS_EXCEPTION("Non zero potential shift in TDHF calculation",1);}

                if(only_fock_ and only_GS_){
                        print("\nWARNING: only_fock and only_GS demanded ...use both");
                        only_fock_ = false;
                        only_GS_ = false;
                }

                // make the truncate thresh
                truncate_thresh_ = FunctionDefaults<3>::get_thresh() * safety_;
                std::cout << "Truncate threshold is set to " << truncate_thresh_ << std::endl;

                // Truncate the current mos
                truncate(world,mos_,truncate_thresh_);
                std::cout << "truncate molecular orbitals to " << truncate_thresh_ << std::endl;

                std::cout << "setup of TDA class ended\n" << std::endl;
        }


        double memwatch(const xfunctionsT &xfunctions,const bool printout)const{
                // sanity_check
                if(xfunctions.empty())return 0.0;
                double allx=0.0; double allVx=0.0; double allr=0.0;
                if(printout)std::cout << "\n\n#" << "  " << "     x " << "     Vx " << "     r " << std::endl;
                if(printout)print("-------------------------------");
                for(size_t i=0;i<xfunctions.size();i++){
                        if(on_the_fly_ and not xfunctions[i].Vx.empty()) MADNESS_EXCEPTION("on the fly calculation used but Vx not empty",1);
                        if(not kain_ and not xfunctions[i].current_residuals.empty()) MADNESS_EXCEPTION("no kain is used but current residuals are not empty",1);

                        // mem information
                        double x_size = get_size(world,xfunctions[i].x);
                        double Vx_size= get_size(world,xfunctions[i].Vx);
                        double r_size=get_size(world,xfunctions[i].current_residuals);
                        allx+=x_size; allVx+=Vx_size; allr=r_size;
                        if(printout)std::cout << i << "  " << x_size <<" "<< Vx_size <<" "<< r_size << " (GB)" <<  std::endl;
                }
                if(printout)print("-------------------------------");
                if(printout)std::cout << "all" << "  " << allx <<" "<< allVx <<" "<< allr << " (GB)\n\n" <<  std::endl;
                return allx+allVx+allr;
        }

        //virtual ~TDA();

        void solve(xfunctionsT &xfunctions);

        void solve_sequential(xfunctionsT xfunctions);

        const SCF &get_calc() const {return calc_;}

        // Print out grid (e.g for Koala or other external program)
        const bool print_grid_TDA() const {return print_grid_;}

        // returns a shallow copy the converged xfunctions
        xfunctionsT get_converged_xfunctions(){return converged_xfunctions_;}

private:

        World & world;

        bool dft_;

        const SCF &calc_;

        vecfuncT mos_;

        bool print_grid_;


        std::string guess_;
        size_t guess_iter_;
        double guess_omega_;

        std::string guess_mode_;

        std::string guess_exop_;
        size_t guess_excitations_;
        std::vector<std::string> custom_exops_;

        size_t excitations_;

        double bsh_eps_;

        size_t iter_max_;
        double econv_;
        double dconv_;
        // Convergence criteria (residual norm) for the high thresh sequential iterations in the end
        double hard_dconv_;
        double hard_econv_;
        size_t nfreeze_;

        bool plot_;

        bool debug_;

        bool only_fock_;
        bool only_GS_;

        double highest_excitation_;

        double lo;

        vecfuncT active_mo_;

        // the projector on the unperturbed density
        Projector<double,3> rho0;

        real_function_3d active_density_;

        real_function_3d density_;

        bool on_the_fly_;

        bool read_;

        bool only_sequential_;

        bool sequential_;

        TDA_DFT xclib_interface_;

        double ipot_;

        bool rydberg_;

        bool kain_;
        kain_solver_helper_struct kain_solvers;

        size_t kain_subspace_;

        double shift_;

        double truncate_thresh_;

        double safety_;

        real_function_3d unperturbed_vxc_;

        vecfuncT lda_intermediate_;

        std::vector<vecfuncT> exchange_intermediate_;

        mutable real_function_3d coulomb_;

        std::vector<xfunction> converged_xfunctions_;

        void print_status(const xfunctionsT & xfunctions)const;

        void print_xfunction(const xfunction &x)const;

        void initialize(xfunctionsT & xfunctions);

        void guess_physical(xfunctionsT & xfunctions);

        void add_diffuse_functions(vecfuncT &mos);

        vecfuncT make_excitation_operators()const;

        // Calls iterate_all with the right settings
        void iterate_guess(xfunctionsT &xfunctions);

        // Calls iterate_all with the right settings
        void iterate(xfunctionsT &xfunctions);

        // guess: guess inititialisation or final iterations
        void iterate_all(xfunctionsT &xfunctions,bool guess);

        void iterate_one(xfunction & xfunction,bool ptfock,bool guess);

        // instead of x = G(-2VPsi) and G parametrized with epsilon+omega
        // here: x= G(-2VPsi + 2omega x) and G aprametrized only with epsilon (always bound)
        void iterate_one_unbound(xfunction & xfunction,bool ptfock,bool guess);

        void iterate_one_yanai(xfunction & xfunction,bool guess);

        void update_energies(xfunctionsT &xfunctions);

        void normalize(xfunctionsT &xfunctions);
        void normalize(xfunction &xfunction);

        void project_out_converged_xfunctions(xfunctionsT & xfunctions);

        void orthonormalize_GS(xfunctionsT &xfunctions);

        // 1. Call make_perturbed_fock_matrix(xfunctions)
        // 2. Diagonalize
        // 3. Update Energy and xfunctions
        bool orthonormalize_fock(xfunctionsT &xfunctions,const bool guess);

        double measure_offdiagonality(const madness::Tensor<double> &U,const size_t size)const;

        std::vector<vecfuncT> transform_vecfunctions(const std::vector<vecfuncT> &xfunctions,const madness::Tensor<double> U)const;

        // 1. Make projector
        // 2. apply
        void project_out_occupied_space(vecfuncT &x);

        double perturbed_fock_matrix_element(const vecfuncT &xr, const vecfuncT &Vxp,const vecfuncT &xp)const;

        // Can also be used to calculate diagonal elements of the fock matrix
        double expectation_value(const xfunction &x,const vecfuncT &Vx);

        // 1. Calculate potential (V0 + Gamma) --> Call perturbed_potential(exfunctions)
        // 2. Calculate Matrix Elements
        Tensor<double> make_perturbed_fock_matrix(const xfunctionsT &xfunctions)const;

        // 1. Call get_V0
        // 2. Call apply_gamma or apply_gamma_dft
        // 3. return V0 + Gamma
        vecfuncT apply_perturbed_potential(const xfunction & xfunction)const;

        vecfuncT apply_gamma(const xfunction &xfunction)const;
        vecfuncT apply_gamma_dft(const xfunction &xfunction)const;

        vecfuncT apply_hartree_potential(const vecfuncT &x)const;

        // This has to be done just one time because only the unperturbed orbitals are needed
        std::vector<vecfuncT> make_exchange_intermediate()const;

        vecfuncT make_lda_intermediate()const;

        vecfuncT get_V0(const vecfuncT &x)const;

        real_function_3d get_coulomb_potential() const;

        real_function_3d get_vxc_potential()const;

        void plot_vecfunction(const vecfuncT &x,std::string msg, bool plot = true)const;

        bool check_convergence(xfunctionsT &xfunctions);

        void print_performance(const xfunctionsT &xfunctions,const std::string prename)const;

        void truncate_xfunctions(xfunctionsT &xfunctions);



        double oscillator_strength_length(const xfunction& xfunction) const;


        double oscillator_strength_velocity(const xfunction& root) const;

        void analyze(xfunctionsT& roots) const;

        void save_xfunctions(const xfunctionsT &xfunctions)const;
        bool read_xfunctions(xfunctionsT &xfunctions);
};

} /* namespace madness */

#endif /* TDA_H_ */