1 /*
2 * This source code is part of
3 *
4 * E R K A L E
5 * -
6 * HF/DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2012
9 * Copyright (c) 2010-2012, 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
17 #include "emd.h"
18 #include "emd_gto.h"
19 #include "emd_similarity.h"
20 #include "emdcube.h"
21 //include "gto_fourier_ylm.h"
22 #include "spherical_expansion.h"
23 #include "../solidharmonics.h"
24 #include "../spherical_harmonics.h"
25 #include "../checkpoint.h"
26 #include "../settings.h"
27 #include "../stringutil.h"
28 #include "../timer.h"
29
30 #include <armadillo>
31 #include <cstdio>
32 #include <cstdlib>
33 #include <cfloat>
34
35 #ifdef _OPENMP
36 #include <omp.h>
37 #endif
38
39 #ifdef SVNRELEASE
40 #include "../version.h"
41 #endif
42
43 Settings settings;
44
main_guarded(int argc,char ** argv)45 int main_guarded(int argc, char **argv) {
46 #ifdef _OPENMP
47 printf("ERKALE - EMD from Hel, OpenMP version, running on %i cores.\n",omp_get_max_threads());
48 #else
49 printf("ERKALE - EMD from Hel, serial version.\n");
50 #endif
51 print_copyright();
52 print_license();
53 #ifdef SVNRELEASE
54 printf("At svn revision %s.\n\n",SVNREVISION);
55 #endif
56 print_hostname();
57
58 if(argc!=1 && argc!=2) {
59 printf("Usage: $ %s (runfile)\n",argv[0]);
60 return 0;
61 }
62
63 Timer t;
64
65 // Parse settings
66 settings.add_string("LoadChk","Checkpoint file to load density from","erkale.chk");
67 settings.add_bool("DoEMD", "Perform calculation of isotropic EMD (moments of EMD, Compton profile)", true);
68 settings.add_double("EMDTol", "Tolerance for the numerical integration of the radial EMD",1e-8);
69 settings.add_string("EMDlm", "Which projection of the radial EMD to compute","");
70 settings.add_bool("EMDAdapt", "Use adaptive grid to compute EMD?", true);
71 settings.add_string("EMDCube", "Calculate EMD on a cube? e.g. -10:.3:10 -5:.2:4 -2:.1:3", "");
72 settings.add_string("EMDOrbitals", "Compute EMD of given orbitals, e.g. 1,2,4:6","");
73 settings.add_string("Similarity", "Compute similarity measure to checkpoint","");
74 settings.add_string("SimilarityGrid", "Grid to use for computing similarity integrals","500 77");
75 settings.add_bool("SimilarityLM", "Seminumerical computation of similarity integrals?", false);
76 settings.add_int("SimilarityLmax", "Maximum angular momentum for seminumerical computation", 6);
77
78 if(argc==2)
79 settings.parse(argv[1]);
80 else
81 printf("Using default settings.\n");
82 settings.print();
83
84 // Get the tolerance
85 double tol=settings.get_double("EMDTol");
86
87 // Load checkpoint
88 Checkpoint chkpt(settings.get_string("LoadChk"),false);
89
90 // Load basis set
91 BasisSet basis;
92 chkpt.read(basis);
93
94 // Load density matrix
95 arma::cx_mat P;
96 arma::mat Pr, Pi;
97 chkpt.read("P",Pr);
98 if(chkpt.exist("P_im")) {
99 chkpt.read("P_im",Pi);
100 P=Pr*COMPLEX1 + Pi*COMPLEXI;
101 } else
102 P=Pr*COMPLEX1;
103
104 // The projection to calculate
105 int l=0, m=0;
106 std::string lmstr=settings.get_string("EMDlm");
107 if(lmstr.size()) {
108 // Get l and m values
109 std::vector<std::string> lmval=splitline(lmstr);
110 if(lmval.size()!=2)
111 throw std::runtime_error("Invalid specification of l and m values.\n");
112 l=readint(lmval[0]);
113 m=readint(lmval[1]);
114 }
115 bool adaptive=settings.get_bool("EMDAdapt");
116
117 // Compute orbital EMDs?
118 if(settings.get_string("EMDOrbitals")!="") {
119 // Get orbitals
120 std::vector<std::string> orbs=splitline(settings.get_string("EMDOrbitals"));
121
122 // Polarized calculation?
123 bool restr;
124 chkpt.read("Restricted",restr);
125 if(restr!= (orbs.size()==1))
126 throw std::runtime_error("Invalid occupancies for spin alpha and beta!\n");
127
128 if(l!=0)
129 printf("\nComputing the (%i %+i) projection of the orbital EMD.\n",l,m);
130
131 for(size_t ispin=0;ispin<orbs.size();ispin++) {
132 // Indices of orbitals to include.
133 std::vector<size_t> idx=parse_range(orbs[ispin]);
134 // Change into C++ indexing
135 for(size_t i=0;i<idx.size();i++)
136 idx[i]--;
137
138 // Read orbitals
139 arma::mat C;
140 if(restr)
141 chkpt.read("C",C);
142 else {
143 if(ispin==0)
144 chkpt.read("Ca",C);
145 else
146 chkpt.read("Cb",C);
147 }
148
149 for(size_t i=0;i<idx.size();i++) {
150 // Names of output files
151 std::ostringstream emdname;
152 std::ostringstream momname;
153 std::ostringstream Jname;
154 std::ostringstream Jintname;
155
156 std::ostringstream suffix;
157 if(l>0)
158 suffix << "_" << l << "_" << m;
159 suffix << ".txt";
160
161 if(restr) {
162 emdname << "emd-" << idx[i]+1 << suffix.str();
163 momname << "moments-" << idx[i]+1 << suffix.str();
164 Jname << "compton-" << idx[i]+1 << suffix.str();
165 Jintname << "compton-interp-" << idx[i]+1 << suffix.str();
166 } else {
167 if(ispin==0) {
168 emdname << "emd-a-" << idx[i]+1 << suffix.str();
169 momname << "moments-a-" << idx[i]+1 << suffix.str();
170 Jname << "compton-a-" << idx[i]+1 << suffix.str();
171 Jintname << "compton-interp-a-" << idx[i]+1 << suffix.str();
172 } else {
173 emdname << "emd-b-" << idx[i]+1 << suffix.str();
174 momname << "moments-b-" << idx[i]+1 << suffix.str();
175 Jname << "compton-b-" << idx[i]+1 << suffix.str();
176 Jintname << "compton-interp-b-" << idx[i]+1 << suffix.str();
177 }
178 }
179
180 // Generate dummy density matrix
181 arma::cx_mat Pdum=C.col(idx[i])*arma::trans(C.col(idx[i]))*COMPLEX1;
182
183 Timer temd;
184
185 GaussianEMDEvaluator *poseval=new GaussianEMDEvaluator(basis,Pdum,l,std::abs(m));
186 GaussianEMDEvaluator *negeval;
187 if(m!=0)
188 negeval=new GaussianEMDEvaluator(basis,Pdum,l,-std::abs(m));
189 else
190 negeval=NULL;
191
192 EMD emd(poseval, negeval, 1, l, m);
193 if(adaptive) {
194 emd.initial_fill();
195 if(l==0 && m==0) emd.find_electrons();
196 emd.optimize_moments(true,tol);
197 } else
198 emd.fixed_fill();
199
200 emd.save(emdname.str());
201 emd.moments(momname.str());
202 if(l==0 && m==0) {
203 emd.compton_profile(Jname.str());
204 emd.compton_profile_interp(Jintname.str());
205 }
206
207 delete poseval;
208 if(m!=0) delete negeval;
209 }
210 }
211 }
212
213 if(settings.get_bool("DoEMD")) {
214 t.print_time();
215
216 printf("\nCalculating EMD properties.\n");
217 printf("Please read and cite the reference:\n%s\n%s\n%s\n", \
218 "J. Lehtola, M. Hakala, J. Vaara and K. Hämäläinen", \
219 "Calculation of isotropic Compton profiles with Gaussian basis sets", \
220 "Phys. Chem. Chem. Phys. 13 (2011), pp. 5630 - 5641.");
221
222 if(l!=0)
223 printf("\nComputing the (%i %+i) projection of the EMD.\n",l,m);
224 else
225 printf("\nComputing the isotropic projection of the EMD.\n");
226
227 // Amount of electrons is
228 int Nel;
229 chkpt.read("Nel",Nel);
230
231 // Construct EMD evaluators
232 Timer temd;
233 GaussianEMDEvaluator *poseval=new GaussianEMDEvaluator(basis,P,l,std::abs(m));
234 GaussianEMDEvaluator *negeval;
235 if(m!=0)
236 negeval=new GaussianEMDEvaluator(basis,P,l,-std::abs(m));
237 else
238 negeval=NULL;
239
240 temd.set();
241 EMD emd(poseval, negeval, Nel, l, m);
242 if(adaptive) {
243 emd.initial_fill();
244 if(l==0 && m==0) emd.find_electrons();
245 emd.optimize_moments(true,tol);
246 } else
247 emd.fixed_fill();
248
249 if(l==0 && m==0) {
250 emd.save("emd.txt");
251 emd.moments("moments.txt");
252 emd.compton_profile("compton.txt");
253 emd.compton_profile_interp("compton-interp.txt");
254 } else {
255 char fname[80];
256 sprintf(fname,"emd_%i_%i.txt",l,m);
257 emd.save(fname);
258
259 sprintf(fname,"moments_%i_%i.txt",l,m);
260 emd.moments(fname);
261 }
262
263 if(l==0 && m==0)
264 printf("Calculating isotropic EMD properties took %s.\n",temd.elapsed().c_str());
265 else
266 printf("Calculating projected EMD properties took %s.\n",temd.elapsed().c_str());
267
268 delete poseval;
269 if(m!=0) delete negeval;
270 }
271
272 // Do EMD on a cube?
273 if(stricmp(settings.get_string("EMDCube"),"")!=0) {
274 t.print_time();
275 Timer temd;
276
277 // Form grid in p space.
278 std::vector<double> px, py, pz;
279 parse_cube(settings.get_string("EMDCube"),px,py,pz);
280
281 // Calculate EMD on cube
282 emd_cube(basis,P,px,py,pz);
283
284 printf("Calculating EMD on a cube took %s.\n",temd.elapsed().c_str());
285 }
286
287 // Compute similarity?
288 if(stricmp(settings.get_string("Similarity"),"")!=0) {
289
290 // Load checkpoint
291 Checkpoint simchk(settings.get_string("Similarity"),false);
292
293 // Get grid size
294 std::vector<std::string> gridsize=splitline(settings.get_string("SimilarityGrid"));
295 if(gridsize.size()!=2) {
296 throw std::runtime_error("Invalid grid size!\n");
297 }
298 int nrad=readint(gridsize[0]);
299 int lmax=readint(gridsize[1]);
300 int radlmax=settings.get_int("SimilarityLmax");
301
302 // Load basis set
303 BasisSet simbas;
304 simchk.read(simbas);
305
306 // Load density matrix
307 arma::mat simPr, simPi;
308 arma::cx_mat simP;
309 simchk.read("P",simPr);
310 if(simchk.exist("P_im")) {
311 simchk.read("P_im",simPi);
312 simP=simPr*COMPLEX1 + simPi*COMPLEXI;
313 } else
314 simP=simPr*COMPLEX1;
315
316 // Compute momentum density overlap
317 arma::cube ovl;
318 if(settings.get_bool("SimilarityLM"))
319 ovl=emd_overlap_semi(basis,P,simbas,simP,nrad,radlmax);
320 else
321 ovl=emd_overlap(basis,P,simbas,simP,nrad,lmax);
322
323 // Amount of electrons
324 int Nela, Nelb;
325 chkpt.read("Nel", Nela);
326 simchk.read("Nel", Nelb);
327
328 // Shape function overlap
329 arma::cube sh=emd_similarity(ovl,Nela,Nelb);
330 sh.slice(0).save("similarity.dat",arma::raw_ascii);
331 sh.slice(1).save("similarity_avg.dat",arma::raw_ascii);
332
333 for(int s=0;s<2;s++) {
334 if(s) {
335 printf("%2s\t%9s\t%9s\t%9s\t%9s\t%9s\t%9s\t%9s\n","k","S0(AA)","S0(BB)","S0(AB)","I0(AA)","I0(BB)","I0(AB)","D0(AB)");
336 for(int k=-1;k<3;k++)
337 // Vandenbussche don't include p^2 in the spherical average
338 printf("%2i\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", k+1, sh(k+1,0,s), sh(k+1,1,s), sh(k+1,2,s), sh(k+1,3,s), sh(k+1,4,s), sh(k+1,5,s), sh(k+1,6,s));
339 } else {
340 printf("%2s\t%9s\t%9s\t%9s\t%9s\t%9s\t%9s\t%9s\n","k","S(AA)","S(BB)","S(AB)","I(AA)","I(BB)","I(AB)","D(AB)");
341 for(int k=-1;k<3;k++)
342 printf("%2i\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", k, sh(k+1,0,s), sh(k+1,1,s), sh(k+1,2,s), sh(k+1,3,s), sh(k+1,4,s), sh(k+1,5,s), sh(k+1,6,s));
343 }
344 printf("\n");
345 }
346 }
347
348 return 0;
349 }
350
main(int argc,char ** argv)351 int main(int argc, char **argv) {
352 try {
353 return main_guarded(argc, argv);
354 } catch (const std::exception &e) {
355 std::cerr << "error: " << e.what() << std::endl;
356 return 1;
357 }
358 }
359