mp2.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$
*/
#ifndef MP2_H_
#define MP2_H_

//#define WORLD_INSTANTIATE_STATIC_TEMPLATES
#include <madness/mra/mra.h>
#include <madness/mra/lbdeux.h>
#include <chem/SCF.h>
#include <examples/nonlinsol.h>
#include <chem/projector.h>
#include <chem/correlationfactor.h>
#include <chem/nemo.h>

#include <iostream>


using namespace madness;

namespace madness {


    struct LBCost {
        double leaf_value;
        double parent_value;
        LBCost(double leaf_value=1.0, double parent_value=1.0)
            : leaf_value(leaf_value)
            , parent_value(parent_value)
        {}

        double operator()(const Key<6>& key, const FunctionNode<double,6>& node) const {
            return node.coeff().rank();
//            if (node.is_leaf()) {
//                return leaf_value;
//            } else {
//                return parent_value;
//            }
        }
    };


    class HartreeFock {
        World& world;
        std::shared_ptr<SCF> calc;
        mutable double coords_sum;     // sum of square of coords at last solved geometry

        // save the Coulomb potential
        mutable real_function_3d coulomb;

        std::vector<real_function_3d> orbitals_;

        std::vector<real_function_3d> R2orbitals_;

    public:

        Nemo nemo_calc;

        HartreeFock(World& world, std::shared_ptr<SCF> calc1) :
                world(world), calc(calc1), coords_sum(-1.0), nemo_calc(world,calc1) {
        }

        bool provides_gradient() const {return true;}

        double value() {
            return value(calc->molecule.get_all_coords());
        }

        double value(const Tensor<double>& x) {

                // fast return if the reference is already solved at this geometry
            double xsq = x.sumsq();
            if (xsq == coords_sum) return calc->current_energy;

            calc->molecule.set_all_coords(x.reshape(calc->molecule.natom(),3));
            coords_sum = xsq;

            // some extra steps if we have a nuclear correlation factor
            if (1) {
//              if (nemo_calc.nuclear_correlation->type()==NuclearCorrelationFactor::GaussSlater) {
                        // converge the nemo equations
                        nemo_calc.value(x);

                } else {
                                // Make the nuclear potential, initial orbitals, etc.
                                calc->make_nuclear_potential(world);
                                calc->potentialmanager->vnuclear().print_size("vnuc");
                                calc->project_ao_basis(world);

                                // read converged wave function from disk if there is one
                                if (calc->param.no_compute) {
                                        calc->load_mos(world);
                                        return calc->current_energy;
                                }

                                if (calc->param.restart) {
                                        calc->load_mos(world);
                                } else {
                                        calc->initial_guess(world);
                                        calc->param.restart = true;
                                }

                                // If the basis for the inital guess was not sto-3g
                                // switch to sto-3g since this is needed for analysis
                                // of the MOs and orbital localization
                                if (calc->param.aobasis != "sto-3g") {
                                        calc->param.aobasis = "sto-3g";
                                        calc->project_ao_basis(world);
                                }


                                calc->solve(world);
                                calc->save_mos(world);

                                // successively tighten threshold
                                if (calc->param.econv<1.1e-6) {
                                        calc->set_protocol<3>(world,1e-6);
                                        calc->make_nuclear_potential(world);
                                        calc->project_ao_basis(world);
                                        calc->project(world);
                                        calc->solve(world);
                                        calc->save_mos(world);
                                }

                                calc->save_mos(world);
                        }

            // compute the full, reconstructed orbitals from nemo
            orbitals_=mul(world,nemo_calc.R,nemo_calc.get_calc()->amo);
            real_function_3d R2=nemo_calc.nuclear_correlation->square();
            R2orbitals_=mul(world,R2,nemo_calc.get_calc()->amo);

            return calc->current_energy;
        }

        Tensor<double> gradient(const Tensor<double>& x) {

            value(x); // Ensures DFT equations are solved at this geometry
            return calc->derivatives(world);
        }

        double coord_chksum() const {return coords_sum;}

        const SCF& get_calc() const {return *calc;}
        SCF& get_calc() {return *calc;}


        real_function_3d orbital(const int i) const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return orbitals_[i];
        }


        std::vector<real_function_3d> orbitals() const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return orbitals_;
        }


        std::vector<real_function_3d> R2orbitals() const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return R2orbitals_;
        }


        real_function_3d R2orbital(const int i) const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return R2orbitals_[i];
        }


        real_function_3d nemo(const int i) const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return calc->amo[i];
        }


        std::vector<real_function_3d> nemos() const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return calc->amo;
        }

        double orbital_energy(const int i) const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return calc->aeps[i];
        }

        real_function_3d get_coulomb_potential() const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            if (coulomb.is_initialized()) return copy(coulomb);
            functionT rho = calc->make_density(world, calc->aocc, orbitals()).scale(2.0);
            coulomb=calc->make_coulomb_potential(rho);
            return copy(coulomb);
        }

        real_function_3d get_nuclear_potential() const {
            return calc->potentialmanager->vnuclear();
        }

        int nocc() const {
            MADNESS_ASSERT(calc->param.spin_restricted);
            return calc->param.nalpha;
        }
    };


    class ElectronPair : public archive::ParallelSerializableObject {

    public:
        ElectronPair()
                : i(-1), j(-1), e_singlet(uninitialized()), e_triplet(uninitialized()),
                        ij_gQf_ij(uninitialized()), ji_gQf_ij(uninitialized()), iteration(0), converged(false) {
        }

        ElectronPair(const int i, const int j)
                : i(i), j(j), e_singlet(uninitialized()), e_triplet(uninitialized()),
                  ij_gQf_ij(uninitialized()), ji_gQf_ij(uninitialized()), iteration(0), converged(false) {
        }

        void print_energy() const {
            if (function.world().rank()==0) {
                printf("final correlation energy %2d %2d %12.8f %12.8f\n",
                                i,j,e_singlet,e_triplet);
            }
        }

        static double uninitialized() {return 1.e10;}

        int i, j;                       
        real_function_6d function;      
        real_function_6d r12phi;        
        real_function_6d constant_term; 

        real_function_6d Uphi0;         
        real_function_6d KffKphi0;      
        std::vector<real_function_3d> phi_k_UK_phi0;    
        std::vector<real_function_3d> phi_l_UK_phi0;    

        double e_singlet;                               
        double e_triplet;                               

        double ij_gQf_ij;                       
        double ji_gQf_ij;                       

        int iteration;                                  
        bool converged;                                 


        template <typename Archive> void serialize (Archive& ar) {
                bool fexist=function.is_initialized();
                bool cexist=constant_term.is_initialized();
                        ar & ij_gQf_ij & ji_gQf_ij & e_singlet & e_triplet & converged
                                & iteration & fexist & cexist;
                        if (fexist) ar & function;
                        if (cexist) ar & constant_term;
        }

        bool load_pair(World& world) {
                std::string name="pair_"+stringify(i)+stringify(j);
                bool exists=archive::ParallelInputArchive::exists(world,name.c_str());
            if (exists) {
                if (world.rank()==0) printf("loading matrix elements %s",name.c_str());
                archive::ParallelInputArchive ar(world, name.c_str(), 1);
                ar & *this;
                if (world.rank()==0) printf(" %s\n",(converged)?" converged":" not converged");
            } else {
                        if (world.rank()==0) print("could not find pair ",i,j," on disk");
            }
            return exists;
        }

        void store_pair(World& world) {
                std::string name="pair_"+stringify(i)+stringify(j);
                if (world.rank()==0) printf("storing matrix elements %s\n",name.c_str());
            archive::ParallelOutputArchive ar(world, name.c_str(), 1);
                ar & *this;
        }
    };


    class MP2 : public OptimizationTargetInterface {

        struct Parameters {
                double thresh_;                 
                double dconv_;                  
                int i;                                  
                int j;                                  

                int freeze;


                bool restart;

                int maxsub;

                int maxiter;

                Parameters(const std::string& input) : thresh_(-1.0), dconv_(-1.0),
                                i(-1), j(-1), freeze(0), restart(false), maxsub(2), maxiter(20) {

                        // get the parameters from the input file
                std::ifstream f(input.c_str());
                position_stream(f, "mp2");
                std::string s;

                while (f >> s) {
                    if (s == "end") break;
                    else if (s == "econv") f >> thresh_;
                    else if (s == "dconv") f >> dconv_;
                    else if (s == "pair") f >> i >> j;
                    else if (s == "maxsub") f >> maxsub;
                    else if (s == "freeze") f >> freeze;
                    else if (s == "restart") restart=true;
                    else continue;
                }
                // set default for dconv if not explicitly given
                if (dconv_<0.0) dconv_=sqrt(thresh_)*0.1;
                }

                void check_input(const std::shared_ptr<HartreeFock> hf) const {
                if (freeze>hf->nocc()) MADNESS_EXCEPTION("you froze more orbitals than you have",1);
                if (i>=hf->nocc()) MADNESS_EXCEPTION("there is no i-th orbital",1);
                if (j>=hf->nocc()) MADNESS_EXCEPTION("there is no j-th orbital",1);
                if (thresh_<0.0) MADNESS_EXCEPTION("please provide the accuracy threshold for MP2",1);
                }
        };


        World& world;                           
        Parameters param;                                               
        std::shared_ptr<HartreeFock> hf;        
        CorrelationFactor corrfac;              
        std::shared_ptr<NuclearCorrelationFactor> nuclear_corrfac;

        std::map<std::pair<int,int>,ElectronPair> pairs;       
        double correlation_energy;                              
        double coords_sum;                                              
        bool do_coupling;

        StrongOrthogonalityProjector<double,3> Q12;

    private:
        struct Intermediates {
            std::string function;      
            std::string r12phi;        
            std::string Kfphi0;        
            std::string Uphi0;         
            std::string KffKphi0;      
            std::string OUKphi0;                

            Intermediates() : r12phi(), Kfphi0(), Uphi0(), KffKphi0() {};

            Intermediates(World& world, const std::string& filename) : function(), r12phi(),
                    Kfphi0(), Uphi0(), KffKphi0() {
                std::ifstream f(filename.c_str());
                position_stream(f, "mp2");
                std::string s;

                while (f >> s) {
                    if (s == "end") break;
                    else if (s == "function") f >> function;
                    else if (s == "r12phi") f >> r12phi;
                    else if (s == "Kfphi0") f >> Kfphi0;
                    else if (s == "Uphi0") f >> Uphi0;
                    else if (s == "KffKphi0") f >> KffKphi0;
                    else {continue;
                    }
                    if (world.rank()==0) print("found intermediate in control file: ",s);
                }
            }

            template <typename Archive> void serialize (Archive& ar) {
                ar & function & r12phi & Kfphi0 & Uphi0 & KffKphi0;
            }

        };


        Intermediates intermediates;
        std::shared_ptr<real_convolution_3d> poisson;

    public:

        MP2(World& world, const std::string& input);


        ElectronPair& pair(const int i, const int j) {
                // since we return a reference the keyval must already exist in the map
                MADNESS_ASSERT(pairs.find(std::make_pair(i, j)) != pairs.end());
                return pairs[std::make_pair(i, j)];
        }


        const ElectronPair& pair(const int i, const int j) const {
                // since we return a reference the keyval must already exist in the map
                MADNESS_ASSERT(pairs.find(std::make_pair(i, j)) != pairs.end());
                return pairs.find(std::make_pair(i, j))->second;
        }

        double coord_chksum() const {return coords_sum;}

        double value();

        double value(const Tensor<double>& x);

        void print_info(World& world) const;

        HartreeFock& get_hf() {return *hf;}

        double thresh() const {return param.thresh_;}

        double zeroth_order_energy(const int i, const int j) const {
            return hf->orbital_energy(i)+hf->orbital_energy(j);
        }

        void solve_residual_equations(ElectronPair& pair) const;

        real_function_6d make_Rpsi(const ElectronPair& pair) const;


        void increment(ElectronPair& pair, real_convolution_6d& green);


        real_function_6d swap_particles(const real_function_6d& f) const;

        double asymmetry(const real_function_6d& f, const std::string s) const;

        void test(const std::string filename);


        double compute_gQf(const int i, const int j, ElectronPair& pair) const;

    private:

        template<typename T, size_t NDIM>
        void save_function(const Function<T,NDIM>& f, const std::string name) const;

        template<typename T, size_t NDIM>
        void load_function(Function<T,NDIM>& f, const std::string name) const;

        real_function_6d make_Uphi0(ElectronPair& pair) const;

        real_function_6d make_KffKphi0(const ElectronPair& pair) const;

        ElectronPair make_pair(const int i, const int j) const;

        void guess_mp1_3(ElectronPair& pair) const;


        double compute_energy(ElectronPair& pair) const;


        real_function_3d K(const real_function_3d& phi, const bool hc=false) const;


        real_function_6d K(const real_function_6d& phi, const bool is_symmetric=false) const;


        real_function_3d J(const real_function_3d& phi) const;


        real_function_6d apply_exchange(const real_function_6d& f,
                        const real_function_3d& orbital_ket,
                        const real_function_3d& orbital_bra, const int particle) const;


        std::vector<real_function_3d> make_chi(const real_function_3d& phi,
                        const real_convolution_3d& op,
                        const bool hc=false) const;


        std::vector<real_function_3d> make_xi(const real_function_3d& phi_i,
                        const real_function_3d& phi_j, const real_convolution_3d& op,
                        const bool hc=false) const;

        void OUKphi0(ElectronPair& pair) const;


        real_function_6d make_fKphi0(const int i, const int j) const;


        real_function_6d JK1phi0_on_demand(const int i, const int j,
                        const bool hc=false) const;

        real_function_6d JK2phi0_on_demand(const int i, const int j,
                        const bool hc=false) const;

        real_function_6d phi0_on_demand(const int i, const int j) const;

        real_function_6d nemo0_on_demand(const int i, const int j) const;



        real_function_6d multiply_with_0th_order_Hamiltonian(const real_function_6d& f,
                        const int i, const int j) const;


        real_function_6d add_local_coupling(const int i, const int j) const;

    };
};

#endif /* MP2_H_ */