1 /*
2 * This source code is part of
3 *
4 * HelFEM
5 * -
6 * Finite element methods for electronic structure calculations on small systems
7 *
8 * Written by Susi Lehtola, 2018-
9 * Copyright (c) 2018- Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16 #include "../general/cmdline.h"
17 #include "../general/checkpoint.h"
18 #include "../general/constants.h"
19 #include "../general/timer.h"
20 #include "utils.h"
21 #include "basis.h"
22 #include <cfloat>
23 #include <climits>
24
25 using namespace helfem;
26
main(int argc,char ** argv)27 int main(int argc, char **argv) {
28 cmdline::parser parser;
29
30 // full option name, no short option, description, argument required
31 parser.add<std::string>("load", 0, "load guess from checkpoint", false, "");
32 parser.add<double>("mumax", 0, "mu max (negative for same as input)", false, -1.0);
33 parser.add<double>("phi", 0, "value of phi angle", false, 0.0);
34 parser.add<int>("Nmu", 0, "number of points in mu (negative for same as input)", false, -1);
35 parser.add<int>("Nnu", 0, "number of points in nu (negative for same as input)", false, -1);
36 parser.add<bool>("offset", 0, "use nu and mu offset?", false, false);
37 parser.add<std::string>("savedens", 0, "save density to file", false, "density.hdf5");
38 parser.parse_check(argc, argv);
39
40 // Get parameters
41 std::string load(parser.get<std::string>("load"));
42 double mumax(parser.get<double>("mumax"));
43 double phi(parser.get<double>("phi"));
44 std::string savedens(parser.get<std::string>("savedens"));
45
46 // Load checkpoint
47 Checkpoint loadchk(load,false);
48 // Basis set
49 diatomic::basis::TwoDBasis basis;
50 loadchk.read(basis);
51
52 // Sanity check
53 if(mumax<0 || mumax>basis.get_mumax())
54 mumax=basis.get_mumax();
55
56 int Nmu(parser.get<int>("Nmu"));
57 if(Nmu<=0)
58 Nmu=2*basis.Nrad();
59
60 int Nnu(parser.get<int>("Nnu"));
61 if(Nnu<=0)
62 Nnu=2*arma::max(basis.get_lval());
63
64 printf("Grid spanning mu = 0 .. %e, Nmu = %i, Nnu = %i\n",mumax,Nmu,Nnu);
65
66 // Density matrix
67 arma::mat P, Pa, Pb;
68 loadchk.read("Pa",Pa);
69 loadchk.read("Pb",Pb);
70 loadchk.read("P",P);
71
72 // mu array
73 arma::vec mu, nu;
74 bool offset(parser.get<bool>("offset"));
75 if(offset) {
76 double dmu=mumax/Nmu;
77 double dnu=M_PI/Nnu;
78 mu=arma::linspace<arma::vec>(dmu,mumax,Nmu);
79 nu=arma::linspace<arma::vec>(0.5*dnu,(Nnu-0.5)*dnu,Nnu);
80 } else {
81 mu=arma::linspace<arma::vec>(0.0,mumax,Nmu);
82 nu=arma::linspace<arma::vec>(0.0,M_PI,Nnu);
83 }
84 double dmudnu=(mu(1)-mu(0))*(nu(1)-nu(0));
85
86 // Density arrays
87 arma::mat dena(Nmu,Nnu);
88 dena.zeros();
89 arma::mat denb(Nmu,Nnu);
90 denb.zeros();
91 arma::mat dV(Nmu,Nnu);
92 dV.zeros();
93 for(size_t inu=0;inu<nu.n_elem;inu++)
94 for(size_t imu=0;imu<mu.n_elem;imu++) {
95 // Evaluate basis functions
96 double muval(mu(imu));
97 double cosnu(cos(nu(inu)));
98 arma::cx_vec bf(basis.eval_bf(muval, cosnu, phi));
99
100 // Evaluate density at the point
101 dena(imu,inu)=std::real(arma::as_scalar(bf.t()*Pa*bf));
102 denb(imu,inu)=std::real(arma::as_scalar(bf.t()*Pb*bf));
103
104 // Volume element is
105 double shmu(sinh(muval));
106 double chmu(cosh(muval));
107 double sinnu(sin(nu(inu)));
108 // Get 2 pi from phi integral (assuming density is symmetric)
109 dV(imu,inu)=2.0*M_PI*std::pow(basis.get_Rhalf(),3)*shmu*(chmu*chmu - cosnu*cosnu)*sinnu*dmudnu;
110 }
111
112 // Total density
113 arma::mat den(dena+denb);
114
115 printf("Norm of Pa on grid is %e\n",arma::sum(arma::sum(dena%dV)));
116 printf("Norm of Pb on grid is %e\n",arma::sum(arma::sum(denb%dV)));
117 printf("Norm of P on grid is %e\n",arma::sum(arma::sum(den%dV)));
118
119 Checkpoint savechk(savedens,true);
120 savechk.write("mu",mu);
121 savechk.write("nu",nu);
122 savechk.write("P",den);
123 savechk.write("Pa",dena);
124 savechk.write("Pb",denb);
125 savechk.write("R",2*basis.get_Rhalf());
126 printf("Saved density to file %s\n",savedens.c_str());
127
128 return 0;
129 }
130