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-2011
9 * Copyright (c) 2010-2011, 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 "bfprod.h"
18 #include "fourierprod.h"
19 #include "lmtrans.h"
20 #include "momentum_series.h"
21 #include "../basis.h"
22 #include "../basislibrary.h"
23 #include "../dftgrid.h"
24 #include "../dftfuncs.h"
25 #include "../elements.h"
26 #include "../density_fitting.h"
27 #include "../linalg.h"
28 #include "../lmgrid.h"
29 #include "../mathf.h"
30 #include "../scf.h"
31 #include "../settings.h"
32 #include "../stringutil.h"
33 #include "../tempered.h"
34 #include "../timer.h"
35 #include "../xyzutils.h"
36 #include "xrsscf.h"
37
38 // Needed for libint init
39 #include "../eriworker.h"
40
41 #include <algorithm>
42 #include <armadillo>
43 #include <cstdio>
44 #include <cstdlib>
45 #include <cfloat>
46 #include <istream>
47
48 #ifdef _OPENMP
49 #include <omp.h>
50 #endif
51
52 #ifdef SVNRELEASE
53 #include "../version.h"
54 #endif
55
56 Settings settings;
57
58
59 /**
60 * Was loading a success?
61 *
62 * LOAD_FAIL: Failed to load. (start from scratch)
63 * LOAD_SUCC: Succesfully loaded completed calculation. (no calculation needed)
64 * LOAD_DIFF: Succesfully loaded calculation of a different type. (use as starting point)
65 */
66 enum loadresult {LOAD_FAIL, LOAD_SUCC, LOAD_DIFF};
67
68
parse_method(const std::string & method)69 enum xrs_method parse_method(const std::string & method) {
70 enum xrs_method met;
71 if(stricmp(method,"TP")==0)
72 met=TP;
73 else if(stricmp(method,"FCH")==0)
74 met=FCH;
75 else if(stricmp(method,"XCH")==0)
76 met=XCH;
77 else throw std::runtime_error("Unrecognized method.\n");
78
79 return met;
80 }
81
82
83 /// Augment the basis set with diffuse functions
augment_basis(const BasisSet & basis)84 BasisSet augment_basis(const BasisSet & basis) {
85 // Get indices of atoms to augment
86 std::vector<size_t> augind=parse_range(splitline(settings.get_string("XRSAugment"))[0]);
87 // Convert to C++ indexing
88 for(size_t i=0;i<augind.size();i++) {
89 if(augind[i]==0)
90 throw std::runtime_error("Error - nuclear index must be positive!\n");
91 augind[i]--;
92 }
93
94 bool verbose=settings.get_bool("Verbose");
95 if(verbose) {
96 printf("\nAugmenting basis set with diffuse functions.\n");
97 fflush(stdout);
98 }
99
100 // Form augmented basis
101 BasisSet augbas(basis);
102
103 // Basis set for augmentation functions
104 BasisSetLibrary augbaslib;
105 augbaslib.load_basis(settings.get_string("XRSDoubleBasis"));
106
107 // Loop over excited atoms
108 for(size_t iaug=0;iaug<augind.size();iaug++) {
109 // Index of atom is
110 const size_t ind=augind[iaug];
111
112 // The symbol of the atom
113 std::string el=basis.get_symbol(ind);
114
115 // The basis to use for the atom.
116 ElementBasisSet elbas;
117
118 try {
119 // Check first if a special set is wanted for given center
120 elbas=augbaslib.get_element(el,ind+1);
121 } catch(std::runtime_error & err) {
122 // Did not find a special basis, use the general one instead.
123 elbas=augbaslib.get_element(el,0);
124 }
125
126 // Get original number of shells
127 size_t Nsh_orig=augbas.get_Nshells();
128 // Add shells, no sorting.
129 augbas.add_shells(ind,elbas,false);
130 // Convert contractions on the added shells
131 for(size_t ish=Nsh_orig;ish<augbas.get_Nshells();ish++)
132 augbas.convert_contraction(ish);
133 }
134
135 // Finalize augmentation basis
136 augbas.finalize();
137
138 return augbas;
139 }
140
141 /**
142 * Double basis set method
143 *
144 * Augment the basis with diffuse functions and diagonalize the Fock
145 * matrix in the unoccupied space. The occupied orbitals and their
146 * energies stay the same in the approximation.
147 */
augmented_solution(const BasisSet & basis,const uscf_t & sol,size_t nocca,size_t noccb,dft_t dft,BasisSet & augbas,arma::mat & Caug,arma::vec & Eaug,bool spin,enum xrs_method method)148 void augmented_solution(const BasisSet & basis, const uscf_t & sol, size_t nocca, size_t noccb, dft_t dft, BasisSet & augbas, arma::mat & Caug, arma::vec & Eaug, bool spin, enum xrs_method method) {
149 Timer ttot;
150
151 augbas=augment_basis(basis);
152 // Need to update pointers in augbas
153 augbas.update_nuclear_shell_list();
154
155
156 // Amount of functions in original basis set is
157 const size_t Nbf=basis.get_Nbf();
158 // Total number of functions in augmented set is
159 const size_t Ntot=augbas.get_Nbf();
160 // Amount of augmentation functions is
161 const size_t Naug=Ntot-Nbf;
162
163 bool verbose=settings.get_bool("Verbose");
164
165 if(verbose) {
166 printf("\nAugmented original basis (%i functions) with %i diffuse functions.\n",(int) Nbf,(int) Naug);
167 printf("Computing unoccupied orbitals.\n");
168 fflush(stdout);
169
170 fprintf(stderr,"Calculating unoccupied orbitals in augmented basis.\n");
171 fflush(stderr);
172 }
173
174
175 // Augmented solution
176 uscf_t augsol(sol);
177
178 // Form density matrix.
179 arma::mat Paaug(Ntot,Ntot), Pbaug(Ntot,Ntot), Paug(Ntot,Ntot);
180
181 augsol.Pa.zeros(Ntot,Ntot);
182 augsol.Pa.submat(0,0,Nbf-1,Nbf-1)=sol.Pa;
183
184 augsol.Pb.zeros(Ntot,Ntot);
185 augsol.Pb.submat(0,0,Nbf-1,Nbf-1)=sol.Pb;
186
187 augsol.P=augsol.Pa+augsol.Pb;
188
189 augsol.Ca=project_orbitals(sol.Ca,basis,augbas);
190 augsol.Cb=project_orbitals(sol.Cb,basis,augbas);
191
192 // Checkpoint
193 std::string augchk=settings.get_string("AugChk");
194 bool delchk=false;
195 if(stricmp(augchk,"")==0) {
196 delchk=true;
197 augchk=tempname();
198 }
199
200 {
201 // Copy base checkpoint to augmented checkpoint
202 std::string chk=settings.get_string("SaveChk");
203 char cmd[chk.size()+augchk.size()+5];
204 sprintf(cmd,"cp %s %s",chk.c_str(),augchk.c_str());
205 int cperr=system(cmd);
206 if(cperr) {
207 ERROR_INFO();
208 throw std::runtime_error("Error copying checkpoint file.\n");
209 }
210 }
211
212 Timer taug;
213
214 // Open augmented checkpoint file in write mode, don't truncate it!
215 Checkpoint chkpt(augchk,true,false);
216 // Update the basis set
217 chkpt.write(augbas);
218
219 {
220 // Initialize solver
221 XRSSCF solver(augbas,chkpt,spin);
222 // The fitting basis set from the non-augmented basis is enough,
223 // since the density is restricted there.
224 BasisSet fitbas=basis.density_fitting();
225 solver.set_fitting(fitbas);
226
227 // Load occupation number vectors
228 std::vector<double> occa, occb;
229 chkpt.read("occa",occa);
230 chkpt.read("occb",occb);
231
232 // Construct Fock matrix
233 DFTGrid grid(&augbas,verbose,dft.lobatto);
234 DFTGrid nlgrid(&augbas,verbose,dft.lobatto);
235 if(dft.adaptive) {
236 // Form DFT quadrature grid
237 grid.construct(augsol.Pa,augsol.Pb,dft.gridtol,dft.x_func,dft.c_func);
238 } else {
239 // Fixed size grid
240 grid.construct(dft.nrad,dft.lmax,dft.x_func,dft.c_func);
241 if(dft.nl)
242 nlgrid.construct(dft.nlnrad,dft.nllmax,true,false,false,true);
243 }
244
245 // Range separated functional?
246 double omega, kfull, kshort;
247 range_separation(dft.x_func,omega,kfull,kshort);
248 // Evaluators for range separated part
249 if(omega!=0.0)
250 solver.fill_rs(omega);
251
252 switch(method) {
253 case(TP):
254 solver.Fock_half_hole(augsol,dft,occa,occb,grid,nlgrid);
255 break;
256
257 case(FCH):
258 case(XCH):
259 solver.Fock_full_hole(augsol,dft,occa,occb,grid,nlgrid,method==XCH);
260 }
261 }
262
263 if(verbose) {
264 printf("Fock operator constructed in augmented space in %s.\n",taug.elapsed().c_str());
265 fprintf(stderr,"Constructed augmented Fock operator (%s).\n",taug.elapsed().c_str());
266 }
267
268 // Diagonalize unoccupied space. Loop over spin
269 for(size_t ispin=0;ispin<2;ispin++) {
270
271 // Form Fock operator
272 arma::mat H;
273 if(!ispin)
274 H=augsol.Ha;
275 else
276 H=augsol.Hb;
277
278 arma::mat AOtoO;
279 if(!ispin)
280 AOtoO=project_orbitals(sol.Ca,basis,augbas);
281 else
282 AOtoO=project_orbitals(sol.Cb,basis,augbas);
283
284 // Amount of occupied orbitals
285 size_t nocc;
286 if(!ispin)
287 nocc=nocca;
288 else
289 nocc=noccb;
290
291 // Convert Fock operator to unoccupied MO basis.
292 taug.set();
293 arma::mat H_MO=arma::trans(AOtoO.cols(nocc,AOtoO.n_cols-1))*H*AOtoO.cols(nocc,AOtoO.n_cols-1);
294 if(verbose) {
295 printf("H_MO formed in %s.\n",taug.elapsed().c_str());
296 fflush(stdout);
297 }
298
299 // Diagonalize Fockian to find orbitals and energies
300 taug.set();
301 arma::vec Eval;
302 arma::mat Evec;
303 eig_sym_ordered(Eval,Evec,H_MO);
304
305 if(verbose) {
306 printf("H_MO diagonalized in unoccupied subspace in %s.\n",taug.elapsed().c_str());
307 fflush(stdout);
308 }
309
310 // Store energies
311 arma::vec Ea;
312 Ea.zeros(AOtoO.n_cols);
313 // Occupied orbital energies
314 if(!spin)
315 Ea.subvec(0,nocc-1)=sol.Ea.subvec(0,nocc-1);
316 else
317 Ea.subvec(0,nocc-1)=sol.Eb.subvec(0,nocc-1);
318 // Virtuals
319 Ea.subvec(nocc,AOtoO.n_cols-1)=Eval;
320
321 // Back-transform orbitals to AO basis
322 arma::mat Ca;
323 Ca.zeros(Ntot,AOtoO.n_cols);
324 // Occupied orbitals, padded with zeros
325 Ca.submat(0,0,Nbf-1,nocc-1)=AOtoO.submat(0,0,Nbf-1,nocc-1);
326 // Unoccupied orbitals
327 Ca.cols(nocc,AOtoO.n_cols-1)=AOtoO.cols(nocc,AOtoO.n_cols-1)*Evec;
328
329 // Save results
330 if(!ispin) {
331 augsol.Ca=Ca;
332 augsol.Ea=Ea;
333 } else {
334 augsol.Cb=Ca;
335 augsol.Eb=Ea;
336 }
337
338 // Return variables
339 if((bool) ispin==spin) {
340 Eaug=Ea;
341 Caug=Ca;
342 }
343 }
344
345 // Write out updated solution vectors
346 chkpt.write("Ca",augsol.Ca);
347 chkpt.write("Cb",augsol.Cb);
348 chkpt.write("Ea",augsol.Ea);
349 chkpt.write("Eb",augsol.Eb);
350 // and the density matrix
351 chkpt.write("Pa",augsol.Pa);
352 chkpt.write("Pb",augsol.Pb);
353 chkpt.write("P",augsol.P);
354 // and the energy
355 chkpt.write(augsol.en);
356
357 if(delchk) {
358 int err=remove(augchk.c_str());
359 if(err) {
360 ERROR_INFO();
361 throw std::runtime_error("Error removing temporary file.\n");
362 }
363 }
364
365 if(verbose) {
366 printf("Augmentation took %s.\n",ttot.elapsed().c_str());
367 fprintf(stderr,"Augmentation done (%s).\n",taug.elapsed().c_str());
368 fflush(stdout);
369 fflush(stderr);
370 }
371 }
372
373 typedef struct {
374 // Transition energy
375 double E;
376 // Total oscillator strength
377 double w;
378 // Decomposition of strength
379 std::vector<double> wdec;
380 } spectrum_t;
381
382 /// Compute transitions to unoccupied states.
compute_transitions(const BasisSet & basis,const arma::mat & C,const arma::vec & E,size_t iat,size_t ixc,size_t nocc)383 std::vector<spectrum_t> compute_transitions(const BasisSet & basis, const arma::mat & C, const arma::vec & E, size_t iat, size_t ixc, size_t nocc) {
384 // Returned array
385 std::vector<spectrum_t> ret;
386
387 // Coordinates of excited atom
388 coords_t xccen=basis.get_nuclear_coords(iat);
389 // Dipole moment matrix
390 std::vector<arma::mat> mom1=basis.moment(1,xccen.x,xccen.y,xccen.z);
391 // Compute RHS of transition
392 arma::vec rhs_x=mom1[getind(1,0,0)]*C.col(ixc);
393 arma::vec rhs_y=mom1[getind(0,1,0)]*C.col(ixc);
394 arma::vec rhs_z=mom1[getind(0,0,1)]*C.col(ixc);
395
396 for(size_t ix=nocc;ix<C.n_cols;ix++) {
397 spectrum_t sp;
398 // Transition energy is
399 sp.E=(E[ix]-E[ixc]);
400
401 // Compute oscillator strengths in x, y and z
402 double wx=arma::dot(rhs_x,C.col(ix));
403 double wy=arma::dot(rhs_y,C.col(ix));
404 double wz=arma::dot(rhs_z,C.col(ix));
405
406 // Spherical average of oscillator strength in XRS is
407 sp.w=2.0/3.0*(wx*wx + wy*wy + wz*wz);
408 // Store decomposition
409 sp.wdec.push_back(wx);
410 sp.wdec.push_back(wy);
411 sp.wdec.push_back(wz);
412
413 ret.push_back(sp);
414 }
415
416 return ret;
417 }
418
compute_qdep_transitions_series(const BasisSet & basis,const arma::mat & C,const arma::vec & E,size_t ixc,size_t nocc,std::vector<double> qval)419 std::vector< std::vector<spectrum_t> > compute_qdep_transitions_series(const BasisSet & basis, const arma::mat & C, const arma::vec & E, size_t ixc, size_t nocc, std::vector<double> qval) {
420 Timer t;
421
422 // Get the grid for computing the spherical averages.
423 std::vector<angular_grid_t> grid=form_angular_grid(2*basis.get_max_am());
424 // We normalize the weights so that for purely dipolar transitions we
425 // get the same output as with using the dipole matrix.
426 for(size_t i=0;i<grid.size();i++) {
427 // Dipole integral is only wrt theta - divide off phi part.
428 grid[i].w/=2.0*M_PI;
429 }
430
431 // Amount of transitions is
432 size_t Ntrans=C.n_cols-nocc;
433
434 // Initialize return array
435 std::vector< std::vector<spectrum_t> > ret(qval.size());
436 for(size_t iq=0;iq<qval.size();iq++) {
437 ret[iq].resize(Ntrans);
438 for(size_t ix=nocc;ix<C.n_cols;ix++) {
439 // Transition energy is
440 ret[iq][ix-nocc].E=(E[ix]-E[ixc]);
441 // Transition speed is
442 ret[iq][ix-nocc].w=0.0;
443 }
444 }
445
446 // Copy the orbitals to complex form.
447 arma::cx_mat Cm(C.n_rows,Ntrans);
448 for(size_t ir=0;ir<C.n_rows;ir++)
449 for(size_t ix=nocc;ix<C.n_cols;ix++)
450 Cm(ir,ix-nocc)=C(ir,ix);
451
452 printf("Computing transitions using series method.\n");
453 arma::cx_vec tmp;
454
455 // Evaluator
456 momentum_transfer_series ser(&basis);
457
458 // Loop over q values.
459 for(size_t iq=0;iq<qval.size();iq++) {
460 printf("\tq = %8.3f ... ",qval[iq]);
461 fflush(stdout);
462 Timer tq;
463
464 // Loop over angular mesh
465 for(size_t ig=0;ig<grid.size();ig++) {
466
467 // Current q is
468 arma::vec q(3);
469 q(0)=qval[iq]*grid[ig].r.x;
470 q(1)=qval[iq]*grid[ig].r.y;
471 q(2)=qval[iq]*grid[ig].r.z;
472 // and the weight is
473 double w=grid[ig].w;
474
475 // Get momentum transfer matrix
476 arma::cx_mat momtrans=ser.get(q,10*DBL_EPSILON,10*DBL_EPSILON);
477 // Contract with initial state
478 tmp=momtrans*C.col(ixc);
479
480 // Compute matrix elements for all transitions
481 #ifdef _OPENMP
482 #pragma omp parallel for
483 #endif
484 // Loop over transitions
485 for(size_t it=0;it<Ntrans;it++) {
486 // The matrix element is
487 std::complex<double> el=arma::dot(tmp,Cm.col(it));
488 // so the transition velocity is increased by
489 ret[iq][it].w+=w*std::norm(el);
490 }
491 }
492
493 printf("done (%s)\n",tq.elapsed().c_str());
494 }
495
496 return ret;
497 }
498
499
compute_qdep_transitions_fourier(const BasisSet & basis,const arma::mat & C,const arma::vec & E,size_t ixc,size_t nocc,std::vector<double> qval)500 std::vector< std::vector<spectrum_t> > compute_qdep_transitions_fourier(const BasisSet & basis, const arma::mat & C, const arma::vec & E, size_t ixc, size_t nocc, std::vector<double> qval) \
501 {
502 Timer t;
503
504 printf("Computing transitions using Fourier method.\n");
505
506 // Form products of basis functions.
507 const size_t Nbf=basis.get_Nbf();
508
509 std::vector<prod_gaussian_3d> bfprod=compute_products(basis);
510
511 printf("Products done in %s.\n",t.elapsed().c_str());
512 t.set();
513
514 // Form Fourier transforms of the products.
515 std::vector<prod_fourier> bffour=fourier_transform(bfprod);
516 printf("Fourier transform done in %s.\n",t.elapsed().c_str());
517 t.set();
518
519 // Get the grid for computing the spherical averages.
520 std::vector<angular_grid_t> grid=form_angular_grid(2*basis.get_max_am());
521 // We normalize the weights so that for purely dipolar transitions we
522 // get the same output as with using the dipole matrix.
523 for(size_t i=0;i<grid.size();i++) {
524 // Dipole integral is only wrt theta - divide off phi part.
525 grid[i].w/=2.0*M_PI;
526 }
527
528 // Amount of transitions is
529 size_t Ntrans=C.n_cols-nocc;
530
531 // Initialize return array
532 std::vector< std::vector<spectrum_t> > ret(qval.size());
533 for(size_t iq=0;iq<qval.size();iq++) {
534 ret[iq].resize(Ntrans);
535 for(size_t ix=nocc;ix<C.n_cols;ix++) {
536 // Transition energy is
537 ret[iq][ix-nocc].E=(E[ix]-E[ixc]);
538 // Transition speed is
539 ret[iq][ix-nocc].w=0.0;
540 }
541 }
542
543 // Copy the orbitals to complex form.
544 arma::cx_mat Cm(C.n_rows,Ntrans);
545 for(size_t ir=0;ir<C.n_rows;ir++)
546 for(size_t ix=nocc;ix<C.n_cols;ix++)
547 Cm(ir,ix-nocc)=C(ir,ix);
548
549 // Loop over q values.
550 arma::cx_vec tmp;
551
552 for(size_t iq=0;iq<qval.size();iq++) {
553 printf("\tq = %8.3f ... ",qval[iq]);
554 fflush(stdout);
555 Timer tq;
556
557 // Loop over angular mesh
558 for(size_t ig=0;ig<grid.size();ig++) {
559
560 // Current q is
561 arma::vec q(3);
562 q(0)=qval[iq]*grid[ig].r.x;
563 q(1)=qval[iq]*grid[ig].r.y;
564 q(2)=qval[iq]*grid[ig].r.z;
565 // and the weight is
566 double w=grid[ig].w;
567
568 // Get momentum transfer matrix
569 arma::cx_mat momtrans=momentum_transfer(bffour,Nbf,q);
570
571 // Contract with initial state
572 tmp=momtrans*C.col(ixc);
573
574 // Compute matrix elements for all transitions
575 #ifdef _OPENMP
576 #pragma omp parallel for
577 #endif
578 // Loop over transitions
579 for(size_t it=0;it<Ntrans;it++) {
580 // The matrix element is
581 std::complex<double> el=arma::dot(tmp,Cm.col(it));
582 // so the transition velocity is increased by
583 ret[iq][it].w+=w*std::norm(el);
584 }
585 }
586
587 printf("done (%s)\n",tq.elapsed().c_str());
588 }
589
590 return ret;
591 }
592
593
compute_qdep_transitions_local(const BasisSet & basis,const arma::mat & C,const arma::vec & E,size_t iat,size_t ixc,size_t nocc,std::vector<double> q)594 std::vector< std::vector<spectrum_t> > compute_qdep_transitions_local(const BasisSet & basis, const arma::mat & C, const arma::vec & E, size_t iat, size_t ixc, size_t nocc, std::vector<double> q) {
595 // Get wanted quadrature info
596 int Nrad=settings.get_int("XRSNrad");
597 int Lmax=settings.get_int("XRSLmax");
598 int Lquad=settings.get_int("XRSLquad");
599
600 if(Lquad<Lmax) {
601 Lquad=Lmax;
602 fprintf(stderr,"Increasing the quadrature order to Lmax.\n");
603 }
604
605 // Do lm expansion of orbitals
606 lmtrans lm(C,basis,basis.get_nuclear_coords(iat),Nrad,Lmax,Lquad);
607
608 printf("\n");
609 lm.print_info();
610 printf("\n");
611
612 // Save radial distribution of excited orbital
613 lm.write_prob(ixc,"excited_orb.dat");
614
615 // Returned array
616 std::vector< std::vector<spectrum_t> > ret(q.size());
617 for(size_t iq=0;iq<q.size();iq++)
618 ret[iq].resize(C.n_cols-nocc);
619
620 // Loop over q
621 printf("Computing transitions.\n");
622 for(size_t iq=0;iq<q.size();iq++) {
623 printf("\tq = %8.3f ... ",q[iq]);
624 fflush(stdout);
625 Timer tq;
626
627 // Compute Bessel functions
628 bessel_t bes=lm.compute_bessel(q[iq]);
629
630 // Loop over transitions
631 #ifdef _OPENMP
632 #pragma omp parallel for schedule(dynamic,1)
633 #endif
634 for(size_t ix=nocc;ix<C.n_cols;ix++) {
635 Timer ttrans;
636
637 spectrum_t sp;
638 // Transition energy is
639 sp.E=(E[ix]-E[ixc]);
640
641 // Get oscillator strength
642 std::vector<double> osc=lm.transition_velocity(ixc,ix,bes);
643 // Total strength is
644 sp.w=osc[0];
645 // Store the rest
646 osc.erase(osc.begin());
647 sp.wdec=osc;
648
649 // Store result
650 ret[iq][ix-nocc]=sp;
651
652 // printf("\t%3i -> %3i (%s)\n",(int) ixc+1,(int) ix+1,ttrans.elapsed().c_str());
653 }
654
655 printf("done (%s)\n",tq.elapsed().c_str());
656 }
657
658
659 return ret;
660 }
661
save_spectrum(const std::vector<spectrum_t> & sp,const char * fname="dipole.dat")662 void save_spectrum(const std::vector<spectrum_t> & sp, const char *fname="dipole.dat") {
663 // Compute transitions
664 FILE *out=fopen(fname,"w");
665
666 for(size_t i=0;i<sp.size();i++) {
667 // Print out energy in eV oscillator strength and its components
668 fprintf(out,"%e %e",sp[i].E*HARTREEINEV,sp[i].w);
669 for(size_t j=0;j<sp[i].wdec.size();j++)
670 fprintf(out," % e",sp[i].wdec[j]);
671 fprintf(out,"\n");
672 }
673 fclose(out);
674 }
675
save_xas_spectrum(const std::vector<spectrum_t> & sp,const char * fname="dipole_xas.dat")676 void save_xas_spectrum(const std::vector<spectrum_t> & sp, const char *fname="dipole_xas.dat") {
677 // Compute transitions
678 FILE *out=fopen(fname,"w");
679
680 for(size_t i=0;i<sp.size();i++) {
681 // Print out energy in eV and transition amplitude
682 fprintf(out,"%e %e\n",sp[i].E*HARTREEINEV,sp[i].E*sp[i].w);
683 }
684 fclose(out);
685 }
686
load(const BasisSet & basis,Checkpoint & chkpt,uscf_t & sol,arma::vec & core)687 enum loadresult load(const BasisSet & basis, Checkpoint & chkpt, uscf_t & sol, arma::vec & core) {
688 // Was the load a success?
689 enum loadresult ok=LOAD_SUCC;
690
691 // Basis set used in the checkpoint file
692 BasisSet loadbas;
693
694 try {
695 chkpt.read("Ca",sol.Ca);
696 chkpt.read("Cb",sol.Cb);
697 chkpt.read("Ea",sol.Ea);
698 chkpt.read("Eb",sol.Eb);
699 chkpt.read("Pa",sol.Pa);
700 chkpt.read("Pb",sol.Pb);
701
702 chkpt.read(loadbas);
703 } catch(std::runtime_error & err) {
704 ok=LOAD_FAIL;
705 // fprintf(stderr,"Loading failed due to \"%s\".\n",err.what());
706 }
707
708 if(ok) {
709 // Check consistency
710 if(!(basis==loadbas)) {
711 ok=LOAD_FAIL;
712 // fprintf(stderr,"Basis sets differ!\n");
713 }
714 }
715
716 if(ok) {
717 // Get number of basis functions
718 size_t Nbf=basis.get_Nbf();
719
720 if(sol.Ca.n_rows != Nbf)
721 ok=LOAD_FAIL;
722
723 if(sol.Cb.n_rows != Nbf)
724 ok=LOAD_FAIL;
725
726 if(sol.Ea.n_elem != sol.Ca.n_cols)
727 ok=LOAD_FAIL;
728
729 if(sol.Eb.n_elem != sol.Cb.n_cols)
730 ok=LOAD_FAIL;
731
732 if(sol.Ea.n_elem != sol.Eb.n_elem)
733 ok=LOAD_FAIL;
734
735 if(sol.Pa.n_rows != Nbf || sol.Pa.n_cols != Nbf)
736 ok=LOAD_FAIL;
737
738 if(sol.Pb.n_rows != Nbf || sol.Pb.n_cols != Nbf)
739 ok=LOAD_FAIL;
740
741 // if(!ok) fprintf(stderr,"Dimensions do not match!\n");
742 }
743
744 // Was the calculation converged?
745 if(ok) {
746 bool conv;
747 chkpt.read("Converged",conv);
748
749 if(!conv) {
750 // fprintf(stderr,"Calculation was not converged.\n");
751 ok=LOAD_FAIL;
752 }
753 }
754
755 // Check consistency of spin
756 if(ok) {
757 bool spin;
758 chkpt.read("XRSSpin",spin);
759 if(spin!=settings.get_bool("XRSSpin")) {
760 // fprintf(stderr,"Excited spin does not match.\n");
761 ok=LOAD_FAIL;
762 }
763 }
764
765 // Check consistency of method
766 if(ok) {
767 std::string method;
768 chkpt.read("XRSMethod",method);
769 if(stricmp(method,settings.get_string("XRSMethod"))!=0) {
770 ok=LOAD_DIFF;
771 // fprintf(stderr,"Calculation methods do not match.\n");
772 }
773 }
774
775 if(ok) {
776 // Everything is OK, finish by loading initial state core orbital
777 if(chkpt.exist("Ccore")) {
778 printf("Loaded initial state orbital from checkpoint.\n");
779 chkpt.read("Ccore",core);
780 } else
781 // Failed to read core orbital
782 ok=LOAD_FAIL;
783
784 } else {
785 // Failed to load or solution was not consistent.
786 sol.Ca=arma::mat();
787 sol.Cb=arma::mat();
788 sol.Ea=arma::vec();
789 sol.Eb=arma::vec();
790 sol.Pa=arma::mat();
791 sol.Pb=arma::mat();
792 }
793
794 return ok;
795 }
796
print_header()797 void print_header() {
798 #ifdef _OPENMP
799 printf("ERKALE - XRS from Hel, OpenMP version, running on %i cores.\n",omp_get_max_threads());
800 #else
801 printf("ERKALE - XRS from Hel, serial version.\n");
802 #endif
803 print_copyright();
804 print_license();
805 #ifdef SVNRELEASE
806 printf("At svn revision %s.\n\n",SVNREVISION);
807 #endif
808 print_hostname();
809 }
810
main_guarded(int argc,char ** argv)811 int main_guarded(int argc, char **argv) {
812 print_header();
813
814 if(argc!=2) {
815 printf("Usage: $ %s runfile\n",argv[0]);
816 return 0;
817 }
818
819 // Initialize libint
820 init_libint_base();
821
822 Timer t;
823 t.print_time();
824
825 // Parse settings
826 settings.add_scf_settings();
827
828 // Change default log file
829 settings.set_string("Logfile","erkale_xrs.log");
830
831 // Use convergence settings similar to StoBe. XCH calculations
832 // are hard to converge to the otherwise default settings.
833 settings.set_double("ConvThr",1e-5);
834 settings.set_double("DFTDelta",100);
835
836 // Add xrs specific settings
837 settings.add_string("LoadChk","Initialize with ground state calculation from file","");
838 settings.add_string("SaveChk","Try initializing with and save results to file","erkale_xrs.chk");
839 settings.add_string("AugChk","Save augmented calculation to file (if applicable)","erkale_xrs_aug.chk");
840
841 settings.add_string("XRSDoubleBasis","The augmentation basis to use for double-basis set calculations","X-AUTO");
842
843 settings.add_bool("XRSSpin","Spin to excite (false for alpha, true for beta)",false);
844 settings.add_string("XRSInitialState","Initial atomic state to excite","1s");
845 settings.add_int("XRSInitialOrbital","Index of orbital in state to excite", 1);
846 settings.add_string("XRSMethod", "Which kind of calculation to perform: TP, XCH or FCH","TP");
847 settings.add_string("XRSAugment","Which atoms to augment with diffuse functions? E.g. 1,3:5,10","");
848 settings.add_string("XRSQval","List or range of Q values to compute","");
849 settings.add_string("XRSQMethod","Method of computing momentum transfer matrix: Local, Fourier or Series","Fourier");
850
851 settings.add_int("XRSNrad","Local: how many point to use in radial integration",200);
852 settings.add_int("XRSLmax","Local: expand orbitals up to Lmax",5);
853 settings.add_int("XRSLquad","Local: perform angular expansion using quadrature of Lquad order",30);
854
855 settings.parse(std::string(argv[1]),true);
856
857 // Redirect output?
858 std::string logfile=settings.get_string("Logfile");
859 if(stricmp(logfile,"stdout")!=0) {
860 // Redirect stdout to file
861 FILE *outstream=freopen(logfile.c_str(),"w",stdout);
862 if(outstream==NULL) {
863 ERROR_INFO();
864 throw std::runtime_error("Unable to redirect output!\n");
865 }
866
867 // Print out the header in the logfile too
868 print_header();
869 fprintf(stderr,"\n");
870 }
871
872 // Get used settings
873 const bool verbose=settings.get_bool("Verbose");
874 const bool spin=settings.get_bool("XRSSpin");
875
876 // Parse method
877 enum xrs_method method=parse_method(settings.get_string("XRSMethod"));
878
879 // Print out settings
880 if(verbose)
881 settings.print();
882
883 // Read in atoms.
884 std::vector<atom_t> atoms;
885 std::string atomfile=settings.get_string("System");
886 if(file_exists(atomfile))
887 atoms=load_xyz(atomfile,!settings.get_bool("InputBohr"));
888 else {
889 // Check if a directory has been set
890 char * libloc=getenv("ERKALE_SYSDIR");
891 if(libloc) {
892 std::string filename=std::string(libloc)+"/"+atomfile;
893 if(file_exists(filename))
894 atoms=load_xyz(filename,!settings.get_bool("InputBohr"));
895 else
896 throw std::runtime_error("Unable to open xyz input file!\n");
897 } else
898 throw std::runtime_error("Unable to open xyz input file!\n");
899 }
900
901 // Get index of excited atom
902 size_t xcatom=get_excited_atom_idx(atoms);
903
904 // Read in basis set
905 BasisSetLibrary baslib;
906 std::string basfile=settings.get_string("Basis");
907 baslib.load_basis(basfile);
908
909 // Construct basis set
910 BasisSet basis;
911 construct_basis(basis,atoms,baslib);
912
913 // Get exchange and correlation functionals and grid settings
914 dft_t dft_init(parse_dft(true));
915 dft_t dft(parse_dft(false));
916
917 // Final convergence settings
918 double convthr(settings.get_double("ConvThr"));
919
920 // Make initialization parameters more relaxed
921 double initfac=settings.get_double("DFTDelta");
922 double init_convthr(convthr*initfac);
923
924 // TP solution
925 uscf_t sol;
926 sol.en.E=0.0;
927
928 // Index of excited orbital
929 size_t xcorb;
930
931 // Number of occupied states
932 int nocca, noccb;
933 get_Nel_alpha_beta(basis.Ztot()-settings.get_int("Charge"),settings.get_int("Multiplicity"),nocca,noccb);
934
935 // Ground state energy
936 energy_t gsen;
937 bool didgs=false;
938
939 // Initial state
940 std::string state=settings.get_string("XRSInitialState");
941 // Initial orbital
942 int iorb=settings.get_int("XRSInitialOrbital")-1;
943
944 // Try to load orbitals and energies
945 arma::vec core;
946 enum loadresult loadok=LOAD_FAIL;
947 if(file_exists(settings.get_string("SaveChk"))) {
948 Checkpoint testload(settings.get_string("SaveChk"),false);
949 loadok=load(basis,testload,sol,core);
950
951 if(loadok==LOAD_SUCC)
952 fprintf(stderr,"Loaded converged calculation from checkpoint file.\n");
953 else if(loadok==LOAD_DIFF) {
954 std::string met;
955 testload.read("XRSMethod",met);
956 fprintf(stderr,"Initializing with previous %s calculation.\n",met.c_str());
957 } else
958 fprintf(stderr,"No converged calculation found in SaveChk.\n");
959
960 // Can we get a ground state energy from the checkpoint file?
961 if(loadok==LOAD_DIFF && testload.exist("E_GS")) {
962 testload.read("E_GS",gsen.E);
963 didgs=true;
964 }
965 }
966
967 // No existing calculation found or system was different => perform calculation
968 if(loadok==LOAD_DIFF || loadok==LOAD_FAIL) {
969
970 // Create checkpoint
971 Checkpoint chkpt(settings.get_string("SaveChk"),true);
972 chkpt.write(basis);
973 chkpt.write("XRSSpin",settings.get_bool("XRSSpin"));
974 chkpt.write("XRSMethod",settings.get_string("XRSMethod"));
975
976 // Index of core orbital
977 int icore=-1;
978
979 // Initialize calculation with ground state if necessary
980 if(loadok==LOAD_FAIL) {
981 if(settings.get_string("LoadChk").size()==0)
982 throw std::runtime_error("Need a ground-state calculation in LoadChk to do calculation!\n");
983
984 printf("Initializing with ground-state calculation from %s.\n",settings.get_string("LoadChk").c_str());
985
986 // Read checkpoint file
987 Checkpoint load(settings.get_string("LoadChk"),false);
988
989 // Restricted calculation?
990 bool restr;
991 load.read("Restricted",restr);
992
993 // Load ground-state energy
994 load.read(gsen);
995 didgs=true;
996
997 // Load basis
998 BasisSet oldbas;
999 load.read(oldbas);
1000
1001 // Check that basis sets are the same
1002 if(!(basis == oldbas))
1003 throw std::runtime_error("Must use the same basis for all phases of XRS calculation!\n");
1004
1005 // Read orbitals
1006 if(restr) {
1007 load.read("C",sol.Ca);
1008 load.read("E",sol.Ea);
1009 load.read("H",sol.Ha);
1010 load.read("P",sol.P);
1011 sol.Eb=sol.Ea;
1012 sol.Cb=sol.Ca;
1013 sol.Hb=sol.Ha;
1014 sol.Pa=sol.Pb=sol.P/2.0;
1015 } else {
1016 // Load energies and orbitals
1017 load.read("Ca",sol.Ca);
1018 load.read("Ea",sol.Ea);
1019 load.read("Ha",sol.Ha);
1020 load.read("Pa",sol.Pa);
1021 load.read("Cb",sol.Cb);
1022 load.read("Eb",sol.Eb);
1023 load.read("Hb",sol.Hb);
1024 load.read("Pb",sol.Pb);
1025 load.read("P",sol.P);
1026 }
1027
1028 // Determine orbitals, update energies and densities
1029 if(spin) {
1030 icore=localize(basis,noccb,xcatom,sol.Cb,state,iorb);
1031 sol.Eb=arma::diagvec(arma::trans(sol.Cb)*sol.Hb*sol.Cb);
1032
1033 std::vector<double> occb;
1034 if(method==XCH)
1035 occb=xch_occ(icore,noccb);
1036 else if(method==FCH)
1037 occb=fch_occ(icore,noccb);
1038 else if(method==TP)
1039 occb=tp_occ(icore,noccb);
1040 sol.Pb=form_density(sol.Cb,occb);
1041
1042 } else {
1043 icore=localize(basis,nocca,xcatom,sol.Ca,state,iorb);
1044 sol.Ea=arma::diagvec(arma::trans(sol.Ca)*sol.Ha*sol.Ca);
1045
1046 std::vector<double> occa;
1047 if(method==XCH)
1048 occa=xch_occ(icore,nocca);
1049 else if(method==FCH)
1050 occa=fch_occ(icore,nocca);
1051 else if(method==TP)
1052 occa=tp_occ(icore,nocca);
1053 sol.Pa=form_density(sol.Ca,occa);
1054 }
1055 sol.P=sol.Pa+sol.Pb;
1056 }
1057
1058 // Proceed with TP calculation. Initialize solver
1059 XRSSCF solver(basis,chkpt,spin);
1060
1061 // Set initial state core orbital
1062 if(loadok==LOAD_FAIL) { // LOAD_FAIL
1063 if(spin)
1064 solver.set_core(sol.Cb.col(icore));
1065 else
1066 solver.set_core(sol.Ca.col(icore));
1067 } else
1068 // LOAD_DIFF
1069 solver.set_core(core);
1070
1071 // Write number of electrons to file
1072 chkpt.write("Nel",nocca+noccb);
1073 chkpt.write("Nel-a",nocca);
1074 chkpt.write("Nel-b",noccb);
1075
1076 // Store core orbital
1077 chkpt.write("Ccore",solver.get_core());
1078
1079 // Do we have the ground state energy?
1080 if(didgs)
1081 chkpt.write("E_GS",gsen.E);
1082
1083 // Do calculation
1084 if(method==FCH || method==XCH) {
1085 if(dft.adaptive)
1086 solver.full_hole(sol,init_convthr,dft_init,method==XCH);
1087 xcorb=solver.full_hole(sol,convthr,dft,method==XCH);
1088
1089 // Get excited state energy
1090 energy_t excen;
1091 chkpt.read(excen);
1092
1093 // Do we have a ground state energy?
1094 if(didgs) {
1095 if(method==XCH) {
1096 printf("\nAbsolute energy correction: first peak should be at %.2f eV.\n",(excen.E-gsen.E)*HARTREEINEV);
1097 fprintf(stderr,"\nAbsolute energy correction: first peak should be at %.2f eV.\n",(excen.E-gsen.E)*HARTREEINEV);
1098 } else {
1099 printf("Vertical photoionization energy is %.2f eV.\n",(excen.E-gsen.E)*HARTREEINEV);
1100 fprintf(stderr,"Vertical photoionization energy is %.2f eV.\n",(excen.E-gsen.E)*HARTREEINEV);
1101 }
1102 } else {
1103 printf("Cannot estimate excitation energy without a ground state calculation.\n");
1104 fprintf(stderr,"Cannot estimate excitation energy without a ground state calculation.\n");
1105 }
1106 } else {
1107 if(dft.adaptive)
1108 solver.half_hole(sol,init_convthr,dft_init);
1109 xcorb=solver.half_hole(sol,convthr,dft);
1110 }
1111 printf("\n\n");
1112
1113 } else {
1114 // Loaded calculation, just find excited state
1115
1116 if(spin)
1117 xcorb=find_excited_orb(basis,core,sol.Cb,noccb);
1118 else
1119 xcorb=find_excited_orb(basis,core,sol.Ca,nocca);
1120 }
1121
1122 // Augment the solutions if necessary
1123 BasisSet augbas;
1124 arma::mat C_aug;
1125 arma::vec E_aug;
1126
1127 if(stricmp(settings.get_string("XRSAugment"),"")!=0)
1128 augmented_solution(basis,sol,nocca,noccb,dft,augbas,C_aug,E_aug,spin,method);
1129 else {
1130 // No augmentation necessary, just copy the solutions
1131 augbas=basis;
1132 if(spin) {
1133 C_aug=sol.Cb;
1134 E_aug=sol.Eb;
1135 } else {
1136 C_aug=sol.Ca;
1137 E_aug=sol.Ea;
1138 }
1139 }
1140
1141 // Number of occupied states
1142 size_t nocc;
1143 if(spin)
1144 nocc=noccb;
1145 else
1146 nocc=nocca;
1147
1148
1149 // Compute dipole transitions
1150 std::vector<spectrum_t> sp=compute_transitions(augbas,C_aug,E_aug,xcatom,xcorb,nocc);
1151 // Save spectrum
1152 save_spectrum(sp);
1153 save_xas_spectrum(sp);
1154
1155 // Get values of q to compute for
1156 std::vector<double> qvals=parse_range_double(settings.get_string("XRSQval"));
1157
1158 // The q-dependent spectra
1159 std::vector< std::vector<spectrum_t> > qsp;
1160 // The filename
1161 std::string spname;
1162
1163 if(qvals.size()) {
1164 // Series method
1165 if(stricmp(settings.get_string("XRSQMethod"),"Series")==0) {
1166 qsp=compute_qdep_transitions_series(augbas,C_aug,E_aug,xcorb,nocc,qvals);
1167 spname="trans_ser";
1168 }
1169
1170 // Fourier method
1171 if(stricmp(settings.get_string("XRSQMethod"),"Fourier")==0) {
1172 qsp=compute_qdep_transitions_fourier(augbas,C_aug,E_aug,xcorb,nocc,qvals);
1173 spname="trans_four";
1174 }
1175
1176 // Local method (Sakko et al)
1177 if(stricmp(settings.get_string("XRSQMethod"),"Local")==0) {
1178 qsp=compute_qdep_transitions_local(augbas,C_aug,E_aug,xcatom,xcorb,nocc,qvals);
1179 spname="trans_loc";
1180 }
1181
1182 // Save transitions
1183 for(size_t i=0;i<qvals.size();i++) {
1184 char fname[80];
1185 sprintf(fname,"%s-%.2f.dat",spname.c_str(),qvals[i]);
1186 save_spectrum(qsp[i],fname);
1187 }
1188 }
1189
1190 if(verbose) {
1191 printf("\nRunning program took %s.\n",t.elapsed().c_str());
1192 t.print_time();
1193 }
1194
1195 return 0;
1196 }
1197
main(int argc,char ** argv)1198 int main(int argc, char **argv) {
1199 try {
1200 return main_guarded(argc, argv);
1201 } catch (const std::exception &e) {
1202 std::cerr << "error: " << e.what() << std::endl;
1203 return 1;
1204 }
1205 }
1206