psi4-data 1:1.1-5 / usr / share / psi4 / plugin / scf / scf.cc.template

/*
 * @BEGIN LICENSE
 *
 * @plugin@ by Psi4 Developer, a plugin to:
 *
 * Psi4: an open-source quantum chemistry software package
 *
 * Copyright (c) 2007-2017 The Psi4 Developers.
 *
 * The copyrights for code used from other parties are included in
 * the corresponding files.
 *
 * This file is part of Psi4.
 *
 * Psi4 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, version 3.
 *
 * Psi4 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 Psi4; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 *
 * @END LICENSE
 */

#include "scf.h"

#include "psi4/liboptions/liboptions.h"
#include "psi4/libfock/jk.h"
#include "psi4/libmints/integral.h"
#include "psi4/libmints/vector.h"
#include "psi4/libmints/molecule.h"
#include "psi4/libmints/basisset.h"

namespace psi{ namespace @plugin@{

SCF::SCF(SharedWavefunction ref_wfn, Options &options)
   : Wavefunction(options)
{

    // Shallow copy useful objects from the passed in wavefunction
    shallow_copy(ref_wfn);

    print_   = options_.get_int("PRINT");
    maxiter_ = options_.get_int("SCF_MAXITER");
    e_convergence_ = options_.get_double("E_CONVERGENCE");
    d_convergence_ = options_.get_double("D_CONVERGENCE");

    init_integrals();
}


SCF::~SCF()
{
}

void SCF::init_integrals()
{

    // The basisset object contains all of the basis information and is formed in the new_wavefunction call
    // The integral factory oversees the creation of integral objects
    std::shared_ptr<IntegralFactory> integral(new IntegralFactory
            (basisset_, basisset_, basisset_, basisset_));

    // Determine the number of electrons in the system
    // The molecule object is built into all wavefunctions
    int charge = molecule_->molecular_charge();
    int nelec  = 0;
    for(int i = 0; i < molecule_->natom(); ++i)
        nelec += (int)molecule_->Z(i);
    nelec -= charge;
    if(nelec % 2)
        throw PSIEXCEPTION("This is only an RHF code, but you gave it an odd number of electrons.  Try again!");
    ndocc_ = nelec / 2;

    outfile->Printf("    There are %d doubly occupied orbitals\n", ndocc_);
    molecule_->print();
    if(print_ > 1){
         basisset_->print_detail();
         options_.print();
    }

    nso_ =  basisset_->nbf();

    e_nuc_ = molecule_->nuclear_repulsion_energy();
    outfile->Printf("\n    Nuclear repulsion energy: %16.8f\n\n", e_nuc_);

    // These don't need to be declared, because they belong to the class
    S_ = SharedMatrix(new Matrix("Overlap matrix", nso_, nso_));
    H_ = SharedMatrix(new Matrix("Core Hamiltonian matrix", nso_, nso_));
    // These don't belong to the class, so we have to define them as having type SharedMatrix
    SharedMatrix T = SharedMatrix(new Matrix("Kinetic integrals matrix", nso_, nso_));
    SharedMatrix V = SharedMatrix(new Matrix("Potential integrals matrix", nso_, nso_));

    // Form the one-electron integral objects from the integral factory
    std::shared_ptr<OneBodyAOInt> sOBI(integral->ao_overlap());
    std::shared_ptr<OneBodyAOInt> tOBI(integral->ao_kinetic());
    std::shared_ptr<OneBodyAOInt> vOBI(integral->ao_potential());
    // Compute the one electron integrals, telling each object where to store the result
    sOBI->compute(S_);
    tOBI->compute(T);
    vOBI->compute(V);

    // Form h = T + V by first cloning T and then adding V
    H_->copy(T);
    H_->add(V);

    if(print_ > 3){
        S_->print();
        T->print();
        V->print();
        H_->print();
    }

    outfile->Printf("    Forming JK object\n\n");
    // Construct a JK object that compute J and K SCF matrices very efficiently
    jk_ = JK::build_JK(basisset_, std::shared_ptr<BasisSet>(), options_);

    // This is a very heavy compute object, lets give it 80% of our total memory
    jk_->set_memory(Process::environment.get_memory() * 0.8);
    jk_->initialize();
    jk_->print_header();

}


double SCF::compute_electronic_energy()
{
    SharedMatrix HplusF = H_->clone();
    HplusF->add(F_);
    return D_->vector_dot(HplusF);
}

void SCF::update_Cocc()
{
    for(int p = 0; p < nso_; ++p){
        for(int i = 0; i < ndocc_; ++i){
            Cocc_->set(p, i, C_->get(p, i));
        }
    }
}


double SCF::compute_energy()
{
   // Allocate some matrices
   X_  = SharedMatrix(new Matrix("S^-1/2", nso_, nso_));
   F_  = SharedMatrix(new Matrix("Fock matrix", nso_, nso_));
   Ft_ = SharedMatrix(new Matrix("Transformed Fock matrix", nso_, nso_));
   C_  = SharedMatrix(new Matrix("MO Coefficients", nso_, nso_));
   Cocc_ = SharedMatrix(new Matrix("Occupied MO Coefficients", nso_, ndocc_));
   D_  = SharedMatrix(new Matrix("The Density Matrix", nso_, nso_));
   SharedMatrix Temp1(new Matrix("Temporary Array 1", nso_, nso_));
   SharedMatrix Temp2(new Matrix("Temporary Array 2", nso_, nso_));
   SharedMatrix FDS(new Matrix("FDS", nso_, nso_));
   SharedMatrix SDF(new Matrix("SDF", nso_, nso_));
   SharedMatrix Evecs(new Matrix("Eigenvectors", nso_, nso_));
   SharedVector Evals(new Vector("Eigenvalues", nso_));

   // Form the X_ matrix (S^-1/2)
   X_->copy(S_);
   X_->power(-0.5, 1.e-14);

   F_->copy(H_);
   Ft_->transform(F_, X_);
   Ft_->diagonalize(Evecs, Evals, ascending);

   C_->gemm(false, false, 1.0, X_, Evecs, 0.0);
   update_Cocc();
   D_->gemm(false, true, 1.0, Cocc_, Cocc_, 0.0);

   if(print_ > 1){
       outfile->Printf("MO Coefficients and density from Core Hamiltonian guess:\n");
       C_->print();
       D_->print();
   }

   int iter = 1;
   bool converged = false;
   double e_old;
   double e_new = e_nuc_ + compute_electronic_energy();

   outfile->Printf("    Energy from core Hamiltonian guess: %20.16f\n\n", e_new);

   outfile->Printf("    *=======================================================*\n");
   outfile->Printf("    * Iter       Energy            delta E    ||gradient||  *\n");
   outfile->Printf("    *-------------------------------------------------------*\n");

   while(!converged && iter < maxiter_){
       e_old = e_new;

       // Add the core Hamiltonian term to the Fock operator
       F_->copy(H_);


       // The JK object handles all of the two electron integrals
       // To enhance efficiency it does use the density, but the orbitals themselves
       // D_uv = C_ui C_vj
       // J_uv = I_uvrs D_rs
       // K_uv = I_urvs D_rs

       // Here we clear the old Cocc and push_back our new one
       std::vector<SharedMatrix>& Cl = jk_->C_left();
       Cl.clear();
       Cl.push_back(Cocc_);
       jk_->compute();

       // Obtain the new J and K matrices
       const std::vector<SharedMatrix>& J = jk_->J();
       const std::vector<SharedMatrix>& K = jk_->K();

       // Proceede as normal
       J[0]->scale(2.0);
       F_->add(J[0]);
       F_->subtract(K[0]);

       // Compute the energy
       e_new = e_nuc_ + compute_electronic_energy();
       double dE = e_new - e_old;

       // Compute the orbital gradient, FDS-SDF
       Temp1->gemm(false, false, 1.0, D_, S_, 0.0);
       FDS->gemm(false, false, 1.0, F_, Temp1, 0.0);
       Temp1->gemm(false, false, 1.0, D_, F_, 0.0);
       SDF->gemm(false, false, 1.0, S_, Temp1, 0.0);
       Temp1->copy(FDS);
       Temp1->subtract(SDF);
       double dRMS = Temp1->rms();

       if(print_ > 1){
           Ft_->print();
           Evecs->print();
           Evals->print();
           C_->print();
           D_->print();
           FDS->print();
           SDF->print();
           Temp1->set_name("Orbital gradient");
           Temp1->print();
       }

       converged = (fabs(dE) < e_convergence_) && (dRMS < d_convergence_);

       outfile->Printf("    * %3d %20.14f    %9.2e    %9.2e    *\n", iter, e_new, dE, dRMS);

       // Transform the Fock operator and diagonalize it
       Ft_->transform(F_, X_);
       Ft_->diagonalize(Evecs, Evals, ascending);

       // Form the orbitals from the eigenvectors of the transformed Fock matrix
       C_->gemm(false, false, 1.0, X_, Evecs, 0.0);

       // Update our occupied orbitals
       update_Cocc();
       D_->gemm(false, true, 1.0, Cocc_, Cocc_, 0.0);
       iter++;

   }
   outfile->Printf("    *=======================================================*\n");

   if(!converged)
       throw PSIEXCEPTION("The SCF iterations did not converge.");

   Evals->set_name("Orbital Energies");
   Evals->print();
   energy_ = e_new;

   return e_new;
}

}} // End namespaces