lmp2.h

 
/*
 * Copyright 2009 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This file is a part of the MPQC LMP2 library.
 *
 * The MPQC LMP2 library is free software: you can redistribute it
 * and/or modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation, either
 * version 3 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
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this program.  If not, see
 * <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef _chemistry_qc_lmp2_lmp2_h
#define _chemistry_qc_lmp2_lmp2_h
 
#include <math.h>
 
#include <math/optimize/diis.h>
#include <math/scmat/repl.h>
#include <math/scmat/matrix.h>
 
#include <chemistry/qc/basis/basis.h>
#include <chemistry/qc/basis/tbint.h>
#include <chemistry/qc/scf/scf.h>
 
#include <chemistry/qc/lmp2/sma.h>
#include <chemistry/qc/lmp2/domain.h>
#include <chemistry/qc/lmp2/util.h>
#include <chemistry/qc/lmp2/lcorr.h>
#include <chemistry/qc/lmp2/parallel.h>
 
namespace sc {
 
typedef std::vector<std::pair<int,int> > my_occ_pairs_t;
 
typedef std::set<std::pair<int,int> > k_2occ_local_pairs_t;
 
typedef std::set<std::pair<int,int> > k_3_4_occ_local_pairs_t;
 
typedef std::map<sma2::triplet<int,int,int>, domainmapvirbs_t > domainmap_triple;
 
class LMP2: public LCorr {
 
    sc::Ref<sc::MessageGrp> msg_;
 
    sc::Ref<sc::SCF> ref_;
    
    int max_iter_;
 
    int nfzc_;
    
    double bound_;
    double S_threshold_;
    double integral_threshold_;
    double q1_threshold_;
    double q2_threshold_;
    double q3_threshold_;
    double q4_threshold_;
    double threshold_factor_;
    double distance_threshold_;
    double completeness_threshold_;
 
    // "canonical" or "pipek-mezey"
    std::string occ_orbitals_;
    // "canonical" or "projected_atomic"
    std::string vir_orbitals_;
    
    int parallel_transform123a_;
    int parallel_transform123b_;
    int parallel_transform4a_;
    int parallel_transform4b_;
    int parallel_iterations_;
    int completedomains_;
    bool singlevirb_;
 
    bool i_ge_j_;
    bool m_ge_n_;
    bool minimize_q2_;
    bool always_use_dist_t_;
 
    sc::Ref<sc::SelfConsistentExtrapolation> extrap_T_;
 
    sma2::Range ao_;
    sma2::Range vir_;
    sma2::Range occ_act_;
    sma2::Array<2> P_;
    sma2::Array<2> L_;
    domainmap_t domainmap_;
    my_occ_pairs_t my_occ_pairs_;
    std::vector<std::set<int> > paired_occ_;
 
    // Maps occ pairs to node numbers.
    sc::Ref<sma2::PairMapping> pair_mapping_;
 
    sma2::Array<4> K_2occ_;
 
    sma2::Array<4> T_;
    sma2::Array<4> T_jirs_;
    sma2::Array<4> T_local_;
    double T_n_element_;
 
    sma2::Array<2> F_occ_;
    sma2::Array<2> F_vir_;
    sma2::Array<2> S_;
 
    std::vector<double> F_diag_;
 
    void clear();
 
    bool analyze_occ_orbs_;
    std::vector<SCVector3> r_i_;      
    std::vector<double> dist_ij_;     
    std::vector<double> emp2_ij_;     
    void analyze_occ_orbs(const RefSCMatrix& scf_local);
 
    double compute_lmp2_energy();
    double compute_ecorr_lmp2();
 
    void rearrange_q2_all_ij(sma2::Array<4> &q2_K_oo);
    void rearrange_q2_i_ge_j(sma2::Array<4> &q2_K_oo);
 
    void compute_K_2occ(
        sc::RefSCMatrix &T_schwarz, sc::RefSCVector &T_schwarz_maxvec,
        double T_schwarz_maxval, sc::RefSCMatrix &Pmax,
        sc::RefSCVector &Pmaxvec,
        std::vector<std::vector<double> > &Lmax,
        sc::RefSCMatrix &Lshellmax,
        std::vector<std::vector<double> > &PPmax,
        const std::vector<std::vector<double> > &Dmax, double &Dmax_maxval, std::vector<double> &Lmaxvec,
        std::vector<std::multimap<double,int,std::greater<double> > > &L_map);
 
    void compute_Schwarz_screening_quantities(sc::Ref<sc::TwoBodyInt> &tbint,
                                              sc::RefSCMatrix &T_schwarz,
                                              sc::RefSCVector &T_schwarz_maxvec,
                                              double &T_schwarz_maxval,
                                              int nshell);
    
    void compute_P_screening_quantities(sc::RefSCMatrix &Pmax,
                                        sc::RefSCMatrix &Pmatrix,
                                        sc::RefSCVector &Pmaxvec,
                                        double &Pmax_maxval);
 
    void compute_L_and_D_screening_quantities(std::vector<std::vector<double> >  &Dmax,
                                              double &Dmax_maxval,
                                              std::vector<std::vector<double> > &Lmax,
                                              std::vector<double> &Lmaxvec,
                                              sc::RefSCMatrix &Lshellmax,
                                              std::vector<std::vector<int> > &L_blocks,
                                              std::vector<std::multimap<double,int,std::greater<double> > > &L_map,
                                              double &Lmax_maxval,
                                              double T_schwarz_maxval,
                                              double Pmax_maxval);
 
    void compute_PPmax_screening_matrix(std::vector<std::vector<double > > &PPmax,
                                        sc::RefSCMatrix &Pmax,
                                        const std::vector<std::set<int> >
                                           &virb_united_pair_domains,
                                        double T_schwarz_maxval,
                                        double Pmax_maxval,
                                        double Lmax_maxval);
 
    void compute_doubles_W();
 
    void compute_LMP2_residual(sma2::Array<4> &Res);
    void compute_LMP2_residual_SFT(sma2::Array<4> &R);
    void compute_LMP2_residual_SFTS(sma2::Array<4> &R);
    void compute_LMP2_residual_SFTS_i(sma2::Array<4> &R);
    void compute_LMP2_residual_SFTS_ij(sma2::Array<4> &R);
 
    double iterate_LMP2_equations(double energy_tolerance, double rms_tolerance);
 
    void compute_delta_T(sma2::Array<4> &Res,
                         sma2::Array<4> &delta_T);
    
  public:
    LMP2(const sc::Ref<sc::KeyVal> &);
    LMP2(sc::StateIn &);
    ~LMP2();
    void save_data_state(sc::StateOut &);
    void compute(void);
    int nelectron(void);
    sc::RefSymmSCMatrix density(void);
    double magnetic_moment() const;
    int value_implemented(void) const;
 
};
 
}
 
#endif