mcpdft_solver.h

Go to the documentation of this file.

/*
 * @BEGIN LICENSE
 *
 * mcpdft by Psi4 Developer, a plugin to:
 *
 * Psi4: an open-source quantum chemistry software package
 *
 * Copyright (c) 2007-2016 The Psi4 Developers.
 *
 * The copyrights for code used from other parties are included in
 * the corresponding files.
 *
 * 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.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 * @END LICENSE
 */


#ifndef MCPDFT_SOLVER_H
#define MCPDFT_SOLVER_H

#define PSIF_DCC_QMO          268
#define PSIF_V2RDM_CHECKPOINT 269
#define PSIF_V2RDM_D2AA       270
#define PSIF_V2RDM_D2AB       271
#define PSIF_V2RDM_D2BB       272
#define PSIF_V2RDM_D3AAA      273
#define PSIF_V2RDM_D3AAB      274
#define PSIF_V2RDM_D3BBA      275
#define PSIF_V2RDM_D3BBB      276
#define PSIF_V2RDM_D1A        277
#define PSIF_V2RDM_D1B        278

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string>

#include "psi4/libmints/wavefunction.h"

// for reading integrals from disk
#include <psi4/libiwl/iwl.h>

// for dft
#include "psi4/libfock/v.h"
#include "psi4/libfunctional/superfunctional.h"

// for grid
#include "psi4/libfock/points.h"
#include "psi4/libfock/cubature.h"

#include "psi4/psi4-dec.h"
#include <psi4/psifiles.h>
#include <psi4/libpsio/psio.hpp>
#include <psi4/libpsi4util/PsiOutStream.h>

namespace psi{ namespace mcpdft{

struct tpdm {
    int i;
    int j;
    int k;
    int l;
    double val;
};

struct opdm {
    int i;
    int j;
    double val;
};

class MCPDFTSolver: public Wavefunction{

  public:

    MCPDFTSolver(std::shared_ptr<psi::Wavefunction> reference_wavefunction,Options & options);
    ~MCPDFTSolver();
    void common_init();
    double compute_energy();
    void polyradical_analysis();
    virtual bool same_a_b_orbs() const { return same_a_b_orbs_; }
    virtual bool same_a_b_dens() const { return same_a_b_dens_; }

    void print_banner();

//  protected:

    double coulomb_energy_;
    double hf_ex_energy_;
    double lr_ex_energy_;


    bool is_low_memory_;

    opdm * opdm_a_;

    opdm * opdm_b_;

    std::shared_ptr<VBase> potential_;

    std::shared_ptr<PointFunctions> points_func_;

    int max_functions_;

    int max_points_;

    int * symmetry_;

    std::vector < std::vector < std::pair<int,int> > > gems_;

    // geminal / orbital maps

    int *** bas_;
    int *** ibas_;

    // read erfc integrals from disk
    void ReadAllIntegrals(iwlbuf *Buf);

    void ReadRangeSeparatedIntegrals();

    double * erfc_tei_;

    double reference_energy_;

    double two_electron_energy_;

    double lr_Vee_energy_;

    double sr_Vee_energy_;

    double mp2_corr_energy_;

    int * pitzer_offset_;

    int deriv_;

    bool is_gga_;

    bool is_df_;

    bool is_meta_;

    long int available_memory_;

    bool is_unpolarized_;

    long int phi_points_;

    int nQ_;

//     std::shared_ptr<Matrix> super_phi_ao_;
// 
//     /// d phi / dx matrix (AO)
//     std::shared_ptr<Matrix> super_phi_x_ao_;
// 
//     /// d phi / dy matrix (AO)
//     std::shared_ptr<Matrix> super_phi_y_ao_;
// 
//     /// d phi / dz matrix (AO)
//     std::shared_ptr<Matrix> super_phi_z_ao_;

    std::shared_ptr<Matrix> super_phi_;

    std::shared_ptr<Matrix> super_phi_x_;

    std::shared_ptr<Matrix> super_phi_y_;

    std::shared_ptr<Matrix> super_phi_z_;

    std::shared_ptr<Matrix> super_phi_xx_;

    std::shared_ptr<Matrix> super_phi_xy_;

    std::shared_ptr<Matrix> super_phi_xz_;

    std::shared_ptr<Matrix> super_phi_yy_;

    std::shared_ptr<Matrix> super_phi_yz_;

    std::shared_ptr<Matrix> super_phi_zz_;

    std::shared_ptr<Matrix> super_gamma_aa_; 

    std::shared_ptr<Matrix> super_gamma_ab_;   

    std::shared_ptr<Matrix> super_gamma_bb_;   

    std::shared_ptr<Matrix> super_tau_a_;

    std::shared_ptr<Matrix> super_tau_b_;   

    std::shared_ptr<Vector> grid_x_;

    std::shared_ptr<Vector> grid_y_;

    std::shared_ptr<Vector> grid_z_;

    std::shared_ptr<Vector> grid_w_;
    
    std::shared_ptr<Vector> sigma_aa_;

    std::shared_ptr<Vector> sigma_bb_;

    std::shared_ptr<Vector> sigma_ab_;

    void GetGridInfo();

    void BuildPhiMatrixAO(std::string phi_type, std::shared_ptr<Matrix> myphi);

    void TransformPhiMatrixAOMO(std::shared_ptr<Matrix> phi_in, std::shared_ptr<Matrix> phi_out);

    void PrintTPDM(double* D);

    void ReadTPDM();

    void ReadOPDM();

    void ReadCITPDM(double* D, const char* fileName); 

    void ReadCIOPDM(std::shared_ptr<Matrix> D, const char* fileName); 

    std::vector< std::shared_ptr<Matrix> > BuildJK();

    double RangeSeparatedTEE(std::string range_separation_type);

    std::shared_ptr<Vector> rho_a_;

    std::shared_ptr<Vector> rho_b_;

    std::shared_ptr<Vector> rho_;
 
    std::shared_ptr<Vector> rho_a_x_;

    std::shared_ptr<Vector> rho_b_x_;

    std::shared_ptr<Vector> rho_a_y_;

    std::shared_ptr<Vector> rho_b_y_;

    std::shared_ptr<Vector> rho_a_z_;

    std::shared_ptr<Vector> rho_b_z_;
    
    std::shared_ptr<Vector> tr_rho_a_;

    std::shared_ptr<Vector> tr_rho_b_;

    std::shared_ptr<Vector> tr_rho_a_x_;

    std::shared_ptr<Vector> tr_rho_b_x_;

    std::shared_ptr<Vector> tr_rho_a_y_;

    std::shared_ptr<Vector> tr_rho_b_y_;

    std::shared_ptr<Vector> tr_rho_a_z_;

    std::shared_ptr<Vector> tr_rho_b_z_;
    
    std::shared_ptr<Vector> tr_sigma_aa_;

    std::shared_ptr<Vector> tr_sigma_bb_;

    std::shared_ptr<Vector> tr_sigma_ab_;
    
    std::shared_ptr<Vector> tr_sigma_;

    std::shared_ptr<Vector> R_;

    std::shared_ptr<Vector> zeta_;
   
    std::shared_ptr<Vector> m_;
    
    std::shared_ptr<Vector> tr_m_;

    std::shared_ptr<Vector> tr_zeta_;

    std::shared_ptr<Vector> rs_;

    std::shared_ptr<Vector> tr_rs_;

    std::shared_ptr<Vector> tau_a_;

    std::shared_ptr<Vector> tau_b_;

    std::shared_ptr<Vector> pi_;

    std::shared_ptr<Vector> pi_x_;

    std::shared_ptr<Vector> pi_y_;

    std::shared_ptr<Vector> pi_z_;

    void BuildRho();

    void BuildRhoFast(int na, int nb);

    void BuildPi(double * D2ab);
 
    void BuildPiFast(tpdm * D2ab, int nab);

    void BuildPiLowMemory(tpdm * D2ab, int nab);

    void BuildExchangeCorrelationHole(int p, tpdm * D2ab, int nab, tpdm * D2aa, int naa, tpdm * D2bb, int nbb);

    std::shared_ptr<Vector> Dr_;

    std::shared_ptr<Vector> Ur_;

    std::shared_ptr<Vector> Nr_;

    void Build_Grad_Pi();

    void Build_R();
    
    void Translate();

    void Fully_Translate();

    double Gfunction(double r, double A, double a1, double b1, double b2, double b3, double b4, double p);

    //############################################################
    //############ Exchange functions' declarations ##############
    //############################################################

    double EX_LDA(std::shared_ptr<Vector> rho_a, std::shared_ptr<Vector> rho_b);

    double EX_LSDA_Sigma(std::shared_ptr<Vector> rho_sigma);

    double EX_LSDA(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> ZETA);

    double EX_LSDA(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B);

    double EX_B86_MGC();

    double EX_B88(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
    double EX_B88_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);

    double EX_PBE(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);
    double EX_PBE_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);

    double EX_wPBE_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);

    double EX_RPBE();

    double EX_UPBE();

    //############################################################
    //########### Correlation functions' declarations ############
    //############################################################

    double EC_B88_OP(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_BB);

    double EC_LYP_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, 
                  std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_AB, std::shared_ptr<Vector> SIGMA_BB);

    double EC_PBE(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, 
                  std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_AB, std::shared_ptr<Vector> SIGMA_BB);
    double EC_PBE_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, 
                  std::shared_ptr<Vector> SIGMA_AA, std::shared_ptr<Vector> SIGMA_AB, std::shared_ptr<Vector> SIGMA_BB);

    double EC_VWN3_RPA(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> ZETA, std::shared_ptr<Vector> RS);
    double EC_VWN3_RPA_III(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B);
    // double EC_VWN3_RPA_III(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B, std::shared_ptr<Vector> ZETTA);

    double EC_PW92_I(std::shared_ptr<Vector> RHO_A, std::shared_ptr<Vector> RHO_B);

};

}} // end of namespaces

#endif