MADNESS  0.10.1
InitParameters.h
Go to the documentation of this file.
1 
2 /// \file InitParameters
3 /// \brief Input parameters for a Dirac Fock calculation, read from a specified archive resulting from a nonrelativistic moldft calculation or restarted from a previous Dirac Fock calculation.
4 
5 
6 #ifndef MADNESS_APPS_DFGUESSPARAMS_H_INCLUDED
7 #define MADNESS_APPS_DFGUESSPARAMS_H_INCLUDED
8 
9 #include "fcwf.h"
10 #include <madness/chem/NWChem.h>
11 
13 double myxfunc(const madness::coord_3d& r);
14 double myyfunc(const madness::coord_3d& r);
15 
16 namespace madness{
17 
18 
20  // Ground state parameters that are read in from archive
21  std::string inFile; ///< Name of input archive to read in
22  double Init_total_energy; ///< Total energy of the nonrelativistic ground state
23  bool spinrestricted; ///< Indicates if input calc. was spin-restricted
25  unsigned int num_occupied; ///< Number of orbitals
26  Tensor<double> energies; ///< Energies of input orbitals
27  Tensor<double> occ; ///< Occupancy of input orbitals
28  double L; ///< Box size of input - Dirac Fock calcluation is in same box
29  int order; ///< Order of polynomial used in input
30  Molecule molecule; ///< The molecule used in input calculation
31  std::vector<Fcwf> orbitals; ///< The occupied orbitals
32 
33  // Default constructor
35 
36  // Initializes InitParameters using the contents of file \c filename
37  void read(World& world, const std::string& filename, bool restart, bool Krestricted){
38  // Save the filename
39  inFile = filename;
40 
41  //First check to see if we're starting the job from a saved DF calculation rather than a moldft calculation
42  if(restart){
43  if(world.rank()==0) print("\n Reading initial data from restarted DF calculation");
44  archive::ParallelInputArchive input(world, filename.c_str());
45  input & Init_total_energy;
46  input & spinrestricted;
47  input & closed_shell;
48  input & num_occupied;
49  input & energies;
50  input & L;
51  input & order;
52  input & molecule;
53 
54  //Code breaks if spinrestricted (state of archive) and Krestricted (requested by user) don't match
55  //This functionality could probably be added at some point, but it's a bit of work
56  MADNESS_CHECK(spinrestricted == Krestricted);
57 
58  // Set this so we can read in whats
59  // written in the archive
62 
63  //Now we just have to unpack the orbitals
64  for(unsigned int i=0; i < num_occupied; i++){
65  Fcwf reader(world);
66  for(int j=0; j < 4; j++){
67  input & reader[j];
68  }
69  orbitals.push_back(copy(reader));
70  }
71  }
72  else{ //If we're not reading in from DF, then we're reading in from moldft
73 
74  //some dummy variables for reading/computing
75  std::vector<int> dummy2;
76  Tensor<double> temp_energies;
77 
78  //read in what's in the archive. See SCF::save_mos for how these archives are stored
79  archive::ParallelInputArchive input(world, filename.c_str());
80  unsigned int version=0;
81  std::string xc;
82  std::string localize_method;
83 
84  input & version;
85  input & Init_total_energy; // double
86  input & spinrestricted; // bool
87  input & L; // double box size
88  input & order; // int wavelet order
89  input & molecule; // Molecule
90  input & xc;
91  input & localize_method;
92 
93  input & num_occupied; // int
94  input & temp_energies; // Tensor<double> orbital energies
95  input & occ; // Tensor<double> orbital occupations
96  input & dummy2; // std::vector<int> sets of orbitals(?)
97 
98  //For now assume spin-restricted means closed shell in moldft
100 
101  // Check that order is positive and less than 30
102  if (order < 1 or order > 30){
103  if(world.rank() == 0) print("\n ***PLEASE NOTE***\n Invalid wavelet order read from archive, setting to 8.\n This seems to happen when the default wavelet order is used in moldft.");
104  order = 8;
105  }
106 
107  // Set this so we can read in whats
108  // written in the archive
111 
112  //Now read in the orbitals and construct Fcwfs from them.
113  complex_derivative_3d Dx(world,0);
114  complex_derivative_3d Dy(world,1);
115  complex_derivative_3d Dz(world,2);
116  //double myc = 137.0359895; //speed of light in atomic units
117  std::complex<double> myi(0,1);
118  if(spinrestricted){
119  //If the calculation was spin-restricted in moldft, then we only have "spin-up" orbitals
120 
121  //Initialize some functions for reading in the orbitals
122  real_function_3d reader;
123  complex_function_3d complexreader;
124  Fcwf spinup(world);
125  Fcwf spindown(world); //used if !Krestricted
126  real_function_3d xfunc = real_factory_3d(world).f(myxfunc);
127  real_function_3d yfunc = real_factory_3d(world).f(myyfunc);
128 
129  //Handle Kramers-restricted and unrestricted cases differently
130  if(Krestricted){
131  //Loop over the occupied orbitals and convert
132  for(unsigned int i = 0; i < num_occupied; i++){
133  //read in orbital
134  input & reader;
135 
136  //change to a complex function
137  complexreader = function_real2complex(reader);
138 
139  //build up the corresponding fcwfs, using kinetic balance to define the small component
140  spinup[0] = complexreader;
141  spinup[1] = complex_factory_3d(world);
142  spinup[2] = (-myi) * Dz(complexreader);
143  spinup[2].scale(0.5);
144  spinup[3] = (-myi) * (Dx(complexreader) + myi * Dy(complexreader));
145  spinup[3].scale(0.5);
146  spinup.normalize();
147  orbitals.push_back(spinup);
148  }
149 
150  //Update energies
151  energies = temp_energies;
152  }
153  else{
154  //Loop over the occupied orbitals and convert
155  for(unsigned int i = 0; i < num_occupied; i++){
156  //read in orbital
157  input & reader;
158 
159  //change to a complex function
160  complexreader = function_real2complex(reader);
161 
162  //build up the corresponding fcwfs, using kinetic balance to define the small component
163  spinup[0] = complexreader;
164  spinup[1] = complex_factory_3d(world);
165  spinup[2] = (-myi) * Dz(complexreader);
166  spinup[2].scale(0.5);
167  spinup[3] = (-myi) * (Dx(complexreader) + myi * Dy(complexreader));
168  spinup[3].scale(0.5);
169  spinup.normalize();
170  spindown[0] = complex_factory_3d(world);
171  spindown[1] = complexreader;
172  spindown[2] = (-myi) * (Dx(complexreader) - myi * Dy(complexreader));
173  spindown[2].scale(0.5);
174  spindown[3] = (myi) * Dz(complexreader);
175  spindown[3].scale(0.5);
176  spindown.normalize();
177  orbitals.push_back(spinup);
178  orbitals.push_back(spindown);
179  }
180 
181  //Double length of energies tensor and fill in as needed.
183  for(unsigned int i = 0; i < num_occupied; i++){
184  energies[2*i] = temp_energies[i];
185  energies[2*i+1] = temp_energies[i];
186  }
187  num_occupied *= 2;
188  }
189  }
190  else{
191 
192  if(world.rank()==0) print("number of alpha read in from moldft is:" ,num_occupied);
193 
194  // Read in alpha ground state orbitals
195  real_function_3d reader;
196  complex_function_3d complexreader;
197  Fcwf fcwfreader(world);
198  for(unsigned int i = 0; i < num_occupied; i++){
199  input & reader;
200  complexreader = function_real2complex(reader);
201  fcwfreader[0] = complexreader;
202  fcwfreader[1] = complex_factory_3d(world);
203  fcwfreader[2] = (-myi) * Dz(complexreader);
204  fcwfreader[2].scale(0.5);
205  fcwfreader[3] = (-myi) * (Dx(complexreader) + myi * Dy(complexreader));
206  fcwfreader[3].scale(0.5);
207  fcwfreader.normalize();
208  orbitals.push_back(fcwfreader);
209  }
210 
211  if(!Krestricted){
212  // Read in beta quantities
213  unsigned int num_betas=0;
214  input & num_betas;
215 
216  Tensor<double> beta_energies;
217  input & beta_energies;
218 
219  Tensor<double> dummy3;
220  input & dummy3;
221 
222  std::vector<int> dummy4;
223  input & dummy4;
224 
225  //read in beta ground state orbitals
226  for(unsigned int i = 0; i < num_betas; i++){
227  input & reader;
228  complexreader = function_real2complex(reader);
229  fcwfreader[0] = complex_factory_3d(world);
230  fcwfreader[1] = complexreader;
231  fcwfreader[2] = (-myi) * (Dx(complexreader) - myi * Dy(complexreader));
232  fcwfreader[2].scale(0.5);
233  fcwfreader[3] = (myi) * Dz(complexreader);
234  fcwfreader[3].scale(0.5);
235  fcwfreader.normalize();
236  orbitals.push_back(fcwfreader);
237  }
238 
239  //Handle energy tensor and number of occupied orbitals
240  energies = Tensor<double>(num_occupied + num_betas);
241  for(unsigned int i = 0; i < num_occupied; i++){
242  energies(i) = temp_energies(i);
243  }
244  for(unsigned int i = 0; i < num_betas; i++){
245  energies(num_occupied + i) = beta_energies(i);
246  }
247  num_occupied += num_betas;
248 
249  }
250  else{
251  energies = temp_energies;
252  }
253 
254 
255 
256  }
257 
258  //reorder orbitals and energies in ascending order, if necessary.
259  double tempdouble;
260  Fcwf fcwfreader(world);
261  for(unsigned int i = 0; i < num_occupied; i++){
262  for(unsigned int j = i+1; j < num_occupied; j++){
263  if(energies(j) < energies(i)){
264  if(world.rank()==0) print("swapping orbitals", i, " and ", j);
265  tempdouble = energies(j);
266  energies(j) = energies(i);
267  energies(i) = tempdouble;
268  fcwfreader = orbitals[j];
269  orbitals[j] = orbitals[i];
270  orbitals[i] = fcwfreader;
271  }
272  }
273  }
274  }
275  }
276 
277  //This function no longer works
278  //TODO: Update this function before using it
279  void readnw(World& world, const std::string& filename, bool Krestricted){
280  //Called to read in initial parameters from an nwchem output file
281 
282  //For now just use default values for L and order
283  order = 6;
284  L = 50.0;
287 
288  //Need to set this to something...
289  Init_total_energy = 0.0;
290 
291  //Construct interface object from slymer namespace
292  slymer::NWChem_Interface nwchem(filename,std::cout);
293 
294  //For parallel runs, silencing all but 1 slymer instance
295  if(world.rank() != 0) {
296  std::ostream dev_null(nullptr);
297  nwchem.err = dev_null;
298  }
299 
300  //Read in basis set
302 
303  //Read in the molecular orbital coefficients, energies, and occupancies
305 
306  //Need to construct a molecule object by ourselves
307  molecule = Molecule();
308  unsigned int anum;
309  double x,y,z,q;
310  for(unsigned int i=0; i < nwchem.atoms.size(); i++){
311  anum = symbol_to_atomic_number(nwchem.atoms[i].symbol);
312  q = anum*1.0;
313  x = nwchem.atoms[i].position[0];
314  y = nwchem.atoms[i].position[1];
315  z = nwchem.atoms[i].position[2];
316  molecule.add_atom(x,y,z,q,anum);
317  }
318 
319  //Find out how many orbitals we're dealing with by looking at the occupancies
320  unsigned int numalpha(0), numbeta(0);
321 
322  bool have_beta(false);
323  for(unsigned int i = 0; i < nwchem.beta_occupancies.size(); i++){
324  if(nwchem.beta_occupancies[i] > 0.0) have_beta = true;
325  }
326 
327  if(have_beta){
328  //we're reading from an unrestricted calculation
329  //and for now we will assume this is an open shell calculation
330  for(unsigned int i = 0; i < nwchem.occupancies.size(); i++){
331  if(nwchem.occupancies[i] == 1.0) numalpha+=1;
332  }
333  for(unsigned int i = 0; i < nwchem.beta_occupancies.size(); i++){
334  if(nwchem.beta_occupancies[i] == 1.0) numbeta+=1;
335  }
336 
337  //Right now DF can only handle a single unpaired electron
338  MADNESS_CHECK(numalpha-1 == numbeta);
339  }
340  else{
341  for(unsigned int i = 0; i < nwchem.occupancies.size(); i++){
342  if(nwchem.occupancies[i] == 2.0) numalpha += 1;
343  }
344  numbeta = numalpha;
345  }
346  closed_shell = !have_beta;
347 
348  //correctly set the number of occupied orbitals for the DF calculation
349  if(Krestricted){
350  num_occupied = numalpha;
351  }
352  else{
353  num_occupied = numalpha+numbeta;
354  }
355 
356 
357  //Let's print everything so we have a visual check on what we're working with (for now)
358  if(world.rank()==0) print("\nalpha occupancies:\n",nwchem.occupancies);
359  if(world.rank()==0) print("\nbeta occupancies:\n",nwchem.beta_occupancies);
360  if(world.rank()==0) print("\nenergies:\n",nwchem.energies);
361  if(world.rank()==0) print("\nbeta energies:\n",nwchem.beta_energies);
362  if(world.rank()==0) print("num alpha",numalpha);
363  if(world.rank()==0) print("num beta",numbeta);
364 
365 
366  //Now that we know how many orbitals we have. initialize and fill energy tensor
368  if(Krestricted){
369  for(unsigned int i=0; i < numalpha; i++){
370  energies[i] = nwchem.energies[i];
371  }
372  }
373  else{
374  if(closed_shell){
375  for(unsigned int i=0; i < numalpha; i++){
376  energies[2*i] = nwchem.energies[i];
377  energies[2*i+1] = nwchem.energies[i];
378  }
379  }
380  else{
381  for(unsigned int i=0; i < numalpha-1; i++){
382  energies[2*i] = nwchem.energies[i];
383  energies[2*i+1] = nwchem.energies[i];
384  }
385  energies[2*(numalpha-1)] = nwchem.energies[numalpha-1];
386  }
387  }
388 
389  //Cast the 'basis set' into a Gaussian basis and iterate over it
391  int ii = 0;
392  for(auto basis : slymer::cast_basis<slymer::GaussianFunction>(nwchem.basis_set)) {
393  //Get the center of gaussian as its special point
394  std::vector<coord_3d> centers;
395  coord_3d r;
396  r[0] = basis.get().center[0]; r[1] = basis.get().center[1]; r[2] = basis.get().center[2];
397  centers.push_back(r);
398 
399  //Now make the function
400  temp1.push_back(FunctionFactory<double,3>(world).functor(std::shared_ptr<FunctionFunctorInterface<double,3>>(new slymer::Gaussian_Functor(basis.get(), centers))));
401  double norm2 = temp1[ii].norm2();
402  if(world.rank() == 0) print("function", ii, "has norm", norm2);
403  ii++;
404  }
405 
406  //Normalize aos
407  normalize(world, temp1);
408 
409  //Transform aos now to get alpha mos
410  vector_real_function_3d temp = transform(world, temp1, nwchem.MOs , true);
411 
412  //Convert and store alpha occupied MOs.
413  complex_function_3d complexreader(world);
414  Fcwf spinup(world);
415  Fcwf spindown(world);
416  complex_derivative_3d Dx(world,0);
417  complex_derivative_3d Dy(world,1);
418  complex_derivative_3d Dz(world,2);
419  std::complex<double> myi(0,1);
420  for(unsigned int i = 0; i < numalpha-1; i++){
421  complexreader = function_real2complex(temp[i]);
422  spinup[0] = complexreader;
423  spinup[1] = complex_factory_3d(world);
424  spinup[2] = (-myi) * Dz(complexreader);
425  spinup[2].scale(0.5);
426  spinup[3] = (-myi) * (Dx(complexreader) + myi * Dy(complexreader));
427  spinup[3].scale(0.5);
428  spinup.normalize();
429  orbitals.push_back(spinup);
430  if(!Krestricted){
431  spindown[0] = complex_factory_3d(world);
432  spindown[1] = complexreader;
433  spindown[2] = (-myi) * (Dx(complexreader) - myi * Dy(complexreader));
434  spindown[2].scale(0.5);
435  spindown[3] = (myi) * Dz(complexreader);
436  spindown[3].scale(0.5);
437  spindown.normalize();
438  orbitals.push_back(spindown);
439  }
440  }
441  complexreader = function_real2complex(temp[numalpha-1]);
442  spinup[0] = complexreader;
443  spinup[1] = complex_factory_3d(world);
444  spinup[2] = (-myi) * Dz(complexreader);
445  spinup[2].scale(0.5);
446  spinup[3] = (-myi) * (Dx(complexreader) + myi * Dy(complexreader));
447  spinup[3].scale(0.5);
448  spinup.normalize();
449  orbitals.push_back(spinup);
450  if(closed_shell and not Krestricted){
451  spindown[0] = complex_factory_3d(world);
452  spindown[1] = complexreader;
453  spindown[2] = (-myi) * (Dx(complexreader) - myi * Dy(complexreader));
454  spindown[2].scale(0.5);
455  spindown[3] = (myi) * Dz(complexreader);
456  spindown[3].scale(0.5);
457  spindown.normalize();
458  orbitals.push_back(spindown);
459  }
460 
461  //Assure that the numbers line up
463 
464  }
465 
466  // Prints all information
467  void print_params() const
468  {
469  madness::print("\n Input Parameters");
470  madness::print(" -------------------------");
471  madness::print(" Input Archive:", inFile);
472  madness::print(" Spin Restricted:", spinrestricted);
473  madness::print(" No. of occ. orbitals:", num_occupied);
474  madness::print(" L:", L);
475  madness::print(" Wavelet Order:", order);
476  madness::print(" Initial Total Energy:", Init_total_energy);
477  madness::print(" Orbital Energies:", energies);
478  }
479  };
480 }
481 #endif
482 
483 //kthxbye
double q(double t)
Definition: DKops.h:18
double myyfunc(const madness::coord_3d &r)
Definition: DF.cc:26
double myxfunc(const madness::coord_3d &r)
Definition: DF.cc:21
Function< std::complex< double >, 3 > function_real2complex(const Function< double, 3 > &r)
Definition: DF.cc:97
Definition: fcwf.h:12
void scale(std::complex< double > a)
Definition: fcwf.cc:144
void normalize()
Definition: fcwf.cc:203
long size() const
Returns the number of elements in the tensor.
Definition: basetensor.h:138
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition: derivative.h:266
static void set_k(int value)
Sets the default wavelet order.
Definition: funcdefaults.h:273
static void set_cubic_cell(double lo, double hi)
Sets the user cell to be cubic with each dimension having range [lo,hi].
Definition: funcdefaults.h:461
FunctionFactory implements the named-parameter idiom for Function.
Definition: function_factory.h:86
FunctionFactory & functor(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > &f)
Definition: function_factory.h:141
A multiresolution adaptive numerical function.
Definition: mra.h:122
Definition: molecule.h:124
void add_atom(double x, double y, double z, double q, int atn)
Definition: molecule.cc:346
A parallel world class.
Definition: world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
An archive for storing local or parallel data, wrapping a BinaryFstreamInputArchive.
Definition: parallel_archive.h:366
const madness::Tensor< double > & MOs
Publically accessible alpha MO expansions coefficients. Column is the MO, row is the basis function.
Definition: ESInterface.h:76
const madness::Tensor< double > & beta_occupancies
Publically accessible list of beta MO occupancies (in eV).
Definition: ESInterface.h:80
const madness::Tensor< double > & beta_energies
Publically accessible list of beta MO energies (in eV).
Definition: ESInterface.h:78
const madness::Tensor< double > & energies
Publically accessible list of alpha MO energies.
Definition: ESInterface.h:75
const Atoms & atoms
Publically accessible list of atoms.
Definition: ESInterface.h:73
const BasisSet & basis_set
Publicly accessible basis set.
Definition: ESInterface.h:72
std::reference_wrapper< std::ostream > err
Output stream for messages.
Definition: ESInterface.h:70
const madness::Tensor< double > & occupancies
Publically accessible list of alpha MO occupancies (in eV).
Definition: ESInterface.h:77
Definition: gaussian.h:376
Class for interfacing with NWChem (tested on version 6.6).
Definition: NWChem.h:22
virtual void read(Properties::Properties props) override
Read the specified properties and store them in the member variables.
Definition: NWChem.cc:26
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition: madness_exception.h:190
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
const char * version()
Get the MADNESS version number.
Definition: info.cc:17
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
static const char * filename
Definition: legendre.cc:96
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:851
std::vector< real_function_3d > vector_real_function_3d
Definition: functypedefs.h:79
unsigned int symbol_to_atomic_number(const std::string &symbol)
Definition: atomutil.cc:173
Function< std::complex< Q >, NDIM > function_real2complex(const Function< Q, NDIM > &r)
Definition: complexfun.h:241
FunctionFactory< double, 3 > real_factory_3d
Definition: functypedefs.h:93
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition: print.h:225
FunctionFactory< double_complex, 3 > complex_factory_3d
Definition: functypedefs.h:100
void normalize(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Normalizes a vector of functions — v[i] = v[i].scale(1.0/v[i].norm2())
Definition: vmra.h:1614
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition: vmra.h:689
constexpr Properties Occupancies
MO occupancies.
Definition: ESInterface.h:52
constexpr Properties MOs
The MO vector coefficients.
Definition: ESInterface.h:51
constexpr Properties Energies
The MO energies.
Definition: ESInterface.h:50
constexpr Properties Basis
The basis set.
Definition: ESInterface.h:48
XCfunctional xc
Definition: newsolver_lda.cc:53
Definition: InitParameters.h:19
double Init_total_energy
Total energy of the nonrelativistic ground state.
Definition: InitParameters.h:22
bool closed_shell
Definition: InitParameters.h:24
int order
Order of polynomial used in input.
Definition: InitParameters.h:29
std::string inFile
Name of input archive to read in.
Definition: InitParameters.h:21
Molecule molecule
The molecule used in input calculation.
Definition: InitParameters.h:30
unsigned int num_occupied
Number of orbitals.
Definition: InitParameters.h:25
double L
Box size of input - Dirac Fock calcluation is in same box.
Definition: InitParameters.h:28
void readnw(World &world, const std::string &filename, bool Krestricted)
Definition: InitParameters.h:279
Tensor< double > occ
Occupancy of input orbitals.
Definition: InitParameters.h:27
void read(World &world, const std::string &filename, bool restart, bool Krestricted)
Definition: InitParameters.h:37
Tensor< double > energies
Energies of input orbitals.
Definition: InitParameters.h:26
void print_params() const
Definition: InitParameters.h:467
std::vector< Fcwf > orbitals
The occupied orbitals.
Definition: InitParameters.h:31
InitParameters()
Definition: InitParameters.h:34
bool spinrestricted
Indicates if input calc. was spin-restricted.
Definition: InitParameters.h:23