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
13double myxfunc(const madness::coord_3d& r);
14double myyfunc(const madness::coord_3d& r);
15
16namespace 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
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 myyfunc(const madness::coord_3d &r)
Definition DF.cc:26
double myxfunc(const madness::coord_3d &r)
Definition DF.cc:21
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
Abstract base class interface required for functors used as input to Functions.
Definition function_interface.h:68
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 tensor is a multidimension array.
Definition tensor.h:317
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:182
#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
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
static const char * filename
Definition legendre.cc:96
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
Function< std::complex< Q >, NDIM > function_real2complex(const Function< Q, NDIM > &r)
Definition complexfun.h:241
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
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:1671
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
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