chem/SCF.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 MADNESS_CHEM_SCF_H__INCLUDED
#define MADNESS_CHEM_SCF_H__INCLUDED

//#define WORLD_INSTANTIATE_STATIC_TEMPLATES

#include <madness/mra/mra.h>

#include <chem/molecule.h>
#include <chem/molecularbasis.h>
#include <chem/corepotential.h>
#include <chem/xcfunctional.h>
#include <chem/potentialmanager.h>
#include <chem/gth_pseudopotential.h> 

#include <madness/tensor/solvers.h>
#include <madness/tensor/distributed_matrix.h>


namespace madness {

typedef std::shared_ptr< WorldDCPmapInterface< Key<3> > > pmapT;
typedef Vector<double,3> coordT;
typedef std::shared_ptr< FunctionFunctorInterface<double,3> > functorT;
typedef Function<double,3> functionT;
typedef std::vector<functionT> vecfuncT;
typedef std::pair<vecfuncT,vecfuncT> pairvecfuncT;
typedef std::vector<pairvecfuncT> subspaceT;
typedef Tensor<double> tensorT;
typedef DistributedMatrix<double> distmatT;
typedef FunctionFactory<double,3> factoryT;
typedef SeparatedConvolution<double,3> operatorT;
typedef std::shared_ptr<operatorT> poperatorT;
typedef Function<std::complex<double>,3> complex_functionT;
typedef std::vector<complex_functionT> cvecfuncT;
typedef Convolution1D<double_complex> complex_operatorT;

    extern distmatT distributed_localize_PM(World & world,
                                const vecfuncT & mo,
                                const vecfuncT & ao,
                                const std::vector<int> & set,
                                const std::vector<int> & at_to_bf,
                                const std::vector<int> & at_nbf,
                                const double thresh = 1e-9,
                                const double thetamax = 0.5,
                                const bool randomize = true,
                                       const bool doprint = false);

inline double mask1(double x) {
    /* Iterated first beta function to switch smoothly
       from 0->1 in [0,1].  n iterations produce 2*n-1
       zero derivatives at the end points. Order of polyn
       is 3^n.

       Currently use one iteration so that first deriv.
       is zero at interior boundary and is exactly representable
       by low order multiwavelet without refinement */

    x = (x*x*(3.-2.*x));
    return x;
}

static double mask3(const coordT& ruser) {
    coordT rsim;
    user_to_sim(ruser, rsim);
    double x= rsim[0], y=rsim[1], z=rsim[2];
    double lo = 0.0625, hi = 1.0-lo, result = 1.0;
    double rlo = 1.0/lo;

    if (x<lo)
        result *= mask1(x*rlo);
    else if (x>hi)
        result *= mask1((1.0-x)*rlo);
    if (y<lo)
        result *= mask1(y*rlo);
    else if (y>hi)
        result *= mask1((1.0-y)*rlo);
    if (z<lo)
        result *= mask1(z*rlo);
    else if (z>hi)
        result *= mask1((1.0-z)*rlo);

    return result;
}

class MolecularGuessDensityFunctor : public FunctionFunctorInterface<double,3> {
private:
    const Molecule& molecule;
    const AtomicBasisSet& aobasis;
public:
    MolecularGuessDensityFunctor(const Molecule& molecule, const AtomicBasisSet& aobasis)
        : molecule(molecule), aobasis(aobasis) {}

    double operator()(const coordT& x) const {
        return aobasis.eval_guess_density(molecule, x[0], x[1], x[2]);
    }

    std::vector<coordT> special_points() const {return molecule.get_all_coords_vec();}
};


class AtomicBasisFunctor : public FunctionFunctorInterface<double,3> {
private:
    const AtomicBasisFunction aofunc;

public:
    AtomicBasisFunctor(const AtomicBasisFunction& aofunc)
        : aofunc(aofunc)
    {}

    double operator()(const coordT& x) const {
        return aofunc(x[0], x[1], x[2]);
    }

    std::vector<coordT> special_points() const {
        return std::vector<coordT>(1,aofunc.get_coords_vec());
    }
};

class MolecularDerivativeFunctor : public FunctionFunctorInterface<double,3> {
private:
    const Molecule& molecule;
    const int atom;
    const int axis;

public:
    MolecularDerivativeFunctor(const Molecule& molecule, int atom, int axis)
        : molecule(molecule), atom(atom), axis(axis)
    {}

    double operator()(const coordT& x) const {
        return molecule.nuclear_attraction_potential_derivative(atom, axis, x[0], x[1], x[2]);
    }

    std::vector<coordT> special_points() const {
        return std::vector<coordT>(1,molecule.get_atom(atom).get_coords());
    }
};

class CorePotentialDerivativeFunctor : public FunctionFunctorInterface<double,3> {
private:
    const Molecule& molecule;
    const int atom;
    const int axis;
    std::vector<coordT> specialpt;
public:
    CorePotentialDerivativeFunctor(const Molecule& molecule, int atom, int axis)
        : molecule(molecule), atom(atom), axis(axis) {}

    double operator()(const coordT& r) const {
        return molecule.core_potential_derivative(atom, axis, r[0], r[1], r[2]);
    }
};

class DipoleFunctor : public FunctionFunctorInterface<double,3> {
private:
    const int axis;
public:
    DipoleFunctor(int axis) : axis(axis) {}
    double operator()(const coordT& x) const {
        return x[axis];
    }
};



class MomentFunctor : public FunctionFunctorInterface<double,3> {
private:
    const int i, j, k;
public:
    MomentFunctor(int i, int j, int k) : i(i), j(j), k(k) {}
    MomentFunctor(const std::vector<int>& x) : i(x[0]), j(x[1]), k(x[2]) {}
    double operator()(const coordT& r) const {
        double xi=1.0, yj=1.0, zk=1.0;
        for (int p=0; p<i; ++p) xi *= r[0];
        for (int p=0; p<j; ++p) yj *= r[1];
        for (int p=0; p<k; ++p) zk *= r[2];
        return xi*yj*zk;
    }
};

template<typename T>
class VextCosFunctor {
  double _omega;
  Function<T,3> _f;
public:
    VextCosFunctor(World& world,
//        const std::shared_ptr<FunctionFunctorInterface<T,3> >& functor,
        const FunctionFunctorInterface<T,3>* functor,
        double omega) : _omega(omega)
    {
//      _f = factoryT(world).functor(functor);
      _f = factoryT(world).functor(functorT(new DipoleFunctor(2)));
    }
    Function<T,3> operator()(const double t) const {
        return std::cos(_omega * t) * _f;
    }
};


template <typename T, int NDIM>
struct lbcost {
    double leaf_value;
    double parent_value;
    lbcost(double leaf_value=1.0, double parent_value=0.0) : leaf_value(leaf_value), parent_value(parent_value) {}
    double operator()(const Key<NDIM>& key, const FunctionNode<T,NDIM>& node) const {
        if (key.level() < 1) {
            return 100.0*(leaf_value+parent_value);
        }
        else if (node.is_leaf()) {
            return leaf_value;
        }
        else {
            return parent_value;
        }
    }
};


struct CalculationParameters {
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    // !!!                                                                   !!!
    // !!! If you add more data don't forget to add them to serialize method !!!
    // !!!                                                                   !!!
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    // First list input parameters
    double charge;              
    double smear;               
    double econv;               
    double dconv;               
    int k;                                              
    double L;                   
    double maxrotn;             
    int nvalpha;                
    int nvbeta;                 
    int nopen;                  
    int maxiter;                
    int nio;                    
    bool spin_restricted;       
    int plotlo,plothi;          
    bool plotdens;              
    bool plotcoul;              
    bool localize;              
    bool localize_pm;           
    bool restart;               
    bool no_compute;            
    bool no_orient;                             
    bool save;                  
    unsigned int maxsub;        
    double orbitalshift;                
    int npt_plot;               
    tensorT plot_cell;          
    std::string aobasis;        
    std::string core_type;      
    bool derivatives;           
    bool dipole;                
    bool conv_only_dens;        
    // Next list inferred parameters
    int nalpha;                 
    int nbeta;                  
    int nmo_alpha;              
    int nmo_beta;               
    double lo;                  
    std::string xc_data;         
    std::vector<double> protocol_data;  
    bool gopt;                  
    double gtol;                
    bool gtest;                 
    double gval;                
    double gprec;               
    int  gmaxiter;               
    std::string algopt;         
    bool tdksprop;               
    std::string nuclear_corrfac;        
    bool psp_calc;                
    //bool absolvent;             ///< If true calculate solvation effects
    //double epsilon_2;           ///< dielectric constant of solvent
    //double Gamma;               ///< surface tension of solvent
    //double beta;                ///switching parameter controles boundary conditions of solvent cavity
    //double rho_0;               /// threshold density--determines size of molecular cavity
    //double sigma;               ///switching parameter controles boundary conditions of solvent cavity-(SVPE)

    template <typename Archive>
    void serialize(Archive& ar) {
        ar & charge & smear & econv & dconv & k & L & maxrotn & nvalpha & nvbeta
           & nopen & maxiter & nio & spin_restricted;
        ar & plotlo & plothi & plotdens & plotcoul & localize & localize_pm
           & restart & save & no_compute &no_orient & maxsub & orbitalshift & npt_plot & plot_cell & aobasis;
        ar & nalpha & nbeta & nmo_alpha & nmo_beta & lo;
        ar & core_type & derivatives & conv_only_dens & dipole;
        ar & xc_data & protocol_data;
        ar & gopt & gtol & gtest & gval & gprec & gmaxiter & algopt & tdksprop & psp_calc;
        //ar & absolvent & epsilon_2 & Gamma & Gamma & beta & rho_0 & sigma;
    }

    CalculationParameters()
        : charge(0.0)
        , smear(0.0)
        , econv(1e-5)
        , dconv(1e-4)
        , k(-1)
        , L(0.0)
        , maxrotn(0.25)
        , nvalpha(0)
        , nvbeta(0)
        , nopen(0)
        , maxiter(20)
        , nio(1)
        , spin_restricted(true)
        , plotlo(0)
        , plothi(-1)
        , plotdens(false)
        , plotcoul(false)
        , localize(true)
        , localize_pm(true)
        , restart(false)
        , no_compute(false)
        , no_orient(false)
        , save(true)
        , maxsub(8)
        , orbitalshift(0.0)
        , npt_plot(101)
        , aobasis("6-31g")
        , core_type("")
        , derivatives(false)
        , dipole(false)
        , conv_only_dens(false)
        , nalpha(0)
        , nbeta(0)
        , nmo_alpha(0)
        , nmo_beta(0)
        , lo(1e-10)
        , xc_data("lda")
        , protocol_data(madness::vector_factory(1e-4, 1e-6))
        , gopt(false)
        , gtol(1e-3)
        , gtest(false)
        , gval(1e-5)
        , gprec(1e-4)
        , gmaxiter(20)
        , algopt("BFGS")
        , tdksprop(false)
        , nuclear_corrfac("none")
        , psp_calc(false)
          //, absolvent(false)
          //, epsilon_2(78.304)
          //, Gamma(0.0719)
          //, beta(1.3)
          //, rho_0(0.00048)
          //, sigma(0.3)
    {}


    void read_file(const std::string& filename) {
        std::ifstream f(filename.c_str());
        position_stream(f, "dft");
        std::string s;
        xc_data = "lda";
        protocol_data = madness::vector_factory(1e-4, 1e-6);

        while (f >> s) {
            if (s == "end") {
                break;
            }
            else if (s == "charge") {
                f >> charge;
            }
            else if (s == "smear") {
                f >> smear;
            }
            else if (s == "econv") {
                f >> econv;
            }
            else if (s == "dconv") {
                f >> dconv;
            }
            else if (s == "k") {
                f >> k;
            }
            else if (s == "L") {
                f >> L;
            }
            else if (s == "maxrotn") {
                f >> maxrotn;
            }
            else if (s == "nvalpha") {
                f >> nvalpha;
            }
            else if (s == "nvbeta") {
                f >> nvbeta;
            }
            else if (s == "nopen") {
                f >> nopen;
            }
            else if (s == "unrestricted") {
                spin_restricted = false;
            }
            else if (s == "restricted") {
                spin_restricted = true;
            }
            else if (s == "maxiter") {
                f >> maxiter;
            }
            else if (s == "nio") {
                f >> nio;
            }
            else if (s == "xc") {
                char buf[1024];
                f.getline(buf,sizeof(buf));
                xc_data = buf;
            }
            else if (s == "protocol") {
                std::string buf;
                std::getline(f,buf);
                protocol_data = std::vector<double>();
                double d;
                std::stringstream s(buf);
                while (s >> d) protocol_data.push_back(d);
            }
            else if (s == "plotmos") {
                f >> plotlo >> plothi;
            }
            else if (s == "plotdens") {
                plotdens = true;
            }
            //            else if (s == "absolvent") {
            //                absolvent = true;
            //            }
            //            else if (s == "dielec") {
            //                f >> epsilon_2;
            //            }
            //            else if (s == "Gamma") {
            //                f >> Gamma;
            //            }
            //            else if (s == "beta") {
            //                f >> beta;
            //            }
            //            else if (s == "rho_0") {
            //                f >> rho_0;
            //            }
            //            else if (s == "sigma") {
            //                f >> sigma;
            //            }
            else if (s == "plotcoul") {
                plotcoul = true;
            }
            else if (s == "plotnpt") {
                f >> npt_plot;
            }
            else if (s == "plotcell") {
                plot_cell = tensorT(3L,2L);
                f >> plot_cell(0,0) >> plot_cell(0,1) >> plot_cell(1,0) >> plot_cell(1,1) >> plot_cell(2,0) >> plot_cell(2,1);
            }
            else if (s == "aobasis") {
                f >> aobasis;
                if (aobasis!="sto-3g" && aobasis!="sto-6g" && aobasis!="6-31g") {
                    std::cout << "moldft: unrecognized aobasis (sto-3g or sto-6g or 6-31g only): " << aobasis << std::endl;
                    MADNESS_EXCEPTION("input_error", 0);
                }
            }
            else if (s == "canon") {
                localize = false;
            }
            else if (s == "local") {
                localize = true;
            }
            else if (s == "pm") {
                localize_pm = true;
            }
            else if (s == "boys") {
                localize_pm = false;
            }
            else if (s == "restart") {
                restart = true;
            }
            else if (s == "save") {
                //can't redirect true/false as with other variables so create temporary variable
                std::string tmp_save;
                f >> tmp_save;
                if (tmp_save=="true"){
                    save=true;
                }
                else if (tmp_save=="false"){
                    save=false;
                }
                else{
                    std::cout << "moldft: unrecognized value for save (true or false only): " << tmp_save << std::endl;
                    MADNESS_EXCEPTION("input_error", 0);
                }
            }
            else if (s == "no_compute") {
                no_compute = true;
            }
            else if (s == "no_orient") {
                no_orient = true;
            }
            else if (s == "maxsub") {
                f >> maxsub;
                if (maxsub <= 0) maxsub = 1;
                if (maxsub > 20) maxsub = 20;
            }
            else if (s == "orbitalshift") {
                f >> orbitalshift;
            }
            else if (s == "core_type") {
                f >> core_type;
            }
            else if (s == "derivatives") {
                derivatives = true;
            }
            else if (s == "dipole") {
                dipole = true;
            }
            else if (s == "convonlydens") {
                conv_only_dens = true;
            }
            else if (s == "gopt") {
               gopt = true;
            }
            else if (s == "gtol") {
                f >> gtol;
            }
            else if (s == "gtest") {
               gtest = true;
            }
            else if (s == "gval") {
                f >> gval;
            }
            else if (s == "gprec") {
                f >> gprec;
            }
            else if (s == "gmaxiter") {
                f >> gmaxiter;
            }
            else if (s == "algopt") {
                char buf[1024];
                f.getline(buf,sizeof(buf));
                algopt = buf;
            }
            else if (s == "tdksprop") {
              tdksprop = true;
            }
            else if (s == "nuclear_corrfac") {
                std::string str;
                std::getline(f,str);
                nuclear_corrfac=str;
            }
            else if (s == "psp_calc") {
              psp_calc = true;
            }
            else {
                std::cout << "moldft: unrecognized input keyword " << s << std::endl;
                MADNESS_EXCEPTION("input error",0);
            }
            if (nopen != 0) spin_restricted = false;
        }
    }

    void set_molecular_info(const Molecule& molecule, const AtomicBasisSet& aobasis, unsigned int n_core) {
        double z = molecule.total_nuclear_charge();
        int nelec = int(z - charge - n_core*2);
        if (fabs(nelec+charge+n_core*2-z) > 1e-6) {
            error("non-integer number of electrons?", nelec+charge+n_core*2-z);
        }
        nalpha = (nelec + nopen)/2;
        nbeta  = (nelec - nopen)/2;
        if (nalpha < 0) error("negative number of alpha electrons?", nalpha);
        if (nbeta < 0) error("negative number of beta electrons?", nbeta);
        if ((nalpha+nbeta) != nelec) error("nalpha+nbeta != nelec", nalpha+nbeta);
        nmo_alpha = nalpha + nvalpha;
        nmo_beta = nbeta + nvbeta;
        if (nalpha != nbeta) spin_restricted = false;

        // Ensure we have enough basis functions to guess the requested
        // number of states ... a minimal basis for a closed-shell atom
        // might not have any functions for virtuals.
        int nbf = aobasis.nbf(molecule);
        nmo_alpha = std::min(nbf,nmo_alpha);
        nmo_beta = std::min(nbf,nmo_beta);
        if (nalpha>nbf || nbeta>nbf) error("too few basis functions?", nbf);
        nvalpha = nmo_alpha - nalpha;
        nvbeta = nmo_beta - nbeta;

        // Unless overridden by the user use a cell big enough to
        // have exp(-sqrt(2*I)*r) decay to 1e-6 with I=1ev=0.037Eh
        // --> need 50 a.u. either side of the molecule

        if (L == 0.0) {
            L = molecule.bounding_cube() + 50.0;
        }

        lo = molecule.smallest_length_scale();
    }

    void print(World& world) const {
        //time_t t = time((time_t *) 0);
        //char *tmp = ctime(&t);
        //tmp[strlen(tmp)-1] = 0; // lose the trailing newline

        //madness::print(" date of calculation ", tmp);
        madness::print("             restart ", restart);
        madness::print(" number of processes ", world.size());
        madness::print("   no. of io servers ", nio);
        madness::print("     simulation cube ", -L, L);
        madness::print("        total charge ", charge);
        madness::print("            smearing ", smear);
        madness::print(" number of electrons ", nalpha, nbeta);
        madness::print("  number of orbitals ", nmo_alpha, nmo_beta);
        madness::print("     spin restricted ", spin_restricted);
        madness::print("       xc functional ", xc_data);
#ifdef MADNESS_HAS_LIBXC
        madness::print("         xc library  ", "libxc");
#else
        madness::print("         xc library  ", "default (lda only)");
#endif
        if (core_type != "")
            madness::print("           core type ", core_type);
        madness::print(" initial guess basis ", aobasis);
        madness::print(" max krylov subspace ", maxsub);
        madness::print("    compute protocol ", protocol_data);
        madness::print("  energy convergence ", econv);
        madness::print(" density convergence ", dconv);
        madness::print("    maximum rotation ", maxrotn);
        madness::print("    polynomial order ", k);
        madness::print("  maximum iterations ", maxiter);
        if (conv_only_dens)
            madness::print(" Convergence criterion is only density delta.");
        else
            madness::print(" Convergence criteria are density delta & BSH residual.");
        madness::print("        plot density ", plotdens);
        madness::print("        plot coulomb ", plotcoul);
        madness::print("        plot orbital ", plotlo, plothi);
        madness::print("        plot npoints ", npt_plot);
        if (plot_cell.size() > 0)
            madness::print("        plot  volume ", plot_cell(0,0), plot_cell(0,1),
                           plot_cell(1,0), plot_cell(1,1), plot_cell(2,0), plot_cell(2,1));
        else
            madness::print("        plot  volume ", "default");

        std::string loctype = "pm";
        if (!localize_pm) loctype = "boys";
        if (localize)
            madness::print("  localized orbitals ", loctype);
        else
            madness::print("  canonical orbitals ");
        if (derivatives)
            madness::print("    calc derivatives ");
        if (dipole)
            madness::print("         calc dipole ");
        //        if(absolvent){
        //            madness::print("       isodensity solvation ", absolvent);
        //            madness::print("       surface tension      ", Gamma);
        //            madness::print("       switching param(beta)", beta);
        //            madness::print("       dielectric constant  ", epsilon_2);
        //            madness::print("       threshold density    ", rho_0);
        //        }
    }
//};
    void gprint(World& world) const {
        madness::print(" Optimizer parameters:           ");
        madness::print("   Maximum iterations (gmaxiter) ", gmaxiter);
        madness::print("                Tolerance (gtol) ", gtol);
        madness::print("           Gradient value (gval) ", gval);
        madness::print("      Gradient precision (gprec) ", gprec);
        madness::print(" Optimization algorithm (algopt) ", algopt);
        madness::print(" Gradient numerical test (gtest) ", gtest);
    }
};

class SCF {
public:
    std::shared_ptr<PotentialManager> potentialmanager;
    std::shared_ptr<GTHPseudopotential <double> > gthpseudopotential;
    Molecule molecule;
    CalculationParameters param;
    XCfunctional xc;
    AtomicBasisSet aobasis;
    //functionT vacuo_rho;
    //functionT rhoT;
    //functionT rho_elec;
    //functionT rhon;
    //functionT mol_mask;
    //functionT Uabinit;
    functionT mask;

    vecfuncT amo, bmo;

    std::vector<int> aset, bset;

    vecfuncT ao;


    std::vector<int> at_to_bf, at_nbf;

    tensorT aocc, bocc;

    tensorT aeps, beps;
    poperatorT coulop;
    std::vector< std::shared_ptr<real_derivative_3d> > gradop;
    double vtol;
    double current_energy;
    //double esol;//etot;
    //double vacuo_energy;
    static const int vnucextra = 12; // load balance parameter for nuclear pot.

    SCF(World & world, const char *filename);

    template<std::size_t NDIM>
    void set_protocol(World & world, double thresh)
    {
        int k;
        // Allow for imprecise conversion of threshold
        if(thresh >= 0.9e-2)
            k = 4;
        else if(thresh >= 0.9e-4)
            k = 6;
        else if(thresh >= 0.9e-6)
            k = 8;
        else if(thresh >= 0.9e-8)
            k = 10;
        else
            k = 12;

        // k defaults to make sense with thresh, override by providing k in input file
        if (param.k == -1) {
                FunctionDefaults<NDIM>::set_k(k);
//              param.k=k;
        } else {
                FunctionDefaults<NDIM>::set_k(param.k);
        }
        // don't forget to adapt the molecular smoothing parameter!!
        molecule.set_eprec(std::min(thresh,molecule.get_eprec()));
        FunctionDefaults<NDIM>::set_thresh(thresh);
        FunctionDefaults<NDIM>::set_refine(true);
        FunctionDefaults<NDIM>::set_initial_level(2);
        FunctionDefaults<NDIM>::set_truncate_mode(1);
        FunctionDefaults<NDIM>::set_autorefine(false);
        FunctionDefaults<NDIM>::set_apply_randomize(false);
        FunctionDefaults<NDIM>::set_project_randomize(false);
        FunctionDefaults<NDIM>::set_cubic_cell(-param.L, param.L);
        GaussianConvolution1DCache<double>::map.clear();
        double safety = 0.1;
        vtol = FunctionDefaults<NDIM>::get_thresh() * safety;
        coulop = poperatorT(CoulombOperatorPtr(world, param.lo, thresh));
        gradop = gradient_operator<double,3>(world);
        mask = functionT(factoryT(world).f(mask3).initial_level(4).norefine());
        if(world.rank() == 0){
            print("\nSolving NDIM=",NDIM," with thresh", thresh, "    k",
                FunctionDefaults<NDIM>::get_k(), "   conv", std::max(thresh, param.dconv), "\n");
        }
    }

    void save_mos(World& world);

    void load_mos(World& world);

    void do_plots(World& world);

    void project(World & world);

    void make_nuclear_potential(World & world);

    void project_ao_basis(World & world);


    std::vector<int> group_orbital_sets(World& world, const tensorT& eps,
                const tensorT& occ, const int nmo) const;


    distmatT localize_PM(World & world, const vecfuncT & mo, const std::vector<int> & set,
                const double thresh = 1e-9, const double thetamax = 0.5,
                const bool randomize = true, const bool doprint = false) const;


    void analyze_vectors(World & world, const vecfuncT & mo, const tensorT & occ = tensorT(),
                const tensorT & energy = tensorT(), const std::vector<int> & set = std::vector<int>());

    inline double DIP(const tensorT & dip, int i, int j, int k, int l) {
        return dip(i, j, 0) * dip(k, l, 0) + dip(i, j, 1) * dip(k, l, 1) + dip(i, j, 2) * dip(k, l, 2);
    }

    // tensorT localize_boys(World & world, const vecfuncT & mo, const std::vector<int> & set,
    //          const double thresh = 1e-9, const double thetamax = 0.5, const bool randomize = true);


    distmatT kinetic_energy_matrix(World & world, const vecfuncT & v) const;


    vecfuncT core_projection(World & world, const vecfuncT & psi, const bool include_Bc = true);

    double core_projector_derivative(World & world, const vecfuncT & mo,
                        const tensorT & occ, int atom, int axis);

    void initial_guess(World & world);

    void initial_load_bal(World & world);

    functionT make_density(World & world, const tensorT & occ, const vecfuncT & v) const;

    functionT make_density(World & world, const tensorT & occ, const cvecfuncT & v);

    std::vector<poperatorT> make_bsh_operators(World & world, const tensorT & evals);


    vecfuncT apply_hf_exchange(World & world, const tensorT & occ,
                const vecfuncT & psi, const vecfuncT & f) const ;

    // Used only for initial guess that is always spin-restricted LDA
    functionT make_lda_potential(World & world, const functionT & arho);


    functionT make_dft_potential(World & world, const vecfuncT& vf, int ispin, int what)
    {
        return multiop_values<double, xc_potential, 3>(xc_potential(xc, ispin, what), vf);
    }

    double make_dft_energy(World & world, const vecfuncT& vf, int ispin)
    {
        functionT vlda = multiop_values<double, xc_functional, 3>(xc_functional(xc, ispin), vf);
        return vlda.trace();
    }

        vecfuncT apply_potential(World & world, const tensorT & occ,
                        const vecfuncT & amo, const vecfuncT& vf, const vecfuncT& delrho,
                        const functionT & vlocal, double & exc, double & enl, int ispin);

    tensorT derivatives(World & world);

    tensorT dipole(World & world);

    void vector_stats(const std::vector<double> & v, double & rms,
                double & maxabsval) const;

        vecfuncT compute_residual(World & world, tensorT & occ, tensorT & fock,
                        const vecfuncT & psi, vecfuncT & Vpsi, double & err);

        tensorT make_fock_matrix(World & world, const vecfuncT & psi,
                        const vecfuncT & Vpsi, const tensorT & occ,
                        double & ekinetic) const;


    functionT make_coulomb_potential(const functionT& rho) const {
        return apply(*coulop, rho);
    }


    Tensor<double> twoint(World& world, const vecfuncT& psi) const;

    tensorT matrix_exponential(const tensorT& A) const ;


    tensorT get_fock_transformation(World& world, const tensorT& overlap,
                tensorT& fock, tensorT& evals, const tensorT& occ,
                const double thresh_degenerate) const;



    tensorT diag_fock_matrix(World& world, tensorT& fock,
                vecfuncT& psi, vecfuncT& Vpsi, tensorT& evals,
                const tensorT& occ, const double thresh) const;


    void loadbal(World & world, functionT & arho, functionT & brho, functionT & arho_old,
                functionT & brho_old, subspaceT & subspace);


    void rotate_subspace(World& world, const tensorT& U, subspaceT& subspace,
                int lo, int nfunc, double trantol) const;

    void rotate_subspace(World& world, const distmatT& U, subspaceT& subspace,
                int lo, int nfunc, double trantol) const;

    void update_subspace(World & world,
                         vecfuncT & Vpsia, vecfuncT & Vpsib,
                         tensorT & focka, tensorT & fockb,
                         subspaceT & subspace, tensorT & Q,
                         double & bsh_residual, double & update_residual);


    double do_step_restriction(World& world, const vecfuncT& mo,
                vecfuncT& mo_new, std::string spin) const;



    void orthonormalize(World& world, vecfuncT& amo_new) const;


    void propagate(World& world, double omega, int step0);


    complex_functionT APPLY(const complex_operatorT* q1d, const complex_functionT& psi) {
        complex_functionT r = psi;  // Shallow copy violates constness !!!!!!!!!!!!!!!!!
        coordT lo, hi;
        lo[2] = -10;
        hi[2] = +10;

        r.reconstruct();
        r.broaden();
        r.broaden();
        r.broaden();
        r.broaden();
        r = apply_1d_realspace_push(*q1d, r, 2); r.sum_down();
        r = apply_1d_realspace_push(*q1d, r, 1); r.sum_down();
        r = apply_1d_realspace_push(*q1d, r, 0); r.sum_down();

        return r;
    }


    void iterate_trotter(World& world, Convolution1D<double_complex>* G,
                        cvecfuncT& camo, cvecfuncT& cbmo, double t, double time_step,
                        double thresh);


    // For given protocol, solve the DFT/HF/response equations
    void solve(World & world);

};


// Computes molecular energy as a function of the geometry
// This is cludgy ... need better factorization of functionality
// between calculation, main program and this ... or just merge it all.
class MolecularEnergy : public OptimizationTargetInterface {
    World& world;
    SCF& calc;
    mutable double coords_sum;     // sum of square of coords at last solved geometry

public:
    MolecularEnergy(World& world, SCF& calc)
        : world(world)
        , calc(calc)
        , coords_sum(-1.0)
    {}

    bool provides_gradient() const {return true;}

    double value(const Tensor<double>& x) {
        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;

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

        // The below is missing convergence test logic, etc.

        // Make the nuclear potential, initial orbitals, etc.
        for (unsigned int proto=0; proto<calc.param.protocol_data.size(); proto++) {
            calc.set_protocol<3>(world,calc.param.protocol_data[proto]);
            calc.make_nuclear_potential(world);
            calc.project_ao_basis(world);

            if (proto == 0) {
                if (calc.param.restart) {
                    calc.load_mos(world);
                }
                else {
                    calc.initial_guess(world);
                    //calc.param.restart = true;
                }
            }
            else {
                calc.project(world);
            }

            // 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);
            }
            //if(proto==0){
            // if(calc.param.absolvent)//||calc.param.svpe)
            //     calc.solve_gas_phase(world);// computes vacuo density and energy
            //     //}//  calc.save_mos(world); //debugging

            calc.solve(world);
            if (calc.param.save)
              calc.save_mos(world);
        }
        return calc.current_energy;
    }

    madness::Tensor<double> gradient(const Tensor<double>& x) {
        value(x); // Ensures DFT equations are solved at this geometry

        return calc.derivatives(world);
    }
};
}

#endif /* SCF_H_ */