tdhf_CIS.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.

  Copyright (C) 2007,2010 Oak Ridge National Laboratory

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

  For more information please contact:

  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367

  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680

  $Id$
*/
/* * tdhfCIS.h
 *
 *  Created on: May 5, 2014
 *      Author: kottmanj
 */

#ifndef TDHFCIS_H_
#define TDHFCIS_H_

#include <chem/projector.h>
//#include <examples/mp2.h>

#include<examples/nonlinsol.h>
#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>




using namespace madness;


typedef std::vector<Function<double,3> > vecfuncT;

// This is not used anymore, but maybe will be useful later
struct gauss_function : public FunctionFunctorInterface<double,3> {

        typedef std::shared_ptr<FunctionFunctorInterface<double,3> > functorT;

public:
        gauss_function(double molecule_box) : sigma(molecule_box){}

private:
        real_function_3d function;
        double sigma;


public:
        double operator()(const coord_3d &r)const {
                double r2=r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
                double prefactor = 1.0/(sqrt(2.0*constants::pi)*sigma);
                return prefactor*exp(-r2/(2.0*sigma*sigma));
        }
};

// Not used anymore but maybe useful later
struct noise : public FunctionFunctorInterface<double,3> {

typedef std::shared_ptr<FunctionFunctorInterface<double,3> > functorT;

public:
        noise(double size){noise(size,1.0);}
        noise(double size,double width)  {
         Random random(wall_time()*1.e8);
         x=size*(random.get()-random.get());
         y=size*(random.get()-random.get());
         z=size*(random.get()-random.get());
         prefactor= x*y*z/(fabs(x*y*z));
         sigma=width;
        // print("Noise constructor x,y,z,pre,sigma are: ",x,y,z,prefactor,sigma);
        }


private:
        real_function_3d gauss_function;
        double x,y,z;
        double sigma;
        double prefactor;

public:
        double operator()(const coord_3d &r) const{
                double r2 = (r[0]-x)*(r[0]-x)+(r[1]-y)*(r[1]-y)+(r[2]-z)*(r[2]-z);
                return prefactor*exp(-r2/(2.0*sigma*sigma));
        }

};


struct root {
        root(World& world) : world(world), omega(0.0),expv(0.0),converged(false),err(10.0),delta(10.0),iter(0) {}
        root(World& world,vecfuncT& x, double omega) :world(world), x(x), omega(omega),converged(false),err(10.0),delta(10.0),iter(0) {}
        root(World& world, const vecfuncT& x1) : world(world), x(x1),converged(false),err(10.0),delta(10.0),iter(1){}
    root(const root& other) : world(other.world), x(other.x),
                omega(other.omega),expv(other.expv), converged(other.converged),
                err(other.err),delta(other.delta),iter(other.iter),number(other.number){}

    // Constructor for the print routine to create a copy without the x functions
    root(World &world,double omega,double expv, double delta, double error, bool converged,int iter, int number) : world(world),
                omega(omega),expv(expv),converged(converged),err(error),delta(delta),iter(iter),number(number) {}

        World &world;
        vecfuncT x;
        double omega;
        double expv; // expectation_value
        bool converged;
        double err;
        double delta;
        int iter;
        int number;
        //solverT solver; // Maybe later
        std::vector<double> amplitudes_;

        // Operators needed by the nonlinear solver
    root& operator=(const root& other) {
        x=other.x;
        // This is new, needed for sort function
        // Not shure if this changes something when KAIN is used again
        omega=other.omega;
        expv=other.expv;
        converged=other.converged;
        err=other.err;
        delta=other.delta;
        amplitudes_=other.amplitudes_;
        number=other.number;
        iter=other.iter;
        return *this;
    }

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

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

    root operator*(double a) { // Scale by a constant necessary
        scale(world,x,a);
        return *this;
    }
    // An operator for the sort function
    bool operator<(const root &other)const{
        return (this->omega<other.omega);
    }
    bool operator>(const root &other)const{
        return (this->omega>other.omega);
    }
};


namespace madness {




class CIS {

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

private:
        #define TRUE  1
        #define FALSE 0

        static double rfunction(const coord_3d& r) {
            return sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]);
        }

        static double rfunction2(const coord_3d& r) {
            return sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2])+1.0;
        }

        // r^2
        static double monopole(const coord_3d &r){
                return r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
        }

        static double x(const coord_3d &r){return r[0];}
        static double x2(const coord_3d &r){return r[0]-0.7;}
        static double x3(const coord_3d &r){return r[0]+0.7;}
        static double y(const coord_3d &r){return r[1];}
        static double z(const coord_3d &r){return r[2];}

        static double b1u(const coord_3d &r){return r[0]*(r[0]*r[0]-3*r[1]*r[1]);}
        static double b2u(const coord_3d &r){return r[1]*(3*r[0]*r[0]-r[1]*r[1]);}
        static double e1g1(const coord_3d &r){return r[0]*r[2];}
        static double e1g2(const coord_3d &r){return r[1]*r[2];}
        static double a2u(const coord_3d &r){return r[2];}
        static double e2u1(const coord_3d &r){return r[2]*(r[0]*r[0]-r[1]*r[1]);}
        static double e2u2(const coord_3d &r){return r[0]*r[1]*r[2];}

        // Random number Generator from Madness lib misc/ran.h
        static double random_number(const coord_3d &r){
                Random random;
                return random.get();
        }



        // For the sorting of the roots
        static bool compare_roots(const root &a,const root &b){
                return a.omega<b.omega;
        }
        static bool compare_roots_error(const root &a,const root &b){
                return a.err<b.err;
        }


public:


        CIS(World& world, const SCF& calc, const std::string input)
: world(world),
  calc_(calc),
  guess_("physical"),
  nroot_(8),
  nfreeze_(0),
  econv_(calc.param.econv),
  dconv_(calc.param.econv),
  guess_thresh_(1.e-4),
  guess_econv_(0.001),
  guess_dconv_(0.03),
  guess_iter_(20),
  guess_iter_fock_(5),
  guess_roots_(8),
  guess_mode_("all_orbitals"),
  thresh_(dconv_*0.01),
  bsh_eps_(1.e-6),
  iter_max_(50),
  exchange_("hf"),
  print_grid_(false),
  plot_(false),
  fixed_point_(false),

  guess_save_(false),
  guess_damp_(false),
  guess_pull_(false),
  guess_damp_iter_(6),
  guess_damp_fock_(true),

  noise_(false),
  noise_box_(get_calc().molecule.bounding_cube()),
  noise_comp_(0.01),
  noise_width_(3.5),
  noise_gaussnumber_(25),
  hf_(false),
  triplet_(false),
  read_and_save_koala_(false),
  analyze_("false"),
  guess_exf_("dipole"){

                size_t nmo = get_calc().amo.size();
                guess_omega_=-0.9*get_calc().aeps(nmo-1);
                omega_=std::vector<double>(9,100.0);
                active_mo_ = nmo-nfreeze_;

                std::ifstream f(input.c_str());
                position_stream(f, "CIS");
                std::string s, tag;
                while (std::getline(f,s)) {
                        std::istringstream ss(s);
                        ss >> tag;
                        if (tag == "end") break;
                        else if (tag == "guess") ss >> guess_;
                        else if (tag == "guess_preopt") ss >> preopt_;
                        else if (tag == "nroot") ss >> nroot_;
                        // dangerous right now else if (tag == "guess_roots") ss >> guess_roots_;
                        else if (tag == "guess_thresh") ss >> guess_thresh_;
                        else if (tag == "guess_econv") ss >> guess_econv_;
                        else if (tag == "guess_dconv") ss >> guess_dconv_;
                        else if (tag == "guess_iter") ss >> guess_iter_;
                        else if (tag == "guess_iter_fock") ss >> guess_iter_fock_;
                        else if (tag == "guess_roots") ss >> guess_roots_;
                        else if (tag == "guess_omega") ss >> guess_omega_;
                        else if (tag == "guess_mode") ss >> guess_mode_;
                        else if (tag == "thresh") ss >> thresh_;
                        else if (tag == "bsh_eps") ss >> bsh_eps_;
                        else if (tag == "iter_max") ss >> iter_max_;
                        else if (tag == "freeze") ss >> nfreeze_;
                        else if (tag == "econv") ss >> econv_;
                        else if (tag == "dconv") ss >> dconv_;
                        else if (tag == "fixed_point") fixed_point_=true;
                        else if (tag == "print_grid") print_grid_=true;
                        else if (tag == "plot") plot_=true;
                        else if (tag == "guess_save") guess_save_=true;
                        else if (tag == "guess_damp") guess_damp_=true;
                        else if (tag == "guess_pull") guess_pull_=true;
                        else if (tag == "guess_damp_iter") ss >> guess_damp_iter_;
                        else if (tag == "guess_damp_fock") guess_damp_fock_ = true;
                        else if (tag == "noise") noise_=true;
                        else if (tag == "noise_box") ss >> noise_box_;
                        else if (tag == "noise_comp") ss >> noise_comp_;
                        else if (tag == "noise_width") ss >> noise_width_;
                        else if (tag == "noise_gaussnumber") ss >> noise_gaussnumber_;
                        else if (tag == "hf") hf_=true;
                        else if (tag == "triplet") triplet_=true;
                        else if (tag == "koala_read_and_save") read_and_save_koala_=true;
                        else if (tag == "guess_exf") ss >> guess_exf_;
                        else if (tag == "omega0") ss >> omega_[0];
                        else if (tag == "omega1") ss >> omega_[1];
                        else if (tag == "omega2") ss >> omega_[2];
                        else if (tag == "omega3") ss >> omega_[3];
                        else if (tag == "omega4") ss >> omega_[4];
                        else if (tag == "omega5") ss >> omega_[5];
                        else if (tag == "omega6") ss >> omega_[6];
                        else if (tag == "omega7") ss >> omega_[7];
                        else if (tag == "omega8") ss >> omega_[8];
                        else if (tag == "active_mo") ss >> active_mo_;
                        else if (tag == "active_mo") guess_mode_="active_space";
                        else if (tag == "analyze") analyze_=true;
                        else continue;
                }






                if (world.rank() == 0) {
                        madness::print("\n ======= CIS info =======\n");
                        if (nfreeze_==0) madness::print("   # frozen orbitals ","none");
                        if (nfreeze_>0) madness::print("   # frozen orbitals ",0, " to ",nfreeze_-1);
                        madness::print("     active orbitals ", nfreeze_," to ",get_calc().param.nalpha-1);

                        madness::print("          guess from ", guess_);
                        madness::print("        threshold 3D ", FunctionDefaults<3>::get_thresh());
                        madness::print("  energy convergence ", econv_);
                        madness::print("max residual (dconv) ", dconv_);
                        madness::print("     number of roots ", nroot_);
                        madness::print(" omega ", omega_[0],omega_[1],omega_[2]);

                }

                lo=get_calc().param.lo;
        }


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

        // print information of root or root vector
        void print_roots(const std::vector<root> &roots) const;
        void print_roots(const std::vector<root> &roots,const int iter) const;
        void print_root(const root &root) const;

        void sort_roots(std::vector<root> & roots,std::string criterium)const;

        // read and analyze roots
        void Analyze();

        void solve();

        // Internal solver for parallel or sequential optimization
        bool solve_internal_par(const std::string mode,std::vector<root> &roots,const int iter_max);

        std::vector<root>& roots();

        bool print_grid() const;

private:

        World& world;

        const SCF& calc_;

        std::vector<root> roots_;


        std::vector<vecfuncT> exchange_intermediate_;

        mutable real_function_3d coulomb_;


        std::string guess_;

        std::string preopt_;

        int nroot_;

        int nfreeze_;

        double econv_;

        double dconv_;

        Tensor<double> guess_phases_;

        // Thresh for guess calculation
        double guess_thresh_;
        double guess_econv_;

        double guess_dconv_;

        int guess_iter_;
        int guess_iter_fock_;

        int guess_roots_;

        double guess_omega_;

        std::string guess_mode_;

        double thresh_;

        double bsh_eps_;

        int iter_max_;

        std::vector<double> omega_;

        std::string exchange_;



        bool print_grid_;

        bool plot_;

        bool fixed_point_;

        bool guess_save_;

        // Not used anymore, but maybe again with DFT
        bool guess_damp_;
        bool guess_pull_;
        int guess_damp_iter_;
        bool guess_damp_fock_;

        bool noise_;
        double noise_box_;
        double noise_comp_;
        double noise_width_;
        double noise_gaussnumber_;

        bool hf_;

        // Calculating triplets does not work right now
        bool triplet_;

        bool read_and_save_koala_;

        int active_mo_;

        // just read and alalyze the roots
        bool analyze_;

        std::string guess_exf_;

        double lo;

        void initialize_roots(World &world,std::vector<root> &roots);

        void guess_MO(World &world,std::vector<root> &roots);
        void guess_read(World &world,std::vector<root> &roots);
        // you need to print out the grid before and then use koala
        // 1. Madness tdhf calculation with keyword print_grid in CIS
        // 2. Koala calculation using this grid
        // 3. Madness tdhf calculation with keyword guess koala in CIS
        void guess_koala(World &world,std::vector<root> &roots);
        void guess_physical(World &world,std::vector<root> &roots);

        void guess_noise(World &world,std::vector<root> & roots);
        void guess_aspace(World &world,std::vector<root> & roots);
        void guess_forced(World &world,std::vector<root> & roots);
        void guess_benzene_custom(World &world,std::vector<root> & roots);

        // Add noise to a vector of roots (noise are random gauss functions)
        void add_noise(World & world,std::vector<root> &roots)const;

        std::vector<root> guess();

        std::vector<root> guess_big(const std::string exf);

        // Excitation functions for the guess (x,y,z,r,x^2,xy,xz ...)
        std::vector<real_function_3d> excitation_functions(const std::string exf) const;

        bool orthogonalize_fock(World& world,std::vector<root> &roots) const;
        bool orthogonalize_fock(World& world,std::vector<root> &roots,int iter) const;


        template<typename solverT>
        bool iterate_all_CIS_roots(World& world, std::vector<solverT>& solver,
                        std::vector<root>& roots,const std::string mode,const int iter_max) const;

        bool check_convergence(std::vector<root> &roots)const;


        template<typename solverT>
        double iterate_one_CIS_root(World& world, solverT& solver, root& thisroot,const std::string mode)const;

        template<typename solverT>
        double iterate_one_CIS_root(World& world, solverT& solver, root& thisroot,const std::string mode,vecfuncT &Vphi)const;

        // Update the Gamma Potential
        vecfuncT gamma_update(World &world, const root root)const;



        vecfuncT apply_gamma(const vecfuncT& x, const vecfuncT& act,
                        const Projector<double,3>& rho0) const ;

        vecfuncT apply_fock_potential(const vecfuncT& x) const;

        real_function_3d get_coulomb_potential() const;

        std::vector<vecfuncT> make_exchange_intermediate(const vecfuncT& active_mo,
                        const vecfuncT& amo) const;

        double expectation_value(World &world,const root &thisroot,const vecfuncT &Vx)const;


        Tensor<double> make_perturbed_fock_matrix(const std::vector<root>& roots, const std::vector<vecfuncT> &V) const;
                        //const vecfuncT& act, const Projector<double,3>& rho0) const;


        bool load_root(World& world, const int i, root& root)  const;
        bool load_root(World& world, const int i, root& root,const std::string filename_end)  const;


        void save_root(World& world, const int i, const root& root) const;
        void save_root(World& world, const int i, const root& root, const std::string filename_end) const;


        void normalize(World& world, root& x) const;

        // For the calculation of the guess energy
        void orthonormalize_fock(World &world,std::vector<root> &roots)const;

        void orthonormalize_fock(World &world,std::vector<root> &roots, const std::vector<vecfuncT> &Vphi)const;

        void orthonormalize(World& world, std::vector<root>& roots) const;

        Tensor<double> overlap(const std::vector<root>& r1,
                        const std::vector<root>& r2) const;



        double oscillator_strength_length(const root& root) const;


        double oscillator_strength_velocity(const root& root) const;

        void analyze(const std::vector<root>& roots) const;

        const vecfuncT active_mo() const;

};

}




#endif /* TDHFCIS_H_ */

---

Generated on Fri Sep 26 2014 10:38:58 for MADNESS by   1.8.8