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/constants.h"
18 #include "../general/timer.h"
19 #include "../general/elements.h"
20 #include "../general/scf_helpers.h"
21 #include "../general/model_potential.h"
22 #include "utils.h"
23 #include "../atomic/basis.h"
24 #include "basis.h"
25 #include "twodquadrature.h"
26 #include <cfloat>
27
28 using namespace helfem;
29
eval(int Z1,int Z2,double Rrms1,double Rrms2,double Rbond,const polynomial_basis::PolynomialBasis * poly,int Nquad,int Nelem,double Rmax,const arma::ivec & lmmax,int igrid,double zexp,double Ez,double Qzz,double Bz,int norb,double & E,arma::uword & nang,arma::uword & nrad,arma::vec & Eval,int imodel)30 void eval(int Z1, int Z2, double Rrms1, double Rrms2, double Rbond, const polynomial_basis::PolynomialBasis * poly, int Nquad, int Nelem, double Rmax, const arma::ivec & lmmax, int igrid, double zexp, double Ez, double Qzz, double Bz, int norb, double & E, arma::uword & nang, arma::uword & nrad, arma::vec & Eval, int imodel) {
31
32 int lpad=0;
33 int symm=1;
34
35 double Rhalf=0.5*Rbond;
36 double mumax(utils::arcosh(Rmax/Rhalf));
37 arma::vec bval(atomic::basis::normal_grid(Nelem, mumax, igrid, zexp));
38
39 arma::ivec lval, mval;
40 diatomic::basis::lm_to_l_m(lmmax,lval,mval);
41
42 diatomic::basis::TwoDBasis basis(Z1, Z2, Rbond, poly, Nquad, bval, lval, mval, lpad, false);
43
44 bool diag=true;
45 // Symmetry indices
46 std::vector<arma::uvec> dsym=basis.get_sym_idx(symm);
47
48 // Form overlap matrix
49 arma::mat S(basis.overlap());
50 // Form kinetic energy matrix
51 arma::mat T(basis.kinetic());
52 // Get half-inverse
53 arma::mat Sinvh(basis.Sinvh(!diag,symm));
54 // Form nuclear attraction energy matrix
55 arma::mat Vnuc;
56
57 if(imodel==0) {
58 Vnuc=basis.nuclear();
59 } else {
60 modelpotential::ModelPotential * p1, * p2;
61 if(imodel == 1) {
62 // Use GSZ guess
63 p1 = new modelpotential::GSZAtom(Z1);
64 p2 = new modelpotential::GSZAtom(Z2);
65 } else if(imodel == 2) {
66 // Use SAP guess
67 p1 = new modelpotential::SAPAtom(Z1);
68 p2 = new modelpotential::SAPAtom(Z2);
69 } else if(imodel == 3) {
70 // Use Thomas-Fermi guess
71 p1 = new modelpotential::TFAtom(Z1);
72 p2 = new modelpotential::TFAtom(Z2);
73 } else if(imodel < modelpotential::NOSUCH_NUCLEUS+4) {
74 p1 = modelpotential::get_nuclear_model((modelpotential::nuclear_model_t) (imodel-4),Z1,Rrms1);
75 p2 = modelpotential::get_nuclear_model((modelpotential::nuclear_model_t) (imodel-4),Z2,Rrms2);
76 } else {
77 throw std::logic_error("Unsupported guess\n");
78 }
79
80 // Form model integral
81 int lquad = 4*arma::max(lmmax)+12;
82 helfem::diatomic::twodquad::TwoDGrid qgrid;
83 qgrid=helfem::diatomic::twodquad::TwoDGrid(&basis,lquad);
84 Vnuc=qgrid.model_potential(p1, p2);
85 delete p1;
86 delete p2;
87 }
88
89 // Form Hamiltonian
90 arma::mat H0(T+Vnuc);
91 if(Ez!=0.0)
92 H0+=Ez*basis.dipole_z();
93 // Quadrupole coupling
94 if(Qzz!=0.0)
95 H0+=Qzz*basis.quadrupole_zz()/3.0;
96 // Magnetic field coupling
97 if(Bz!=0.0) {
98 printf("Bz=%e\n",Bz);
99 H0+=basis.Bz_field(Bz)-Bz*S/2.0;
100 }
101 // Use core guess
102 arma::vec Ev;
103 arma::mat C;
104 scf::eig_gsym_sub(Ev,C,H0,Sinvh,dsym);
105 //scf::eig_gsym(Ev,C,H0,Sinvh);
106
107 // Sanity check
108 norb=std::min((int) Ev.n_elem,norb);
109 Eval=Ev.subvec(0,norb-1);
110 E=arma::sum(Eval);
111 nang=basis.Nang();
112 nrad=basis.Nrad();
113 }
114
main(int argc,char ** argv)115 int main(int argc, char **argv) {
116 cmdline::parser parser;
117
118 // full option name, no short option, description, argument required
119 parser.add<std::string>("Z1", 0, "first nuclear charge", true);
120 parser.add<std::string>("Z2", 0, "second nuclear charge", true);
121 parser.add<double>("Rrms1", 0, "atom 1 rms size", false, 0.0);
122 parser.add<double>("Rrms2", 0, "atom 2 rms size", false, 0.0);
123 parser.add<double>("Rbond", 0, "internuclear distance", true);
124 parser.add<bool>("angstrom", 0, "input distances in angstrom", false, false);
125 parser.add<double>("Rmax", 0, "practical infinity in au", false, 40.0);
126 parser.add<int>("grid", 0, "type of grid: 1 for linear, 2 for quadratic, 3 for polynomial, 4 for exponential", false, 4);
127 parser.add<double>("zexp", 0, "parameter in radial grid", false, 1.0);
128 parser.add<int>("nnodes", 0, "number of nodes per element", false, 15);
129 parser.add<int>("primbas", 0, "primitive radial basis", false, 4);
130 parser.add<int>("nquad", 0, "number of quadrature points", false, 0);
131 parser.add<double>("Ez", 0, "electric dipole field", false, 0.0);
132 parser.add<double>("Qzz", 0, "electric quadrupole field", false, 0.0);
133 parser.add<double>("Bz", 0, "magnetic dipole field", false, 0.0);
134 parser.add<int>("thresh", 0, "convergence threshold, 10 corresponds to 1e-10", false, 10);
135 parser.add<int>("nadd", 0, "number of funcs to add", false, 2);
136 parser.add<int>("imodel", 0, "model potential: bare nucleus (0), GSZ (1), SAP (2)", false, 0);
137 parser.parse_check(argc, argv);
138
139 // Get parameters
140 double Rmax(parser.get<double>("Rmax"));
141 int igrid(parser.get<int>("grid"));
142 double zexp(parser.get<double>("zexp"));
143 double Ez(parser.get<double>("Ez"));
144 double Qzz(parser.get<double>("Qzz"));
145 double Bz(parser.get<double>("Bz"));
146 int primbas(parser.get<int>("primbas"));
147 // Number of nodes
148 int Nnodes(parser.get<int>("nnodes"));
149 // Order of quadrature rule
150 int Nquad(parser.get<int>("nquad"));
151 int nadd(parser.get<int>("nadd"));
152 int thresh(parser.get<int>("thresh"));
153 int imodel(parser.get<int>("imodel"));
154
155 if(nadd%2)
156 printf("WARNING - Adding an odd number of functions at a time does not give a balanced description of gerade/ungerade orbitals and may give wrong results.\n");
157
158 // Nuclear charge
159 int Z1(get_Z(parser.get<std::string>("Z1")));
160 int Z2(get_Z(parser.get<std::string>("Z2")));
161 double Rrms1(parser.get<double>("Rrms1"));
162 double Rrms2(parser.get<double>("Rrms2"));
163 double Rbond(parser.get<double>("Rbond"));
164
165 if(parser.get<bool>("angstrom")) {
166 // Convert to atomic units
167 Rbond*=ANGSTROMINBOHR;
168 }
169
170 // Determine orbitals
171 std::vector<int> norbs(4);
172 num_orbs(Z1, Z2, norbs[0], norbs[1], norbs[2], norbs[3]);
173 for(size_t m=norbs.size()-1;m<norbs.size();m--)
174 if(norbs[m]==0)
175 norbs.erase(norbs.begin()+m);
176
177 printf("Determining basis set for %s-%s at distance %e with Rmax=%e.\n",element_symbols[Z1].c_str(),element_symbols[Z2].c_str(),Rbond,Rmax);
178
179 // Get primitive basis
180 polynomial_basis::PolynomialBasis *poly(polynomial_basis::get_basis(primbas,Nnodes));
181
182 if(Nquad==0)
183 // Set default value
184 Nquad=5*Nnodes;
185 else if(Nquad<2*Nnodes)
186 throw std::logic_error("Insufficient radial quadrature.\n");
187
188 printf("Using %i point quadrature rule.\n",Nquad);
189
190 // Determine number of elements
191 int Nelem=1;
192
193 // Final angular grid
194 arma::ivec lmgrid(norbs.size());
195
196 // Loop over threshold
197 int ithr=0;
198 int othr=-1;
199 std::vector<bool> init(norbs.size(),true);
200
201 while(ithr<=thresh) {
202 double thr=1.0/std::pow(10.0,ithr);
203
204 if(ithr!=othr) {
205 printf("**** thr = %e ****\n",thr);
206 othr=ithr;
207 }
208 bool cvd=true;
209
210 // Loop over orbital types
211 for(size_t m=norbs.size()-1;m<norbs.size();m--) {
212 // Angular grid in test calculation
213 arma::ivec lmmax;
214 lmmax.ones(m+1);
215 lmmax*=-1;
216 // Start off value
217 if(init[m]) {
218 if(m<norbs.size()-1) {
219 // Safe assumption: n(sigma) > n(pi) > n(delta) > n(phi)
220 lmmax(m)=lmgrid(m+1);
221 } else {
222 // Start off with one function
223 lmmax(m)=m;
224 }
225 init[m]=false;
226 } else {
227 // Use previous value
228 lmmax(m)=lmgrid(m);
229 }
230
231 // Initial estimate
232 double E;
233 arma::vec Eval;
234 arma::uword Nrad, Nang;
235 eval(Z1, Z2, Rrms1, Rrms2, Rbond, poly, Nquad, Nelem, Rmax, lmmax, igrid, zexp, Ez, Qzz, Bz, norbs[m], E, Nrad, Nang, Eval, imodel);
236
237 Eval.t().print("Initial eigenvalues");
238
239 printf("Initial energy is %e\n",E);
240
241 int iiter=0;
242 double dE;
243 while(true) {
244 iiter++;
245
246 double Ea, Er;
247 arma::vec Eva, Evr;
248 arma::uword Nra, Naa;
249 arma::uword Nrr, Nar;
250
251 arma::ivec lmtr(lmmax);
252 lmtr(m)+=nadd;
253
254 printf("m=%i iteration %i\n",(int) m,iiter);
255
256 eval(Z1, Z2, Rrms1, Rrms2, Rbond, poly, Nquad, Nelem, Rmax, lmtr, igrid, zexp, Ez, Qzz, Bz, norbs[m], Ea, Nra, Naa, Eva, imodel);
257 double dEa=Ea-E;
258 printf("Addition of %i partial waves decreases energy by %e\n",nadd,dEa);
259
260 eval(Z1, Z2, Rrms1, Rrms2, Rbond, poly, Nquad, Nelem+nadd, Rmax, lmmax, igrid, zexp, Ez, Qzz, Bz, norbs[m], Er, Nrr, Nar, Evr, imodel);
261 double dEr=Er-E;
262 printf("Addition of %i radial elements decreases energy by %e\n",nadd,dEr);
263
264 dE=std::min(dEa,dEr);
265 if(dE>-thr)
266 break;
267
268 // Angular loop is not converged
269 cvd=false;
270
271 if(dEa==dE) {
272 lmmax=lmtr;
273 E=Ea;
274 Eval=Eva;
275 Nrad=Nra;
276 Nang=Naa;
277 printf("Basis set has now %i partial waves\n",(int) lmmax(m));
278 } else {
279 Nelem+=nadd;
280 E=Er;
281 Eval=Evr;
282 Nrad=Nrr;
283 Nang=Nar;
284 printf("Basis set has now %i radial elements\n",Nelem);
285 }
286 Eval.t().print("Current eigenvalues");
287 printf("\n");
288 }
289 printf("m=%i is converged with %i elements and %i partial waves\n\n",(int) m,(int) Nelem,(int) lmmax(m));
290 lmgrid(m)=lmmax(m);
291 }
292
293 if(cvd) {
294 printf("\n");
295 printf("An estimated accuracy of %e is achieved with\n",thr);
296 printf("--Z1=%s --Z2=%s --Rbond=%e --angstrom=%i --grid=%i --zexp=%e --primbas=%i --nnodes=%i --nelem=%i --Rmax=%e --lmax=", parser.get<std::string>("Z1").c_str(), parser.get<std::string>("Z2").c_str(), parser.get<double>("Rbond"), parser.get<bool>("angstrom"), igrid, zexp, primbas, Nnodes, (int) Nelem, Rmax);
297
298 for(size_t m=0;m<lmgrid.n_elem;m++) {
299 if(m)
300 printf(",");
301 printf("%i",(int) lmgrid(m));
302 }
303 printf("\n\n\n");
304 ithr++;
305 }
306 }
307
308 return 0;
309 }
310