CalculationParameters.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$
 */
 
/// \file CalculationParameters
/// \brief solution parameters for SCF calculations
 
 
 
#ifndef MADNESS_CHEM_CALCULATIONPARAMETERS_H__INCLUDED
#define MADNESS_CHEM_CALCULATIONPARAMETERS_H__INCLUDED
 
#include<madness/chem/molecule.h>
#include<madness/chem/molecularbasis.h>
#include<madness/mra/QCCalculationParametersBase.h>
#include<madness/mra/commandlineparser.h>
 
 
namespace madness {
 
struct CalculationParameters : public QCCalculationParametersBase {
  static constexpr char const* tag = "dft";
 
    CalculationParameters(const CalculationParameters& other) = default;
 
    CalculationParameters(World& world, const commandlineparser& parser) : CalculationParameters() {
        read_input_and_commandline_options(world, parser, "dft");
        // NB: `--optimize` deliberately does NOT set anything here. It selects a
        // workflow task (WorkflowBuilders::add_optimize_workflow_drivers); the
        // `dft` group has no say in whether a geometry is optimized.
        std::string inputfile=parser.value("input");
        // std::string prefix=commandlineparser::remove_extension(commandlineparser::base_name(inputfile));
        // if (prefix!="input") set_derived_value("prefix",prefix);
    }
 
    /// ctor reading out the input file
    CalculationParameters() {
        initialize<std::string>("prefix","mad","prefixes your output/restart/json/plot/etc files");
        initialize<double>("charge",0.0,"total molecular charge");
        initialize<std::string> ("xc","hf","XC input line");
        initialize<std::string> ("hfexalg","multiworld_row","hf exchange algorithm; multiworld bounds memory, needs one subworld per rank",{"multiworld","multiworld_row","fetch_compute","smallmem","largemem"});
        initialize<long>  ("hfex_granularity",1,"hfexalg=multiworld: orbital batches per rank; >1 balances better");
        initialize<bool>  ("hfex_cost_aware",true,"hfexalg=multiworld: place tasks by measured cost, not by count");
        initialize<int>   ("hfex_accumulation",2,"hfexalg=multiworld: gather tile results 1=per subworld, 2=per node");
        initialize<std::string>("bsh_apply","auto","BSH apply backend; auto=macrotask when multinode or at tight thresh, else tile",{"auto","tile","macrotask","plain"});
        initialize<std::vector<std::string>>("memory",{"storefunction","nodereplicated","distributed"},"memory algorithm for storing functions (storing,cloud,target)");
        initialize<double>("smear",0.0,"smearing parameter");
        initialize<double>("econv",1.e-5,"energy convergence");
        initialize<double>("dconv",1.e-4,"density convergence");
        initialize<std::vector<std::string> >("convergence_criteria",{"bsh_residual","total_energy"},"possible values are: bsh_residual, total_energy, each_energy, density");
        initialize<int>   ("k",-1,"polynomial order");
        initialize<double>("l",20,"user coordinates box size");
        initialize<std::string>("deriv","abgv","derivative method",{"abgv","bspline","ble"});
        initialize<std::string>("dft_deriv","bspline","derivative method for gga potentials",{"abgv","bspline","ble"});
        initialize<bool>  ("xc_weak_gga",false,"weak form of the semilocal xc potential; nemo only");
        initialize<double>("maxrotn",0.25,"step restriction used in autoshift algorithm");
        initialize<int>   ("nvalpha",0,"number of alpha virtuals to compute");
        initialize<int>   ("nvbeta",0,"number of beta virtuals to compute");
        initialize<std::vector<double> >("aocc",std::vector<double>(),"explicit alpha occupations (0 or 1) by orbital index within the occupied span; a 0 is a hole");
        initialize<std::vector<double> >("bocc",std::vector<double>(),"explicit beta occupations (0 or 1) by orbital index within the occupied span; a 0 is a hole");
        initialize<int>   ("nopen",0,"number of unpaired electrons = nalpha-nbeta");
        initialize<int>   ("maxiter",25,"maximum number of iterations");
        initialize<int>   ("nio",1,"no. of io servers to use");
        initialize<bool>  ("spin_restricted",true,"true if spin restricted");
        initialize<int>   ("plotlo",0,"range of MOs to print (for both spins if polarized");
        initialize<int>   ("plothi",-1,"range of MOs to print (for both spins if polarized");
        initialize<bool>  ("plotdens",false,"If true print the density at convergence");
        initialize<bool>  ("plotcoul",false,"If true plot the total coulomb potential at convergence");
        initialize<bool>  ("plotcube",false,"If true also write Gaussian .cube files (for Avogadro/VMD) alongside .dx");
        initialize<std::string> ("localize","new","localization method",{"pm","boys","new","new_sys","cholesky","canon"});
        initialize<std::string> ("pointgroup","c1","use point (sub) group symmetry if not localized",{"c1","c2","ci","cs","c2v","c2h","d2","d2h"});
        initialize<std::string>("restart","auto","where the initial orbitals come from: auto picks "
                "between none/iterate/read_only by what is on disk",
                {"auto","none","iterate","read_only","ao","nwchem"});
        // Retired in favour of the `restart` modes above, but kept registered so a
        // deck that still sets them gets a migration error instead of the single
        // easily-missed warning that ignore_unknown_keys would produce. See the
        // guard at the end of set_derived_values().
        initialize<bool>  ("restartao",false,"RETIRED -- use restart=ao");
        initialize<bool>  ("no_compute",false,"RETIRED -- use restart=read_only");
        initialize<bool>  ("save",true,"if true save orbitals to disk");
        initialize<int>   ("maxsub",10,"size of iterative subspace ... set to 0 or 1 to disable");
        initialize<double> ("orbitalshift",0.0,"scf orbital shift: shift the occ orbitals to lower energies");
        initialize<int>    ("npt_plot",101,"no. of points to use in each dim for plots");
//      initialize<Tensor<double> > ("plot_cell",Tensor<double>(),"lo hi in each dimension for plotting (default is all space)");
        initialize<std::vector<double> > ("plot_cell",std::vector<double>(),"lo hi in each dimension for plotting (default is all space)");
        initialize<std::string> ("aobasis","6-31g","AO basis used for initial guess (6-31gss, 6-31g, 3-21g, sto-6g, sto-3g)");
        initialize<bool> ("derivatives",false,"if true calculate nuclear derivatives");
        initialize<bool> ("dipole",false,"if true calculate dipole moment");
        initialize<bool> ("conv_only_dens",false,"if true remove bsh_residual from convergence criteria (deprecated)");
        initialize<bool> ("psp_calc",false,"pseudopotential calculation for all atoms");
        initialize<std::string> ("pcm_data","none","do a PCM (solvent) calculation");
        initialize<std::string> ("dispersion","none","DFT-D3 dispersion correction",{"none","d3bj","d3zero"});
        initialize<std::string> ("dispersion_functional","none","functional whose D3 damping parameters to use");
        initialize<bool> ("dispersion_atm",false,"include the three-body Axilrod-Teller-Muto dispersion term");
        initialize<bool> ("pure_ae",true,"pure all electron calculation with no pseudo-atoms");
        initialize<int>  ("print_level",3,"0: no output; 1: final energy; 2: iterations; 3: timings; 10: debug");
        initialize<std::vector<std::string> > ("population",{"none"},"atomic charges (and spin populations) of the final "
                "orbitals, printed at the end: none, mulliken, lowdin, iao; several allowed, e.g. [iao,mulliken]");
        initialize<std::string> ("population_basis","6-31g","basis set for the mulliken and lowdin populations "
                "(default: the value of aobasis)");
        initialize<std::string> ("population_minbasis","sto-3g","minimal basis set for the iao populations");
        initialize<std::string>  ("molecular_structure","inputfile","where to read the molecule from: inputfile or name from the library");
 
        // Next list inferred parameters
        initialize<int> ("nalpha",-1,"number of alpha spin electrons");
        initialize<int> ("nbeta",-1,"number of beta  spin electrons");
        initialize<int> ("nmo_alpha",-1,"number of alpha spin molecular orbitals");
        initialize<int> ("nmo_beta",-1,"number of beta spin molecular orbitals");
        initialize<double> ("lo",1.e-10,"smallest length scale we need to resolve");
        initialize<std::vector<double> > ("protocol",{1.e-4,1.e-6},"calculation protocol");
 
        // Retired geometry-optimization knobs. Geometry optimization is no longer
        // something an SCF does to itself; it is a workflow task of its own
        // (`madqc --optimize --wf=<scf|nemo>`, qcapp::OptimizeDriver) and every
        // knob it has lives in the `optimization` group (ParameterManager.hpp).
        //
        // These stay registered, and error, rather than being deleted outright:
        // ignore_unknown_keys is true by default, so a deleted key would let an
        // old deck run as a plain single point after one warning nobody reads --
        // silently returning an unoptimized geometry. See the guard at the end of
        // set_derived_values().
        initialize<bool> ("gopt",false,"RETIRED -- use `madqc --optimize`");
        initialize<double> ("gtol",1.e-4,"RETIRED -- use optimization group `gtol`");
        initialize<bool> ("gtest",false,"RETIRED -- never had an effect");
        initialize<double> ("gval",1.e-5,"RETIRED -- use optimization group `value_precision`");
        initialize<double> ("gprec",1.e-4,"RETIRED -- use optimization group `gradient_precision`");
        initialize<int> ("gmaxiter",20,"RETIRED -- use optimization group `maxiter`");
        initialize<bool> ("ginitial_hessian",false,"RETIRED -- use optimization group `initial_hessian`");
        initialize<std::string> ("algopt","bfgs","RETIRED -- use optimization group `algopt`",{"bfgs","cg"});
        initialize<int> ("nv_factor",1,"RETIRED -- use `nv_extra`/`nv_step`");
        initialize<int> ("nv_extra",0,"extra virtuals converged first and dropped stepwise down to nvalpha");
        initialize<int> ("nv_step",0,"virtuals dropped per step-down stage (0: all extras at once)");
        initialize<int> ("nv_its",5,"maximum iterations per intermediate step-down stage");
        initialize<bool> ("freeze_occupied",false,"iterate only the virtuals in the mean field of restarted occupied orbitals, which stay fixed");
        initialize<int> ("vnucextra",2,"load balance parameter for nuclear pot");
        initialize<int> ("loadbalparts",2,"??");
 
          //Keyword to use nwchem output for initial guess
          initialize<std::string> ("nwfile","none","Base name of nwchem output files (.out and .movecs extensions) to read from");
 
    }
 
    std::string get_tag() const override {
        return tag;
    }
 
    public:
    using QCCalculationParametersBase::read_input_and_commandline_options;
 
    std::string prefix() const {return get<std::string>("prefix");}
 
    double econv() const {return get<double>("econv");}
    double dconv() const {return get<double>("dconv");}
 
    bool converge_density() const {
        std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
        return std::find(criteria.begin(),criteria.end(),"density")!=criteria.end();
    }
    bool converge_bsh_residual() const {
        std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
        return std::find(criteria.begin(),criteria.end(),"bsh_residual")!=criteria.end();
    }
    bool converge_total_energy() const {
        std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
        return std::find(criteria.begin(),criteria.end(),"total_energy")!=criteria.end();
    }
    bool converge_each_energy() const {
        std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
        return std::find(criteria.begin(),criteria.end(),"each_energy")!=criteria.end();
    }
 
    int nopen() const {return get<int>("nopen");}
    int nalpha() const {return get<int>("nalpha");}
    int nbeta() const {return get<int>("nbeta");}
 
    int nvalpha() const {return get<int>("nvalpha");}
    int nvbeta() const {return get<int>("nvbeta");}
    std::vector<double> aocc() const {return get<std::vector<double> >("aocc");}
    std::vector<double> bocc() const {return get<std::vector<double> >("bocc");}
    int nv_extra() const {return get<int>("nv_extra");}
    int nv_step() const {return get<int>("nv_step");}
    int nv_its() const {return get<int>("nv_its");}
    bool freeze_occupied() const {return get<bool>("freeze_occupied");}
 
    int nmo_alpha() const {return get<int>("nmo_alpha");}
    int nmo_beta() const {return get<int>("nmo_beta");}
 
    bool have_beta() const {return (nbeta()>0) and (not spin_restricted());}
 
    bool spin_restricted() const {return get<bool>("spin_restricted");}
 
    double lo() const {return get<double>("lo");}
    double L() const {return get<double>("l");}
    int k() const {return get<int>("k");}
 
    std::string localize_method() const {return get<std::string>("localize");}
    bool do_localize() const {return (localize_method()!="canon");}
    bool localize_pm() const {return (localize_method()=="pm");}
 
    std::string pointgroup() const {return get<std::string>("pointgroup");}
    bool do_symmetry() const {return (pointgroup()!="c1");}
    double charge() const {return get<double>("charge");}
    int print_level() const {return get<int>("print_level");}
    std::vector<std::string> population() const {return get<std::vector<std::string> >("population");}
    std::string population_basis() const {return get<std::string>("population_basis");}
    std::string population_minbasis() const {return get<std::string>("population_minbasis");}
 
    int maxiter() const {return get<int>("maxiter");}
    double orbitalshift() const {return get<double>("orbitalshift");}
 
    std::string deriv() const {return get<std::string>("deriv");}
    std::string dft_deriv() const {return get<std::string>("dft_deriv");}
    bool xc_weak_gga() const {return get<bool>("xc_weak_gga");}
    std::string pcm_data() const {return get<std::string>("pcm_data");}
    std::string dispersion() const {return get<std::string>("dispersion");}
    std::string dispersion_functional() const {return get<std::string>("dispersion_functional");}
    bool dispersion_atm() const {return get<bool>("dispersion_atm");}
    std::string xc() const {return get<std::string>("xc");}
    std::string hfexalg() const {return get<std::string>("hfexalg");}
    long hfex_granularity() const {return get<long>("hfex_granularity");}
    bool hfex_cost_aware() const {return get<bool>("hfex_cost_aware");}
    int hfex_accumulation() const {return get<int>("hfex_accumulation");}
    std::string bsh_apply() const {return get<std::string>("bsh_apply");}
 
    std::vector<std::string> memory() const {return get<std::vector<std::string>>("memory");}
 
    std::string aobasis() const {return get<std::string>("aobasis");}
 
    std::vector<double> protocol() const {return get<std::vector<double> >("protocol");}
    bool save() const {return get<bool>("save");}
    /// the `restart` keyval as a string; parse with restart_mode_from_string()
 
    /// Returns the spelling rather than the RestartMode enum because RestartPlan.h
    /// includes this header, so the dependency cannot run the other way.
    std::string restart() const {return get<std::string>("restart");}
    bool restart_cphf() const {return get<bool>("restart_cphf");}
 
    int maxsub() const {return get<int>("maxsub");}
    double maxrotn() const {return get<double>("maxrotn");}
 
    int vnucextra() const {return get<int>("vnucextra");}
    int loadbalparts() const {return get<int>("loadbalparts");}
 
 
    bool derivatives() const {return get<bool>("derivatives");}
    bool dipole() const {return get<bool>("dipole");}
 
    // NB: no gopt()/gtol()/gval()/gprec()/gmaxiter()/ginitial_hessian()/algopt()
    // accessors. Those keyvals are retired; see the initialize() block above and
    // the guard in set_derived_values(). The optimizer's knobs are in
    // OptimizationParameters (chem/ParameterManager.hpp).
 
     std::string nwfile() const {return get<std::string>("nwfile");}
 
    Tensor<double> plot_cell() const {
        std::vector<double> vcell=get<std::vector<double> >("plot_cell");
        if (vcell.size()==0) return Tensor<double>();
        Tensor<double> cell(3,2);
        cell(0,0)=vcell[0];
        cell(0,1)=vcell[1];
        cell(1,0)=vcell[2];
        cell(1,1)=vcell[3];
        cell(2,0)=vcell[4];
        cell(2,1)=vcell[5];
        return cell;
    }
 
 
    void set_derived_values(const Molecule& molecule) {
 
        for (size_t iatom = 0; iatom < molecule.natom(); iatom++) {
            if (molecule.get_pseudo_atom(iatom)){
                set_derived_value("pure_ae",false);
                continue;
            }
        }
        set_derived_value("aobasis",molecule.guess_file());;
        const int n_core = molecule.n_core_orb_all();
 
 
        std::vector<double> proto=get<std::vector<double> >("protocol");
    // No ... The accuracy of computation is INDEPENDENT of the convergence requirement
    // --- actually need more precision than convergence threshold in order to have
    // variational principle working and for robust convergence
        //proto.back()=get<double>("econv");
    set_derived_value("protocol",proto);
    // No ... the energy is variational!  Don't need more accuracy in dconv --- in fact the opposite is true.
    // set_derived_value("dconv",sqrt(get<double>("econv"))*0.1);
 
        double z = molecule.total_nuclear_charge();
        const double charge=get<double>("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);
        }
 
        set_derived_value("nalpha",(nelec + nopen())/2);
        set_derived_value("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());
        if (nalpha() != nbeta()) set_derived_value("spin_restricted",false);
 
        set_derived_value("nmo_alpha",nalpha() + nvalpha());
        set_derived_value("nmo_beta",nbeta() + nvbeta());
 
        if (nv_extra() < 0 or nv_step() < 0 or nv_its() < 1) error("nv_extra, nv_step >= 0 and nv_its >= 1 required");
        if (nv_extra() > 0 and nvalpha() == 0) error("nv_extra requires nvalpha > 0");
        if (freeze_occupied() and nvalpha() == 0 and nvbeta() == 0) error("freeze_occupied requires nvalpha or nvbeta > 0");
 
        // Explicit occupations address the occupied span. A hole there is only
        // index-stable with canonical orbitals: every iteration re-sorts them by
        // eigenvalue, so the empty index stays on the same orbital.
        bool hole = false;
        for (const auto& [name, n] : std::vector<std::pair<std::string, int>>{{"aocc", nalpha()}, {"bocc", nbeta()}}) {
            const auto occ = get<std::vector<double> >(name);
            if (int(occ.size()) > n) error((name + " has more entries than occupied orbitals").c_str(), occ.size());
            for (const double o : occ) {
                if (o != 0.0 and o != 1.0) error((name + " entries must be 0 or 1").c_str(), o);
                if (o == 0.0) hole = true;
            }
        }
        if (not bocc().empty() and spin_restricted()) error("bocc requires spin_restricted false");
        if (hole) {
            set_derived_value("localize", std::string("canon"));
            if (do_localize()) error("a hole in the explicit occupations requires localize canon");
        }
 
        // 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
        set_derived_value("l",molecule.bounding_cube() + 50.0);
 
        //set_derived_value("lo",molecule.smallest_length_scale());
 
        // set highest possible point group for symmetry
        if (do_localize()) set_derived_value("pointgroup",std::string("c1"));
        else set_derived_value("pointgroup",molecule.get_pointgroup());
 
        // above two lines will not override user input, so check input is sane
        if (do_localize() and do_symmetry()) {
            error("\n\nsymmetry and localization cannot be used at the same time\n"
                    "switch from local to canonical orbitals (keyword canon)\n\n");
        }
 
        // Geometry optimization left the `dft` group entirely. `gopt` in particular
        // must not be waved through: a deck that sets it expects an optimized
        // geometry back, and silently returning a single point at the input
        // geometry is a wrong answer, not a missing feature. They stay registered
        // and erroring rather than being deleted, because ignore_unknown_keys is
        // true by default and a deleted key would run the deck as a plain single
        // point after one warning nobody reads.
        for (const char* key : {"gopt","gtol","gtest","gval","gprec",
                                       "gmaxiter","ginitial_hessian","algopt"}) {
            if (is_user_defined(key))
                error(("\n\n`" + std::string(key) + "` has been retired: geometry optimization is now a task of "
                       "its own.\nUse `madqc --optimize --wf=<scf|nemo>` and the `optimization` "
                       "parameter group\n(see `madqc --print_parameters=optimization`).\n\n").c_str());
        }
 
        // Retired keyvals. A clean break needs to be loud: ignore_unknown_keys is
        // true by default, so deleting these outright would let an old deck run
        // with different semantics after one warning nobody reads. They stay
        // registered and erroring until the migration is old news.
        if (is_user_defined("restartao"))
            error("\n\n`restartao` has been retired: use `restart ao` instead\n\n");
        if (is_user_defined("no_compute"))
            error("\n\n`no_compute` has been retired: use `restart read_only` instead\n\n");
        if (is_user_defined("nv_factor"))
            error("\n\n`nv_factor` has been retired: the virtual step-down is now `nv_extra`, "
                  "`nv_step` and `nv_its`\n\n");
 
        // population analysis: known schemes, `none` alone, and all-electron only (the
        // projected basis functions describe all electrons of an atom)
        {
            const std::vector<std::string> pop = population();
            bool any = false, none = false;
            for (const auto& s : pop) {
                MADNESS_CHECK_THROW(s == "none" or s == "mulliken" or s == "lowdin" or s == "iao",
                                    "population: the schemes are none, mulliken, lowdin and iao");
                any = any or s != "none";
                none = none or s == "none";
            }
            MADNESS_CHECK_THROW(not (any and none), "population: none cannot be combined with a scheme");
            MADNESS_CHECK_THROW(not any or (get<bool>("pure_ae") and not get<bool>("psp_calc") and n_core == 0),
                                "population: needs an all-electron calculation (no pseudopotentials or core potentials)");
            set_derived_value("population_basis", aobasis());
        }
 
        // dispersion correction
        if (dispersion()!="none" and xc()!="hf") {
            set_derived_value("dispersion_functional",xc());
        }
 
 
    }
};
 
 
} // namespace madness
 
#endif /* MADNESS_CHEM_CALCULATIONPARAMETERS_H__INCLUDED */