MADNESS  0.10.1
CalculationParameters.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 
31 
32  $Id$
33  */
34 
35 /// \file CalculationParameters
36 /// \brief solution parameters for SCF calculations
37 
38 
39 
40 #ifndef MADNESS_CHEM_CALCULATIONPARAMETERS_H__INCLUDED
41 #define MADNESS_CHEM_CALCULATIONPARAMETERS_H__INCLUDED
42 
47 
48 
49 namespace madness {
50 
52 
54 
56  read_input_and_commandline_options(world, parser, "dft");
57  // convenience option -- needs to be moved to the MolecularOptimizer class
58  if (parser.key_exists("optimize")) set_user_defined_value("gopt",true);
59  }
60 
61  /// ctor reading out the input file
63 
64  initialize<std::string>("prefix","mad","prefixes your output/restart/json/plot/etc files");
65  initialize<double>("charge",0.0,"total molecular charge");
66  initialize<std::string> ("xc","hf","XC input line");
67  initialize<std::string> ("hfexalg","multiworld","hf exchange algorithm: choose from multiworld (default), smallmem, largemem");
68  initialize<double>("smear",0.0,"smearing parameter");
69  initialize<double>("econv",1.e-5,"energy convergence");
70  initialize<double>("dconv",1.e-4,"density convergence");
71  initialize<std::vector<std::string> >("convergence_criteria",{"bsh_residual","total_energy"},"possible values are: bsh_residual, total_energy, each_energy, density");
72  initialize<int> ("k",-1,"polynomial order");
73  initialize<double>("l",20,"user coordinates box size");
74  initialize<std::string>("deriv","abgv","derivative method",{"abgv","bspline","ble"});
75  initialize<std::string>("dft_deriv","abgv","derivative method for gga potentials",{"abgv","bspline","ble"});
76  initialize<double>("maxrotn",0.25,"step restriction used in autoshift algorithm");
77  initialize<int> ("nvalpha",0,"number of alpha virtuals to compute");
78  initialize<int> ("nvbeta",0,"number of beta virtuals to compute");
79  initialize<int> ("nopen",0,"number of unpaired electrons = nalpha-nbeta");
80  initialize<int> ("maxiter",25,"maximum number of iterations");
81  initialize<int> ("nio",1,"no. of io servers to use");
82  initialize<bool> ("spin_restricted",true,"true if spin restricted");
83  initialize<int> ("plotlo",0,"range of MOs to print (for both spins if polarized");
84  initialize<int> ("plothi",-1,"range of MOs to print (for both spins if polarized");
85  initialize<bool> ("plotdens",false,"If true print the density at convergence");
86  initialize<bool> ("plotcoul",false,"If true plot the total coulomb potential at convergence");
87  initialize<std::string> ("localize","new","localization method",{"pm","boys","new","canon"});
88  initialize<std::string> ("pointgroup","c1","use point (sub) group symmetry if not localized",{"c1","c2","ci","cs","c2v","c2h","d2","d2h"});
89  initialize<bool> ("restart",false,"if true restart from orbitals on disk");
90  initialize<bool> ("restartao",false,"if true restart from orbitals projected into AO basis (STO3G) on disk");
91  initialize<bool> ("no_compute",false,"if true use orbitals on disk, set value to computed");
92  initialize<bool> ("save",true,"if true save orbitals to disk");
93  initialize<int> ("maxsub",10,"size of iterative subspace ... set to 0 or 1 to disable");
94  initialize<double> ("orbitalshift",0.0,"scf orbital shift: shift the occ orbitals to lower energies");
95  initialize<int> ("npt_plot",101,"no. of points to use in each dim for plots");
96 // initialize<Tensor<double> > ("plot_cell",Tensor<double>(),"lo hi in each dimension for plotting (default is all space)");
97  initialize<std::vector<double> > ("plot_cell",std::vector<double>(),"lo hi in each dimension for plotting (default is all space)");
98  initialize<std::string> ("aobasis","6-31g","AO basis used for initial guess (6-31gss, 6-31g, 3-21g, sto-6g, sto-3g)");
99  initialize<bool> ("derivatives",false,"if true calculate nuclear derivatives");
100  initialize<bool> ("dipole",false,"if true calculate dipole moment");
101  initialize<bool> ("conv_only_dens",false,"if true remove bsh_residual from convergence criteria (deprecated)");
102  initialize<bool> ("psp_calc",false,"pseudopotential calculation for all atoms");
103  initialize<std::string> ("pcm_data","none","do a PCM (solvent) calculation");
104  initialize<std::string> ("ac_data","none","do a calculation with asymptotic correction (see ACParameters class in chem/AC.h for details)");
105  initialize<bool> ("pure_ae",true,"pure all electron calculation with no pseudo-atoms");
106  initialize<int> ("print_level",3,"0: no output; 1: final energy; 2: iterations; 3: timings; 10: debug");
107  initialize<std::string> ("molecular_structure","inputfile","where to read the molecule from: inputfile or name from the library");
108 
109  // Next list inferred parameters
110  initialize<int> ("nalpha",-1,"number of alpha spin electrons");
111  initialize<int> ("nbeta",-1,"number of beta spin electrons");
112  initialize<int> ("nmo_alpha",-1,"number of alpha spin molecular orbitals");
113  initialize<int> ("nmo_beta",-1,"number of beta spin molecular orbitals");
114  initialize<double> ("lo",1.e-10,"smallest length scale we need to resolve");
115  initialize<std::vector<double> > ("protocol",{1.e-4,1.e-6},"calculation protocol");
116 
117  // geometry optimization parameters
118  // @TODO: need to be moved to molecular optimizer class
119  initialize<bool> ("gopt",false,"geometry optimizer");
120  initialize<double> ("gtol",1.e-4,"geometry tolerance");
121  initialize<bool> ("gtest",false,"geometry tolerance");
122  initialize<double> ("gval",1.e-5,"value precision");
123  initialize<double> ("gprec",1.e-4,"gradient precision");
124  initialize<int> ("gmaxiter",20,"optimization maxiter");
125  initialize<bool> ("ginitial_hessian",false,"compute inital hessian for optimization");
126  initialize<std::string> ("algopt","bfgs","algorithm used for optimization",{"bfgs","cg"});
127  initialize<int> ("nv_factor",1,"factor to multiply number of virtual orbitals with when automatically decreasing nvirt");
128  initialize<int> ("vnucextra",2,"load balance parameter for nuclear pot");
129  initialize<int> ("loadbalparts",2,"??");
130 
131  //Keyword to use nwchem output for initial guess
132  initialize<std::string> ("nwfile","none","Base name of nwchem output files (.out and .movecs extensions) to read from");
133 
134  }
135 
136  public:
138 
139  std::string prefix() const {return get<std::string>("prefix");}
140 
141  double econv() const {return get<double>("econv");}
142  double dconv() const {return get<double>("dconv");}
143 
144  bool converge_density() const {
145  std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
146  return std::find(criteria.begin(),criteria.end(),"density")!=criteria.end();
147  }
148  bool converge_bsh_residual() const {
149  std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
150  return std::find(criteria.begin(),criteria.end(),"bsh_residual")!=criteria.end();
151  }
152  bool converge_total_energy() const {
153  std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
154  return std::find(criteria.begin(),criteria.end(),"total_energy")!=criteria.end();
155  }
156  bool converge_each_energy() const {
157  std::vector<std::string> criteria=get<std::vector<std::string> >("convergence_criteria");
158  return std::find(criteria.begin(),criteria.end(),"each_energy")!=criteria.end();
159  }
160 
161  int nopen() const {return get<int>("nopen");}
162  int nalpha() const {return get<int>("nalpha");}
163  int nbeta() const {return get<int>("nbeta");}
164 
165  int nvalpha() const {return get<int>("nvalpha");}
166  int nvbeta() const {return get<int>("nvbeta");}
167  int nv_factor() const {return get<int>("nv_factor");}
168 
169  int nmo_alpha() const {return get<int>("nmo_alpha");}
170  int nmo_beta() const {return get<int>("nmo_beta");}
171 
172  bool have_beta() const {return (nbeta()>0) and (not spin_restricted());}
173 
174  bool spin_restricted() const {return get<bool>("spin_restricted");}
175  bool no_compute() const {return get<bool>("no_compute");}
176 
177  double lo() const {return get<double>("lo");}
178  double L() const {return get<double>("l");}
179  int k() const {return get<int>("k");}
180 
181  std::string localize_method() const {return get<std::string>("localize");}
182  bool do_localize() const {return (localize_method()!="canon");}
183  bool localize_pm() const {return (localize_method()=="pm");}
184 
185  std::string pointgroup() const {return get<std::string>("pointgroup");}
186  bool do_symmetry() const {return (pointgroup()!="c1");}
187  double charge() const {return get<double>("charge");}
188  int print_level() const {return get<int>("print_level");}
189 
190  int maxiter() const {return get<int>("maxiter");}
191  double orbitalshift() const {return get<double>("orbitalshift");}
192 
193  std::string deriv() const {return get<std::string>("deriv");}
194  std::string dft_deriv() const {return get<std::string>("dft_deriv");}
195  std::string pcm_data() const {return get<std::string>("pcm_data");}
196  std::string ac_data() const {return get<std::string>("ac_data");}
197  std::string xc() const {return get<std::string>("xc");}
198  std::string hfexalg() const {return get<std::string>("hfexalg");}
199 
200  std::string aobasis() const {return get<std::string>("aobasis");}
201 
202  std::vector<double> protocol() const {return get<std::vector<double> >("protocol");}
203  bool save() const {return get<bool>("save");}
204  bool restart() const {return get<bool>("restart");}
205  bool restartao() const {return get<bool>("restartao");}
206  bool restart_cphf() const {return get<bool>("restart_cphf");}
207 
208  int maxsub() const {return get<int>("maxsub");}
209  double maxrotn() const {return get<double>("maxrotn");}
210 
211  int vnucextra() const {return get<int>("vnucextra");}
212  int loadbalparts() const {return get<int>("loadbalparts");}
213 
214 
215  bool derivatives() const {return get<bool>("derivatives");}
216  bool dipole() const {return get<bool>("dipole");}
217 
218  bool gopt() const {return get<bool>("gopt");}
219  std::string algopt() const {return get<std::string>("algopt");}
220  int gmaxiter() const {return get<int>("gmaxiter");}
221  double gtol() const {return get<double>("gtol");}
222  double gval() const {return get<double>("gval");}
223  double gprec() const {return get<double>("gprec");}
224  bool ginitial_hessian() const {return get<bool>("ginitial_hessian");}
225 
226  std::string nwfile() const {return get<std::string>("nwfile");}
227 
229  std::vector<double> vcell=get<std::vector<double> >("plot_cell");
230  if (vcell.size()==0) return Tensor<double>();
231  Tensor<double> cell(3,2);
232  cell(0,0)=vcell[0];
233  cell(0,1)=vcell[1];
234  cell(1,0)=vcell[2];
235  cell(1,1)=vcell[3];
236  cell(2,0)=vcell[4];
237  cell(2,1)=vcell[5];
238  return cell;
239  }
240 
241 
243  std::string inputfile=parser.value("input");
245  if (prefix!="input") set_derived_value("prefix",prefix);
246 
247  for (size_t iatom = 0; iatom < molecule.natom(); iatom++) {
248  if (molecule.get_pseudo_atom(iatom)){
249  set_derived_value("pure_ae",false);
250  continue;
251  }
252  }
253  set_derived_value("aobasis",molecule.guess_file());;
254  const int n_core = molecule.n_core_orb_all();
255 
256 
257  std::vector<double> proto=get<std::vector<double> >("protocol");
258  // No ... The accuracy of computation is INDEPENDENT of the convergence requirement
259  // --- actually need more precision than convergence threshold in order to have
260  // variational principle working and for robust convergence
261  //proto.back()=get<double>("econv");
262  set_derived_value("protocol",proto);
263  // No ... the energy is variational! Don't need more accuracy in dconv --- in fact the opposite is true.
264  // set_derived_value("dconv",sqrt(get<double>("econv"))*0.1);
265 
266  double z = molecule.total_nuclear_charge();
267  const double charge=get<double>("charge");
268  int nelec = int(z - charge - n_core*2);
269  if (fabs(nelec+charge+n_core*2-z) > 1e-6) {
270  error("non-integer number of electrons?", nelec+charge+n_core*2-z);
271  }
272 
273  set_derived_value("nalpha",(nelec + nopen())/2);
274  set_derived_value("nbeta",(nelec - nopen())/2);
275 
276  if (nalpha() < 0) error("negative number of alpha electrons?", nalpha());
277  if (nbeta() < 0) error("negative number of beta electrons?", nbeta());
278  if ((nalpha()+nbeta()) != nelec) error("nalpha+nbeta != nelec", nalpha()+nbeta());
279  if (nalpha() != nbeta()) set_derived_value("spin_restricted",false);
280 
281  set_derived_value("nmo_alpha",nalpha() + nvalpha());
282  set_derived_value("nmo_beta",nbeta() + nvbeta());
283 
284  // Ensure we have enough basis functions to guess the requested
285  // number of states ... a minimal basis for a closed-shell atom
286  // might not have any functions for virtuals.
287  int nbf = aobasis.nbf(molecule);
288  if ((nmo_alpha()>nbf) or (nmo_beta()>nbf)) error("too few basis functions?", nbf);
289 // nmo_alpha = std::min(nbf,nmo_alpha);
290 // nmo_beta = std::min(nbf,nmo_beta);
291 // if (nalpha>nbf || nbeta>nbf) error("too few basis functions?", nbf);
292 // nvalpha = nmo_alpha - nalpha;
293 // nvbeta = nmo_beta - nbeta;
294 
295  // Unless overridden by the user use a cell big enough to
296  // have exp(-sqrt(2*I)*r) decay to 1e-6 with I=1ev=0.037Eh
297  // --> need 50 a.u. either side of the molecule
298  set_derived_value("l",molecule.bounding_cube() + 50.0);
299 
300  set_derived_value("lo",molecule.smallest_length_scale());
301 
302  // set highest possible point group for symmetry
303  if (do_localize()) set_derived_value("pointgroup",std::string("c1"));
304  else set_derived_value("pointgroup",molecule.get_pointgroup());
305 
306  // above two lines will not override user input, so check input is sane
307  if (do_localize() and do_symmetry()) {
308  error("\n\nsymmetry and localization cannot be used at the same time\n"
309  "switch from local to canonical orbitals (keyword canon)\n\n");
310  }
311 
312  //NWChem interface doesn't support geometry optimization
313  if (get<bool>("gopt") && nwfile() != "none") error("NWchem initialization only supports single point energy calculations.");
314 
315  //NWChem only supports Boys localization (or canonical)
316  if (nwfile() != "none") {
317  set_derived_value("localize",std::string("boys"));
318  //Error if user requested something other than Boys
319  if(localize_method() != "boys" and localize_method() != "canon") error("NWchem initialization only supports Boys localization");
320  }
321  }
322 };
323 
324 
325 } // namespace madness
326 
327 #endif /* MADNESS_CHEM_CALCULATIONPARAMETERS_H__INCLUDED */
Contracted Gaussian basis.
Definition: madness/chem/molecularbasis.h:465
Definition: molecule.h:124
class for holding the parameters for calculation
Definition: QCCalculationParametersBase.h:290
virtual void read_input_and_commandline_options(World &world, const commandlineparser &parser, const std::string tag)
Definition: QCCalculationParametersBase.h:325
void set_user_defined_value(const std::string &key, const T &value)
Definition: QCCalculationParametersBase.h:533
void set_derived_value(const std::string &key, const T &value)
Definition: QCCalculationParametersBase.h:403
A parallel world class.
Definition: world.h:132
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
void error(const char *msg)
Definition: world.cc:139
Definition: test_QCCalculationParametersBase.cc:57
Definition: CalculationParameters.h:51
Tensor< double > plot_cell() const
Definition: CalculationParameters.h:228
double gval() const
Definition: CalculationParameters.h:222
bool converge_each_energy() const
Definition: CalculationParameters.h:156
double charge() const
Definition: CalculationParameters.h:187
void set_derived_values(const Molecule &molecule, const AtomicBasisSet &aobasis, const commandlineparser &parser)
Definition: CalculationParameters.h:242
CalculationParameters(const CalculationParameters &other)=default
bool converge_density() const
Definition: CalculationParameters.h:144
virtual void read_input_and_commandline_options(World &world, const commandlineparser &parser, const std::string tag)
Definition: QCCalculationParametersBase.h:325
int nv_factor() const
Definition: CalculationParameters.h:167
bool do_symmetry() const
Definition: CalculationParameters.h:186
bool localize_pm() const
Definition: CalculationParameters.h:183
bool gopt() const
Definition: CalculationParameters.h:218
int k() const
Definition: CalculationParameters.h:179
bool converge_bsh_residual() const
Definition: CalculationParameters.h:148
CalculationParameters()
ctor reading out the input file
Definition: CalculationParameters.h:62
bool converge_total_energy() const
Definition: CalculationParameters.h:152
std::string nwfile() const
Definition: CalculationParameters.h:226
std::string prefix() const
Definition: CalculationParameters.h:139
std::string ac_data() const
Definition: CalculationParameters.h:196
int nmo_alpha() const
Definition: CalculationParameters.h:169
int nopen() const
Definition: CalculationParameters.h:161
bool dipole() const
Definition: CalculationParameters.h:216
bool save() const
Definition: CalculationParameters.h:203
double dconv() const
Definition: CalculationParameters.h:142
std::string algopt() const
Definition: CalculationParameters.h:219
double econv() const
Definition: CalculationParameters.h:141
int nvalpha() const
Definition: CalculationParameters.h:165
int vnucextra() const
Definition: CalculationParameters.h:211
int print_level() const
Definition: CalculationParameters.h:188
std::string deriv() const
Definition: CalculationParameters.h:193
bool do_localize() const
Definition: CalculationParameters.h:182
std::string hfexalg() const
Definition: CalculationParameters.h:198
double L() const
Definition: CalculationParameters.h:178
int gmaxiter() const
Definition: CalculationParameters.h:220
std::string localize_method() const
Definition: CalculationParameters.h:181
std::vector< double > protocol() const
Definition: CalculationParameters.h:202
bool ginitial_hessian() const
Definition: CalculationParameters.h:224
int loadbalparts() const
Definition: CalculationParameters.h:212
std::string aobasis() const
Definition: CalculationParameters.h:200
int nalpha() const
Definition: CalculationParameters.h:162
std::string dft_deriv() const
Definition: CalculationParameters.h:194
int maxsub() const
Definition: CalculationParameters.h:208
int nvbeta() const
Definition: CalculationParameters.h:166
double gprec() const
Definition: CalculationParameters.h:223
bool derivatives() const
Definition: CalculationParameters.h:215
int nbeta() const
Definition: CalculationParameters.h:163
int nmo_beta() const
Definition: CalculationParameters.h:170
double gtol() const
Definition: CalculationParameters.h:221
bool have_beta() const
Definition: CalculationParameters.h:172
bool no_compute() const
Definition: CalculationParameters.h:175
bool spin_restricted() const
Definition: CalculationParameters.h:174
int maxiter() const
Definition: CalculationParameters.h:190
double orbitalshift() const
Definition: CalculationParameters.h:191
double lo() const
Definition: CalculationParameters.h:177
std::string pcm_data() const
Definition: CalculationParameters.h:195
bool restart() const
Definition: CalculationParameters.h:204
bool restart_cphf() const
Definition: CalculationParameters.h:206
bool restartao() const
Definition: CalculationParameters.h:205
std::string pointgroup() const
Definition: CalculationParameters.h:185
double maxrotn() const
Definition: CalculationParameters.h:209
std::string xc() const
Definition: CalculationParameters.h:197
CalculationParameters(World &world, const commandlineparser &parser)
Definition: CalculationParameters.h:55
very simple command line parser
Definition: commandlineparser.h:15
std::string value(const std::string key) const
Definition: commandlineparser.h:59
bool key_exists(std::string key) const
Definition: commandlineparser.h:55
static std::string base_name(std::string const &path, std::string const &delims="/")
Definition: commandlineparser.h:129
static std::string remove_extension(std::string const &filename)
Definition: commandlineparser.h:134
void e()
Definition: test_sig.cc:75
static Molecule molecule
Definition: testperiodicdft.cc:38