ERKALE
ERKALE - DFT from Hel
 All Classes Functions Variables Friends Pages
co-opt.h
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-2014
9  * Copyright (c) 2010-2014, 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 
46 #ifndef ERKALE_CO_OPT
47 #define ERKALE_CO_OPT
48 
49 #include <cstdio>
50 #include <unistd.h>
51 #include <iostream>
52 #include "../completeness/optimize_completeness.h"
53 #include "../completeness/completeness_profile.h"
54 #include "../external/fchkpt_tools.h"
55 #include "../basislibrary.h"
56 #include "../global.h"
57 #include "../timer.h"
58 #include "../stringutil.h"
59 #include "../unitary.h"
60 #include "../elements.h"
61 #include "../mathf.h"
62 
63 #ifdef _OPENMP
64 #include <omp.h>
65 #endif
66 
68 #define NFMAX 70
69 
71 typedef struct {
73  double start;
75  double end;
77  double tol;
78 
80  arma::vec exps;
81 } coprof_t;
82 
84 int maxam(const std::vector<coprof_t> & cpl);
86 int get_nfuncs(const std::vector<coprof_t> & cpl);
88 void print_limits(const std::vector<coprof_t> & cpl, const char *msg=NULL);
89 
91 std::vector<coprof_t> augment_basis(const std::vector<coprof_t> & cplo, int ndiffuse, int ntight);
92 
93 // Construct elemental library
94 ElementBasisSet get_element_library(const std::string & el, const std::vector<coprof_t> & cpl);
95 // Construct universal library
96 BasisSetLibrary get_library(const std::vector<coprof_t> & cpl);
97 
99 void print_scheme(const BasisSetLibrary & baslib, int len=0);
100 
102 typedef struct {
103  // Completeness
104  double tol;
105  // Width
106  double w;
107  // Exponents
108  arma::vec exps;
109 } co_exps_t;
110 
112 arma::vec maxwidth_exps_table(int am, double tol, size_t nexp, double & width, int n_cpl=1);
114 arma::vec span_width(int am, double tol, double & width, int nx=0, int n_cpl=1);
115 
117 template<typename ValueType>
119  private:
121  int n_tau;
123  int n_full;
124 
126  size_t politer;
127 
129  virtual double compute_energy(const std::vector<coprof_t> & cpl)=0;
130 
132  virtual BasisSetLibrary form_basis(const std::vector<coprof_t> & cpl, std::vector<size_t> contract, bool Porth=true) {
134  (void) contract;
135  (void) Porth;
136  return form_basis(cpl);
137  }
138 
140  virtual void update_reference(const std::vector<coprof_t> & cpl) {
142  (void) cpl;
143  };
144 
146  void mog_decompose(const std::vector<int> & tram, const arma::vec & mogs, std::vector<arma::uvec> & amidx, std::vector<arma::vec> & ammog) const {
147  if(tram.size() != mogs.n_elem) {
148  std::ostringstream oss;
149  oss << "Error: trial am size = " << tram.size() << ", but mog size = " << mogs.n_elem << "!\n";
150  throw std::runtime_error(oss.str());
151  }
152 
153  // Sanity check
154  if(!tram.size()) {
155  amidx.clear();
156  ammog.clear();
157  return;
158  }
159 
160  // Maximum angular momentum
161  int maxam=arma::conv_to<arma::ivec>::from(tram).max();
162 
163  amidx.clear();
164  amidx.resize(maxam+1);
165  ammog.clear();
166  ammog.resize(maxam+1);
167 
168  // Mog of ams
169  std::vector< std::vector<double> > amm(maxam+1);
170  // Indices
171  std::vector< std::vector< arma::uword> > ami(maxam+1);
172  for(arma::uword i=0;i<tram.size();i++) {
173  ami[tram[i]].push_back(i);
174  amm[tram[i]].push_back(mogs(i));
175  }
176 
177  // Store mogs
178  for(int am=0;am<=maxam;am++) {
179  if(ami[am].size()) {
180  amidx[am]=arma::conv_to<arma::uvec>::from(ami[am]);
181  ammog[am]=arma::conv_to<arma::vec>::from(amm[am]);
182  } else {
183  amidx[am].zeros(0);
184  ammog[am].zeros(0);
185  }
186  }
187  }
188 
189  public:
193  n_tau=1;
195  n_full=4;
197  politer=0;
198  }
201  }
202 
204  int get_ntau() const {
205  return n_tau;
206  }
207 
209  void set_ntau(int n) {
210  n_tau=n;
211  }
212 
214  int get_nfull() const {
215  return n_full;
216  }
217 
219  void set_nfull(int n) {
220  n_full=n;
221  }
222 
224  arma::vec optimize_completeness(int am, double min, double max, int Nf, double & mog) {
225  return ::optimize_completeness(am, min, max, Nf, n_tau, false, &mog, n_full);
226  }
227 
229  arma::vec maxwidth_exps_table(int am, double tol, size_t nexp, double & width) {
230  // Optimized exponents
231  static std::vector<co_exps_t> opt(max_am+1);
232 
233  // Check if we already have the exponents in store
234  if(opt[am].tol!=tol || opt[am].exps.size()!=nexp) {
235  opt[am].tol=tol;
236  opt[am].exps=maxwidth_exps(am,tol,nexp,opt[am].w,n_tau,n_full);
237  }
238 
239  width=opt[am].w;
240 
241  if(opt[am].exps.size() != nexp) {
242  std::ostringstream oss;
243  oss << "Required " << nexp << " completeness-optimized primitives but got " << opt[am].exps.size() << "!\n";
244  throw std::runtime_error(oss.str());
245  }
246 
247  return opt[am].exps;
248  }
249 
251  arma::vec span_width(int am, double tol, double & width, int nx) {
252  // Check starting point
253  if(nx<0)
254  nx=0;
255 
256  // Determine necessary amount of exponents
257  arma::vec exps;
258  double w=0.0;
259 
260  for(nx++;nx<=NFMAX;nx++) {
261  exps=maxwidth_exps_table(am,tol,nx,w);
262  if(w>=width)
263  break;
264  }
265 
266  // Store real width
267  width=w;
268 
269  // Return exponents
270  return exps;
271  }
272 
273 
275  virtual BasisSetLibrary form_basis(const std::vector<coprof_t> & cpl) {
276  BasisSetLibrary baslib;
277  for(int Z=1;Z<=maxZ;Z++)
278  baslib.add_element(get_element_library(element_symbols[Z],cpl));
279 
280  return baslib;
281  }
282 
284  void save_limits(const std::vector<coprof_t> & cpl, const std::string & fname) const {
285  FILE *out=fopen(fname.c_str(),"w");
286  if(!out)
287  throw std::runtime_error("Error opening completeness range output file.\n");
288 
289  fprintf(out,"%i\n",maxam(cpl));
290  for(int am=0;am<=maxam(cpl);am++)
291  fprintf(out,"%i % .16e % .16e %.16e %i\n",am,cpl[am].start,cpl[am].end,cpl[am].tol,(int) cpl[am].exps.size());
292  fclose(out);
293  }
294 
296  std::vector<coprof_t> load_limits(const std::string & fname) {
297  FILE *in=fopen(fname.c_str(),"r");
298  if(!in)
299  throw std::runtime_error("Error opening completeness range file.\n");
300 
301  int max;
302  if(fscanf(in,"%i",&max)!=1)
303  throw std::runtime_error("Could not read maximum angular momentum from file.\n");
304  if(max<0 || max>max_am) {
305  std::ostringstream oss;
306  oss << "Error - read in maximum angular momentum " << max << "!\n";
307  throw std::runtime_error(oss.str());
308  }
309 
310  // Allocate memory
311  std::vector<coprof_t> cpl(max+1);
312 
313  // Read in ranges and make exponents
314  for(int am=0;am<=max;am++) {
315  // Supposed angular momentum, and amount of exponents
316  int amt, nexp;
317  // Tolerance
318  double ctol;
319 
320  // Read range
321  if(fscanf(in,"%i %lf %lf %lf %i",&amt,&cpl[am].start,&cpl[am].end,&ctol,&nexp)!=5) {
322  std::ostringstream oss;
323  oss << "Error reading completeness range for " << shell_types[am] << " shell!\n";
324  throw std::runtime_error(oss.str());
325  }
326  // Check am is OK
327  if(am!=amt) {
328  std::ostringstream oss;
329  oss << "Read in am " << amt << " does not match supposed am " << am << "!\n";
330  throw std::runtime_error(oss.str());
331  }
332 
333  // Get exponents
334  arma::vec exps=optimize_completeness(am,0.0,cpl[am].end-cpl[am].start,nexp,cpl[am].tol);
335 
336  // Check that tolerances agree
337  if(fabs(ctol-cpl[am].tol)>1e-3*cpl[am].tol) {
338  std::ostringstream oss;
339  oss << "Obtained tolerance " << cpl[am].tol << " for " << shell_types[am] << " shell does not match supposed tolerance " << ctol << "!\n";
340  // throw std::runtime_error(oss.str());
341  printf("Warning: %s",oss.str().c_str());
342  fflush(stdout);
343 
344  // Find correct width
345  double width=cpl[am].end-cpl[am].start;
346 
347  // Get exponents
348  double w;
349  exps=maxwidth_exps_table(am,cpl[am].tol,nexp,w);
350 
351  // Store exponents
352  double dw=w-width;
353  cpl[am].start-=dw/2.0;
354  cpl[am].end+=dw/2.0;
355  }
356 
357  cpl[am].exps=move_exps(exps,cpl[am].start);
358  }
359  fclose(in);
360 
361  return cpl;
362  }
363 
365  std::vector<coprof_t> load_limits(const BasisSetLibrary & baslib, double tol, int maxam) {
366  // Threshold for completeness in initial set
367  double cplthr=1.0-tol;
368 
369  // Returned array
370  std::vector<coprof_t> cpl(max_am+1);
371  for(size_t i=0;i<cpl.size();i++) {
372  cpl[i].start=0.0;
373  cpl[i].end=0.0;
374  cpl[i].tol=tol;
375  }
376 
377  // Get elements
378  std::vector<ElementBasisSet> els=baslib.get_elements();
379 
380  // Loop over elements.
381  for(size_t iel=0;iel<els.size();iel++) {
382  // Get elemental basis
383  ElementBasisSet elbas=els[iel];
384  // and decontract it
385  elbas.decontract();
386 
387  // Compute completeness profile
388  compprof_t prof=compute_completeness(elbas);
389 
390  // Loop over angular momentum
391  for(size_t am=0;am<prof.shells.size();am++) {
392  // Determine lower limit
393  double low=DBL_MAX;
394  for(size_t ia=0;ia<prof.lga.size();ia++)
395  if(prof.shells[am].Y[ia]>=cplthr) {
396  low=prof.lga[ia];
397  break;
398  }
399 
400  // Determine upper limit
401  double high=-DBL_MAX;
402  for(size_t ia=prof.lga.size()-1;ia<prof.lga.size();ia--)
403  if(prof.shells[am].Y[ia]>=cplthr) {
404  high=prof.lga[ia];
405  break;
406  }
407 
408  if(cpl[am].start>low)
409  cpl[am].start=low;
410  if(cpl[am].end<high)
411  cpl[am].end=high;
412  }
413  }
414 
415  // Clear out extra shells
416  for(int am=maxam;am<max_am;am++) {
417  cpl[am].start=0.0;
418  cpl[am].end=0.0;
419  cpl[am].exps.clear();
420  }
421  // Dummy check
422  for(int am=0;am<std::min(maxam,max_am);am++)
423  if(cpl[am].start==DBL_MAX && cpl[am].end==-DBL_MAX) {
424  cpl[am].start=0.0;
425  cpl[am].end=0.0;
426  }
427 
428  printf("Initial composition:\n");
429 
430  // Form exponents.
431  for(int am=0;am<max_am;am++) {
432  if(cpl[am].start == 0.0 && cpl[am].end == 0.0)
433  continue;
434 
435  // Compute width
436  double width=cpl[am].end-cpl[am].start;
437 
438  // Determine amount of exponents necessary to obtain width
439  std::vector<double> widths;
440  std::vector< std::vector<double> > explist;
441 
442  widths.push_back(0.0);
443  explist.resize(1);
444 
445  for(int nf=1;nf<=NFMAX;nf++) {
446  double w;
447  std::vector<double> exps=maxwidth_exps_table(am,cpl[am].tol,nf,w);
448 
449  widths.push_back(w);
450  explist.push_back(exps);
451 
452  // Have we reached the necessary width?
453  if(w>width)
454  break;
455  }
456 
457  // Which amount of functions is closest?
458  size_t nf=0;
459  double mind=DBL_MAX;
460  for(size_t i=0;i<widths.size();i++)
461  if(fabs(widths[i]-width)<mind) {
462  nf=i;
463  mind=fabs(widths[i]-width);
464  }
465 
466  // Yes, we have. Adjust starting and ending points
467  cpl[am].start-=(widths[nf]-width)/2.0;
468  cpl[am].end+=(widths[nf]-width)/2.0;
469  // and store exponents, moving starting point to the correct location
470  cpl[am].exps=move_exps(explist[nf],cpl[am].start);
471 
472  printf("%c % .2f % .2f %2i\n",shell_types[am],cpl[am].start,cpl[am].end,(int) cpl[am].exps.size());
473  fflush(stdout);
474  }
475 
476  printf("\n");
477 
478  return cpl;
479  }
480 
481 
486  std::vector<coprof_t> initial_profile(const std::vector<int> & nshells, int nexp, double cotol=1e-4) {
487  std::vector<size_t> nfuncs(nshells.begin(),nshells.end());
488  for(size_t i=0;i<nfuncs.size();i++)
489  nfuncs[i]*=nexp;
490  return initial_profile(nfuncs,cotol);
491  }
492 
497  std::vector<coprof_t> initial_profile(const std::vector<size_t> & nfuncs, double cotol=1e-4) {
498  std::vector<double> offsets(nfuncs.size(),0.0);
499  return initial_profile(nfuncs,offsets,cotol);
500  }
501 
507  std::vector<coprof_t> initial_profile(const std::vector<size_t> & nfuncs, std::vector<double> offsets, double cotol=1e-4) {
508  Timer t;
509 
510  // Initialize profile
511  std::vector<coprof_t> cpl(nfuncs.size());
512  for(size_t i=0;i<cpl.size();i++) {
513  cpl[i].start=0.0;
514  cpl[i].end=0.0;
515  cpl[i].tol=cotol;
516  }
517 
518  // Exponents
519  std::vector<arma::vec> exps(cpl.size());
520  arma::vec widths(cpl.size());
521 
522  // Generate initial profile
523  printf("\nStarting composition:\n");
524  for(size_t am=0;am<cpl.size();am++) {
525  exps[am]=maxwidth_exps_table(am,cpl[am].tol,nfuncs[am],widths[am]);
526  cpl[am].start=offsets[am];
527  cpl[am].end=cpl[am].start+widths[am];
528  cpl[am].exps=exps[am];
529 
530  printf("%c % .3f % .3f %e %2i\n",shell_types[am],cpl[am].start,cpl[am].end,cpl[am].tol,(int) cpl[am].exps.size());
531  fflush(stdout);
532  }
533 
534  // Widths normalized by amount of exponents
535  arma::vec relw(widths);
536  for(size_t am=0;am<relw.size();am++)
537  relw(am)/=exps[am].size();
538 
539  // Convergence loop
540  const double h0=1e-3; // Initial step size in relative width
541  size_t iiter=0;
542 
543  // Old gradient
544  arma::vec gold;
545 
546  while(true) {
547  Timer tstep;
548 
549  // Conjugate gradient step?
550  bool cg=iiter%cpl.size()!=0;
551  if(cg)
552  printf("\nConjugate gradient iteration %i\n",(int) iiter+1);
553  else
554  printf("\nSteepest descent iteration %i\n",(int) iiter+1);
555  fflush(stdout);
556 
557  // Compute finite difference gradient
558  arma::vec g(cpl.size());
559  for(size_t am=0;am<cpl.size();am++) {
560  // Step size to use for shell
561  double h=h0*relw(am);
562 
563  std::vector<coprof_t> lcpl(cpl);
564  lcpl[am].start-=h;
565  lcpl[am].end-=h;
566  lcpl[am].exps=move_exps(exps[am],lcpl[am].start);
567 
568  std::vector<coprof_t> rcpl(cpl);
569  rcpl[am].start+=h;
570  rcpl[am].end+=h;
571  rcpl[am].exps=move_exps(exps[am],rcpl[am].start);
572 
573  // Energies
574  double El=compute_energy(lcpl);
575  double Er=compute_energy(rcpl);
576 
577  bool Elok=(std::isfinite(El) && fabs(El)!=DBL_MAX);
578  bool Erok=(std::isfinite(Er) && fabs(Er)!=DBL_MAX);
579  // Finite difference gradient. Absorb relative width here
580  if(Elok && Erok)
581  g(am)=(Er-El)/(2*h) * relw(am);
582  else
583  // Something wrong with gradient.
584  g(am)=0.0;
585  printf("\t-grad(%c) = % e\n", shell_types[am], -g(am));
586  fflush(stdout);
587  }
588 
589  {
590  // Print bar
591  size_t Nm=26;
592  char minus[Nm];
593  memset(minus,'-',Nm);
594  minus[Nm-1]='\0';
595  printf("\t%s\n",minus);
596  }
597  printf("\tgrad norm = % e\n",arma::norm(g,2));
598  fflush(stdout);
599 
600  // Normalize gradient
601  g/=arma::norm(g,2);
602 
603  // Orthogonalize wrt old gradient?
604  if(cg) {
605  g-=arma::dot(g,gold)*gold;
606  g/=arma::norm(g,2);
607  }
608 
609  // Store old gradient
610  gold=g;
611 
612  // Line search
613  const int Ntr=20; // This is ridiculously large to prevent runover
614  arma::vec h(Ntr);
615  arma::vec E(Ntr);
616  // Initial value might be positive, so
617  // initialize to largest possible value
618  E.ones();
619  E*=DBL_MAX;
620 
621  // Step sizes
622  h(0)=0.0;
623  double hs=10*h0;
624  for(int i=1;i<Ntr;i++)
625  h(i)=hs*std::pow(2.0,i-1);
626 
627  printf("\n\tLine search:\n");
628  printf("\t%5s %12s %13s %13s\n","trial","step","E","dE");
629  fflush(stdout);
630 
631  int itr;
632  for(itr=0;itr<Ntr;itr++) {
633  // Profile
634  std::vector<coprof_t> trcpl(cpl);
635  for(size_t am=0;am<g.n_elem;am++) {
636  trcpl[am].start-=h(itr)*g(am);
637  trcpl[am].end-=h(itr)*g(am);
638  trcpl[am].exps=move_exps(exps[am],trcpl[am].start);
639  }
640 
641  // Energy
642  E(itr)=compute_energy(trcpl);
643 
644  if(itr==0)
645  printf("\t%2i/%-2i %e % e\n",(int) (itr+1), (int) Ntr, h(itr), E(itr));
646  else
647  printf("\t%2i/%-2i %e % e % e\n",(int) (itr+1), (int) Ntr, h(itr), E(itr), E(itr)-E(0));
648  fflush(stdout);
649 
650  // Check if energy increased
651  if(itr>0 && E(itr)>E(itr-1))
652  break;
653  }
654 
655  // Find minimal energy
656  arma::uword ind;
657  E.min(ind);
658  double step=h(ind);
659 
660  if(ind!=0) {
661  printf("\tMinimum with step %e, decrease by % e.\n",h(ind),E(ind)-E(0));
662  fflush(stdout);
663 
664  if(ind>=1) {
665  // Perform cubic interpolation.
666 
667  arma::vec x=h.subvec(ind-1,ind+1);
668  arma::vec y=E.subvec(ind-1,ind+1);
669 
670  // Fit polynomial
671  arma::vec c=fit_polynomial(x,y);
672  // Derivative coefficients
673  arma::vec dc=derivative_coefficients(c);
674  // Get smallest positive root
675  double st=smallest_positive(solve_roots(dc));
676 
677  printf("\tInterpolation gives step %e, ",st);
678  fflush(stdout);
679 
680  // Trial profile
681  std::vector<coprof_t> trcpl(cpl);
682  for(size_t am=0;am<g.n_elem;am++) {
683  trcpl[am].start-=st*g(am);
684  trcpl[am].end-=st*g(am);
685  trcpl[am].exps=move_exps(exps[am],trcpl[am].start);
686  }
687 
688  // Compute new trial value
689  double Etr=compute_energy(trcpl);
690 
691  if(Etr<E.min()) {
692  step=st;
693  printf("further decrease by % e, accepted.\n",Etr-E.min());
694  } else
695  printf("further increase by % e, rejected.\n",Etr-E.min());
696  fflush(stdout);
697  }
698 
699  // Update profile limits
700  for(size_t am=0;am<g.n_elem;am++) {
701  cpl[am].start-=step*g(am);
702  cpl[am].end-=step*g(am);
703  cpl[am].exps=move_exps(exps[am],cpl[am].start);
704  }
705 
706  printf("Iteration took %s.\n",tstep.elapsed().c_str());
707  print_limits(cpl,"\nCurrent limits");
708  fflush(stdout);
709 
710  } else
711  break;
712 
713  iiter++;
714  }
715 
716  print_limits(cpl,"Initial composition");
717  printf("Generation of initial profile took %s.\n",t.elapsed().c_str());
718  fflush(stdout);
719 
720  return cpl;
721  }
722 
724  virtual std::vector<ValueType> compute_values(const std::vector<BasisSetLibrary> & baslib)=0;
726  std::vector<ValueType> compute_values(const std::vector< std::vector<coprof_t> > & cpl) {
727  std::vector<BasisSetLibrary> baslibs(cpl.size());
728  for(size_t i=0;i<cpl.size();i++)
729  baslibs[i]=form_basis(cpl[i]);
730  return compute_values(baslibs);
731  }
732 
733  // Compute distance of two values. If calculation of val failed, set mog to failval.
734  virtual double compute_mog(const ValueType & val, const ValueType & ref, double failval)=0;
735 
749  double check_polarization(std::vector<coprof_t> & cpl, ValueType & curval, int am_max, double minpol, double maxpol, double dpol, bool polinterp, double cotol, int nx, bool forcesub=false) {
750 
751  printf("\n\n%s\n",print_bar("POLARIZATION SHELLS").c_str());
752  fflush(stdout);
753 
754  // Check necessity of additional polarization shell.
755  Timer t;
756  const int max=maxam(cpl);
757  if(max>=am_max) {
758  printf("Won't add another polarization shell, MaxAM=%i reached.\n",am_max);
759  fflush(stdout);
760  return 0.0;
761  }
762 
763  Timer tpol;
764 
765  // Width of completeness
766  double pw;
767 
768  // Exponent to add: l=max+1, tolerance is cotol, 1 function.
769  const int addam=max+1;
770  arma::vec pexp=maxwidth_exps_table(addam,cotol,nx,pw);
771 
772  if(dpol<=0) {
773  std::ostringstream oss;
774  oss << "Invalid value " << dpol << " for DPol.\n";
775  throw std::runtime_error(oss.str());
776  }
777 
778  // Compute spacing with 5 exponents; this is already pretty much saturated.
779  double sp;
780  {
781  int nexp=5;
782  double wnext;
783  maxwidth_exps_table(addam,cotol,nexp+1,wnext);
784  double wcur;
785  maxwidth_exps_table(addam,cotol,nexp,wcur);
786  sp=(wnext-wcur)/dpol;
787  }
788  // Wanted width is
789  double ww=maxpol-minpol;
790  // Amount of points to use
791  size_t Np=1+round(ww/sp);
792  // Actual width is
793  double width=sp*Np;
794 
795  // so actual minimum and maximum are
796  double minp=minpol;
797  double maxp=maxpol+(width-ww);
798 
799  // Starting points
800  arma::vec startp=arma::linspace(minp,maxp,Np);
801 
802  printf("Scanning for optimal placement of %i exponents on %c shell.\n",(int) pexp.size(), shell_types[addam]);
803  printf("Using %.3f points per exponent interval, so spacing between points is %.5f.\n\n",dpol,sp);
804  fflush(stdout);
805 
806  // Allocate sufficient memory
807  if(cpl.size() <= (size_t) addam) {
808  cpl.resize(addam+1);
809  // Set tolerance
810  cpl[addam].tol=cotol;
811  }
812 
813  // Form trials
814  std::vector< std::vector<coprof_t> > cpls(startp.n_elem,cpl);
815  for(size_t iexp=0;iexp<startp.n_elem;iexp++) {
816  // Completeness profile
817  cpls[iexp][addam].start=startp(iexp)-pw/2.0;
818  cpls[iexp][addam].end=startp(iexp)+pw/2.0;
819  cpls[iexp][addam].exps=move_exps(pexp,cpls[iexp][addam].start);
820  }
821 
822  // Calculate trial values
823  std::vector<ValueType> vals=compute_values(cpls);
824  if(vals.size()!=cpls.size()) {
825  std::ostringstream oss;
826  oss << "Error - requested computation of " << cpls.size() << " but got only " << vals.size() << "!\n";
827  throw std::runtime_error(oss.str());
828  }
829 
830  // Resulting mogs
831  printf("%11s %8s %12s\n","trial ","exponent","mog");
832  fflush(stdout);
833  arma::vec mogs(startp.n_elem);
834  for(size_t i=0;i<vals.size();i++) {
835  // Set failed values to zero
836  mogs(i)=compute_mog(vals[i],curval,0.0);
837  printf("%5i/%-5i % 7.5f %e\n",(int) i+1, (int) startp.n_elem, startp(i), mogs(i));
838  }
839 
840  {
841  // Save out the results of the mog scan.
842  // Increment counter.
843  politer++;
844 
845  // Filename to use
846  std::string fname;
847  {
848  std::ostringstream oss;
849  oss << "polmog_" << politer << "_" << shell_types[addam] << ".dat";
850  fname=oss.str();
851  }
852 
853  arma::mat savemog(mogs.n_rows,2);
854  savemog.col(0)=startp;
855  savemog.col(1)=mogs;
856  savemog.save(fname,arma::raw_ascii);
857 
858  // Save limits
859  for(int am=0;am<=maxam(cpl);am++) {
860  // Form limits
861  arma::mat lim(2,3);
862 
863  // y from 0 to infinity
864  lim(0,0)=0.0;
865  lim(1,0)=DBL_MAX;
866 
867  // Starting point
868  lim(0,1)=cpl[am].start;
869  lim(1,1)=cpl[am].start;
870 
871  // Ending point
872  lim(0,2)=cpl[am].end;
873  lim(1,2)=cpl[am].end;
874 
875  // Output
876  std::ostringstream oss;
877  oss << "limits_" << politer << "_" << shell_types[am] << ".dat";
878  lim.save(oss.str(),arma::raw_ascii);
879  }
880 
881  // Save basis set
882  {
883  std::ostringstream oss;
884  oss << "polbas_" << politer << ".gbs";
885  form_basis(cpl).save_gaussian94(oss.str());
886  }
887  }
888 
889  // Find maximum mog
890  arma::uword imax;
891  double maxmog=mogs.max(imax);
892 
893  printf("\n");
894  printf("%11s % 7.5f %e\n","max mog",startp(imax),maxmog);
895  fflush(stdout);
896 
897  if(polinterp) {
898  // Interpolate mogs to find more accurate maximum. Spacing to use
899  double spacing=std::min(1e-6, sp/10.0);
900  arma::vec finergrid=arma::linspace(minp,maxp,1+ceil((maxp-minp)/spacing));
901  arma::vec moginterp=spline_interpolation(startp,mogs,finergrid);
902 
903  arma::uword iintmean;
904  double intmean;
905  intmean=moginterp.max(iintmean);
906 
907  printf("\n");
908  printf("%11s % 7.5f %e\n","interp mog",finergrid(iintmean),intmean);
909  fflush(stdout);
910 
911  if(intmean>maxmog) {
912  // Trial profile
913  std::vector< std::vector<coprof_t> > intcpl(1,cpl);
914  intcpl[0][addam].start=finergrid(iintmean)-pw/2.0;
915  intcpl[0][addam].end=finergrid(iintmean)+pw/2.0;
916  intcpl[0][addam].exps=move_exps(pexp,intcpl[0][addam].start);
917 
918  // Calculate real mog
919  double intmog;
920  ValueType intval;
921  try {
922  std::vector<ValueType> intvals(compute_values(intcpl));
923  if(intvals.size()!=1) {
924  std::ostringstream oss;
925  oss << "Error - requested computation of " << 1 << " values but got only " << intvals.size() << "!\n";
926  throw std::runtime_error(oss.str());
927  }
928 
929  // Compute value
930  intval=intvals[0];
931  // and mog
932  intmog=compute_mog(intval,curval,0.0);
933  } catch(...) {
934  // Interpolation calculation failed.
935  intmog=0.0;
936  }
937 
938  printf("%11s %8s %e\n","real value","",intmog);
939 
940  // Did interpolation yield a better value?
941  if(intmog > maxmog) {
942  printf("Interpolation accepted.\n");
943 
944  // Switch values.
945  cpl=intcpl[0];
946  curval=intval;
947  } else {
948  printf("Interpolation rejected.\n");
949 
950  // Update values
951  cpl=cpls[imax];
952  curval=vals[imax];
953  }
954  fflush(stdout);
955 
956  } else {
957  printf("Interpolation wouldn't result in a better value.\n");
958  fflush(stdout);
959 
960  // Update values
961  cpl=cpls[imax];
962  curval=vals[imax];
963  }
964  } else {
965  // Update values
966  cpl=cpls[imax];
967  curval=vals[imax];
968  }
969 
970  printf("\nOptimal %c shell is: % .3f ... % .3f (%i funcs) with tolerance %e, mog = %e. (%s)\n",shell_types[addam],cpl[addam].start,cpl[addam].end,(int) cpl[addam].exps.size(),cpl[addam].tol,maxmog,tpol.elapsed().c_str());
971  print_value(curval,"Polarized value");
972  fflush(stdout);
973 
974  if(forcesub)
975  comply_subset(cpl,curval);
976 
977  return maxmog;
978  }
979 
981  bool comply_subset(std::vector<coprof_t> & cpl, ValueType & curval) {
982 
983  // Need update?
984  bool upd=false;
985  // Message printed?
986  bool printed=false;
987 
988  for(int am=maxam(cpl);am>=0;am--) {
989  // Move ends by
990  double smove=0.0;
991  double emove=0.0;
992  coprof_t Yam(cpl[am]);
993 
994  // Figure out moves
995  for(int ham=am+1;ham<=maxam(cpl);ham++) {
996  double ds=cpl[am].start-cpl[ham].start;
997  double de=cpl[ham].end-cpl[am].end;
998 
999  smove=std::max(ds,smove);
1000  emove=std::max(de,emove);
1001  if(cpl[ham].tol<Yam.tol)
1002  Yam.tol=cpl[ham].tol;
1003  }
1004 
1005  Timer tam;
1006 
1007  // Move ends by
1008  double tmove=smove+emove;
1009  if(tmove==0.0 && Yam.tol==cpl[am].tol)
1010  continue;
1011 
1012  // Print out info bar
1013  if(!printed) {
1014  printf("\n\n%s\n",print_bar("Subset check",'-').c_str());
1015  fflush(stdout);
1016  printed=true;
1017  }
1018 
1019  // Current width is
1020  double curw=cpl[am].end-cpl[am].start;
1021  // Necessary width is
1022  double nw=curw+tmove;
1023 
1024  // Get real width
1025  double realw(nw);
1026  arma::vec exps=span_width(am,Yam.tol,realw,Yam.exps.size());
1027 
1028  // Adjust profile
1029  int nadd=(int) (exps.size()-cpl[am].exps.size());
1030  double sreal=(realw-curw)*smove/tmove;
1031  double ereal=(realw-curw)*emove/tmove;
1032 
1033  if(smove>0.0)
1034  printf("%c lower limit should be moved by % .3f; moved by % .3f.\n",shell_types[am],smove,sreal);
1035  if(emove>0.0)
1036  printf("%c upper limit should be moved by % .3f; moved by % .3f.\n",shell_types[am],emove,ereal);
1037  if(cpl[am].tol != Yam.tol)
1038  printf("%c tolerance reduced from %e to %e.\n",shell_types[am],cpl[am].tol,Yam.tol);
1039  printf("%i exponents added to %c shell (%s)\n\n",nadd,shell_types[am],tam.elapsed().c_str());
1040  fflush(stdout);
1041 
1042  Yam.start-=sreal;
1043  Yam.end+=ereal;
1044  Yam.exps=move_exps(exps,Yam.start);
1045  upd=true;
1046  cpl[am]=Yam;
1047  }
1048 
1049  if(upd) {
1050  // Update current value
1051  std::vector< std::vector<coprof_t> > cplv(1,cpl);
1052  std::vector<ValueType> rval(compute_values(cplv));
1053  curval=rval[0];
1054  print_value(curval,"Updated value");
1055  print_limits(cpl,"Updated limits");
1056 
1057  printf("%s\n\n",print_bar("End check",'-').c_str());
1058  }
1059 
1060  return upd;
1061  }
1062 
1075  double scan_profile(std::vector<coprof_t> & cpl, ValueType & curval, int npoints, double dpol, double tol, bool allam=true, bool bothsd=true, bool forcesub=false) {
1076 
1077  printf("\n\n%s\n",print_bar("SHELL STABILITY").c_str());
1078  fflush(stdout);
1079 
1080  // Get elemental libraries
1081  const std::vector<ElementBasisSet> els=form_basis(cpl).get_elements();
1082 
1083  // Perform scan
1084  Timer t;
1085 
1086  printf("Scanning for shell stability.\n");
1087  printf("Using %.3f points per exponent interval.\n\n",dpol);
1088  fflush(stdout);
1089 
1090  // Trial exponents and am
1091  std::vector<double> trexp;
1092  std::vector<int> tram;
1093 
1094  // Step sizes
1095  arma::vec spacing(maxam(cpl)+1);
1096  for(int scanam=0;scanam<=maxam(cpl);scanam++) {
1097  // Compute step size
1098  double step;
1099  {
1100  // We now have nx exponents. Because the extension has already
1101  // been done, expanding the profile by the step size
1102  // w(nx+1)-w(nx) is not enough to get a notable change in the
1103  // mog. Thus, the profile will get at least two new
1104  // functions. So, we can use a bigger step size w(nx+2)-w(nx+1)
1105  // here.
1106 
1107  arma::uword nexp=cpl[scanam].exps.n_elem+1;
1108 
1109  double nextw;
1110  maxwidth_exps_table(scanam,cpl[scanam].tol,nexp+1,nextw);
1111  double curw;
1112  maxwidth_exps_table(scanam,cpl[scanam].tol,nexp,curw);
1113  step=nextw-curw;
1114 
1115  // Store spacing
1116  spacing(scanam)=step;
1117  }
1118 
1119  printf("Spacing between points on %c shell is %.5f\n",shell_types[scanam],step);
1120  fflush(stdout);
1121 
1122  // Generate points
1123  int nmin=0;
1124  int nmax=(int) std::ceil(npoints*dpol);
1125  for(int n=nmax;n>=nmin;n--) {
1126  trexp.push_back(cpl[scanam].start-n*step/dpol);
1127  tram.push_back(scanam);
1128  }
1129  for(int n=nmin;n<=nmax;n++) {
1130  trexp.push_back(cpl[scanam].end+n*step/dpol);
1131  tram.push_back(scanam);
1132  }
1133  }
1134 
1135  // Trial basis sets
1136  std::vector< BasisSetLibrary > trbas;
1137  for(size_t iexp=0;iexp<trexp.size();iexp++) {
1138  // Trial basis set
1139  BasisSetLibrary baslib;
1140  // Add the function
1141  for(size_t iel=0;iel<els.size();iel++) {
1142  ElementBasisSet elbas(els[iel]);
1143 
1144  FunctionShell sh(tram[iexp]);
1145  sh.add_exponent(1.0,std::pow(10.0,trexp[iexp]));
1146  elbas.add_function(sh);
1147 
1148  // and the element to the new library
1149  baslib.add_element(elbas);
1150  }
1151  trbas.push_back(baslib);
1152  }
1153 
1154  // Sanity check
1155  if(!trbas.size()) {
1156  printf("No shells to scan stability of.\n");
1157  return 0.0;
1158  }
1159 
1160  // Compute trial values
1161  std::vector<ValueType> trvals=compute_values(trbas);
1162  if(trvals.size()!=trbas.size()) {
1163  std::ostringstream oss;
1164  oss << "Error - requested computation of " << trbas.size() << " values but got only " << trvals.size() << "!\n";
1165  throw std::runtime_error(oss.str());
1166  }
1167 
1168  // and the mogs
1169  printf("%11s %2s %8s %12s\n","trial ","am","exponent","mog");
1170  fflush(stdout);
1171  arma::vec mogs(trvals.size());
1172  for(size_t iexp=0;iexp<mogs.n_elem;iexp++) {
1173  mogs(iexp)=compute_mog(trvals[iexp],curval,0.0);
1174  printf("%5i/%-5i %-2c % 7.5f %e\n",(int) iexp+1, (int) mogs.n_elem, shell_types[tram[iexp]], trexp[iexp], mogs(iexp));
1175  fflush(stdout);
1176  }
1177 
1178  // Decompose mogs per angular momentum
1179  std::vector<arma::uvec> amidx;
1180  std::vector<arma::vec> ammog;
1181  mog_decompose(tram,mogs,amidx,ammog);
1182  printf("\n\t%2s %8s %12s\n","am","exponent","mog");
1183  for(size_t am=0;am<amidx.size();am++)
1184  if(amidx[am].n_elem) {
1185  // Find location of maximum
1186  arma::uword maxind;
1187  double maxmog=ammog[am].max(maxind);
1188  printf("\t%-2c % 7.5f %e\n",shell_types[am],trexp[amidx[am](maxind)],maxmog);
1189  } else
1190  printf("\t%-2c %8s %e\n",shell_types[am],"",0.0);
1191  fflush(stdout);
1192 
1193  {
1194  // Save out the results of the mog scan
1195  for(int am=0;am<=maxam(cpl);am++) {
1196  std::ostringstream fname;
1197  fname << "scanmog_" << politer << "_" << shell_types[am] << ".dat";
1198 
1199  // Collect exponent values
1200  arma::vec expval(ammog[am].n_elem);
1201  for(size_t i=0;i<ammog[am].n_elem;i++)
1202  expval(i)=trexp[amidx[am][i]];
1203 
1204  arma::mat savemog(expval.n_elem,2);
1205  savemog.col(0)=expval;
1206  savemog.col(1)=ammog[am];
1207  savemog.save(fname.str(),arma::raw_ascii);
1208  }
1209  }
1210 
1211  // Find maximum mog
1212  arma::uword imax;
1213  double maxmog=mogs.max(imax);
1214 
1215  printf("\n");
1216  printf("%11s %-2c % 7.5f %e\n","max mog",shell_types[tram[imax]],trexp[imax],maxmog);
1217  fflush(stdout);
1218 
1219  double moved;
1220 
1221  if(maxmog>=tol) {
1222  printf("\n");
1223 
1224  // Loop over angular momentum
1225  for(int scanam=0;scanam<=maxam(cpl);scanam++) {
1226  Timer tam;
1227 
1228  if(!allam && ((arma::uword) scanam != imax))
1229  // Don't examine this trial
1230  continue;
1231 
1232  if(!ammog[scanam].n_elem)
1233  continue;
1234 
1235  // Indices of maxima
1236  std::vector<arma::uword> idxlist;
1237 
1238  if(bothsd) {
1239  // Gather trials for steep and diffuse end
1240  std::vector<size_t> stidx, dfidx;
1241  std::vector<double> stmog, dfmog;
1242  for(size_t itr=0;itr<trexp.size();itr++)
1243  if(tram[itr]==scanam) {
1244  if(trexp[itr]<=cpl[scanam].start) {
1245  dfidx.push_back(itr);
1246  dfmog.push_back(mogs[itr]);
1247  } else {
1248  stidx.push_back(itr);
1249  stmog.push_back(mogs[itr]);
1250  }
1251  }
1252 
1253  // Form armadillo vectors
1254  arma::vec stm(arma::conv_to<arma::vec>::from(stmog));
1255  arma::vec dfm(arma::conv_to<arma::vec>::from(dfmog));
1256 
1257  // Get maximum mogs and indices
1258  arma::uword stind, dfind;
1259  double stmax=stm.max(stind);
1260  double dfmax=dfm.max(dfind);
1261  if(stmax>=tol)
1262  idxlist.push_back(stidx[stind]);
1263  if(dfmax>=tol)
1264  idxlist.push_back(dfidx[dfind]);
1265  } else {
1266  // Get maximum mog for this am
1267  arma::uword amind;
1268  double ammax=ammog[scanam].max(amind);
1269 
1270  // Index of trial is
1271  if(ammax>=tol)
1272  idxlist.push_back(amidx[scanam](amind));
1273  }
1274 
1275  // Adjust profile
1276  for(size_t i=0;i<idxlist.size();i++) {
1277  imax=idxlist[i];
1278 
1279  int minam=scanam;
1280  int maxam=scanam;
1281  if(forcesub)
1282  minam=0;
1283 
1284  for(int am=minam;am<=maxam;am++) {
1285  moved=0.0;
1286 
1287  if(trexp[imax] < cpl[am].start) {
1288  // Current width is
1289  double curw=cpl[am].end-cpl[am].start;
1290  // Necessary width is
1291  double nw=cpl[am].end-trexp[imax];
1292  // Get real width
1293  double realw(nw);
1294  arma::vec exps=span_width(am,cpl[am].tol,realw,cpl[am].exps.size());
1295 
1296  // Adjust profile
1297  cpl[am].start-=realw-curw;
1298  cpl[am].exps=move_exps(exps,cpl[am].start);
1299  moved=-(realw-curw);
1300 
1301  } else if(trexp[imax] > cpl[am].end) {
1302  // Current width is
1303  double curw=cpl[am].end-cpl[am].start;
1304  // Necessary width is
1305  double nw=trexp[imax]-cpl[am].start;
1306  // Get real width
1307  double realw(nw);
1308  arma::vec exps=span_width(am,cpl[am].tol,realw,cpl[am].exps.size());
1309 
1310  // Adjust profile
1311  cpl[am].end+=realw-curw;
1312  cpl[am].exps=move_exps(exps,cpl[am].start);
1313  moved=+(realw-curw);
1314 
1315  } else if(am==scanam) {
1316  throw std::runtime_error("Possible bug in scan_limits - maximum inside profile!\n");
1317  }
1318 
1319  if(moved>0.0)
1320  printf("%c upper limit should be moved by % .3f (% .3f spacings). (%s)\n",shell_types[am],moved,moved/spacing(am),tam.elapsed().c_str());
1321  else if(moved<0.0)
1322  printf("%c lower limit should be moved by % .3f (% .3f spacings). (%s)\n",shell_types[am],-moved,-moved/spacing(am),tam.elapsed().c_str());
1323  fflush(stdout);
1324  }
1325  }
1326  }
1327  printf("\n");
1328 
1329  // Update current value
1330  std::vector< std::vector<coprof_t> > hlp(1,cpl);
1331  std::vector<ValueType> hlpvals(compute_values(hlp));
1332  if(hlpvals.size()!=1) {
1333  std::ostringstream oss;
1334  oss << "Error - requested computation of " << 1 << " values but got only " << hlpvals.size() << "!\n";
1335  throw std::runtime_error(oss.str());
1336  }
1337  curval=hlpvals[0];
1338  }
1339 
1340  printf("Stability scan done in %s.\n\n",t.elapsed().c_str());
1341  fflush(stdout);
1342 
1343  return maxmog;
1344  }
1345 
1347  double extend_profile(std::vector<coprof_t> & cpl, ValueType & curval, double tau, int next=3, int nxadd=1, bool forcesub=false) {
1348  printf("\n\n%s\n",print_bar("PROFILE EXTENSION").c_str());
1349  printf("Final tolerance is %e.\n",tau);
1350  fflush(stdout);
1351 
1352  if(tau<=0.0)
1353  throw std::runtime_error("Tolerance must be positive!\n");
1354  if(next<2)
1355  throw std::runtime_error("The NExt parameter must be at least two.\n");
1356 
1357  while(true) {
1358  Timer ttot;
1359 
1360  // Separator
1361  printf("\n");
1362  fflush(stdout);
1363 
1364  // Trial profiles
1365  std::vector< std::vector<coprof_t> > trials;
1366  // Descriptions
1367  std::vector< std::string > descr;
1368  // Trial angular momentum
1369  std::vector< int > tram;
1370 
1371  for(int am=0;am<=maxam(cpl);am++) {
1372  Timer t;
1373 
1374  // Get exponents
1375  printf("Determining exponents for %c shell ... ",shell_types[am]);
1376  fflush(stdout);
1377 
1378  double step;
1379  double width;
1380  arma::vec exps=maxwidth_exps_table(am,cpl[am].tol,cpl[am].exps.size()+nxadd,width);
1381 
1382  step=width-(cpl[am].end-cpl[am].start);
1383  printf("step size is %7.5f (%s).\n",step,t.elapsed().c_str());
1384  fflush(stdout);
1385  t.set();
1386 
1387  // Form trials.
1388  for(int itr=0;itr<next;itr++) {
1389  // Displacement fraction is
1390  double d=itr*1.0/(next-1);
1391 
1392  std::vector<coprof_t> trcpl(cpl);
1393  trcpl[am].start=cpl[am].start-d*step;
1394  trcpl[am].end=cpl[am].end+(1.0-d)*step;
1395  trcpl[am].exps=move_exps(exps,trcpl[am].start);
1396 
1397  // Add the trial to the stack
1398  trials.push_back(trcpl);
1399  tram.push_back(am);
1400 
1401  char msg[200];
1402  if(itr==0)
1403  sprintf(msg,"Moved ending point of %c shell by %.3f",shell_types[am],step);
1404  else if(itr==next-1)
1405  sprintf(msg,"Moved starting point of %c shell by %.3f",shell_types[am],step);
1406  else
1407  sprintf(msg,"Moved starting point of %c shell by %.3f and ending point by %.3f",shell_types[am],d*step,(1.0-d)*step);
1408  descr.push_back(msg);
1409  }
1410  }
1411 
1412  // Sanity check
1413  if(!trials.size()) {
1414  printf("No shells can be extended.\n");
1415  return 0.0;
1416  }
1417 
1418  // Compute values
1419  std::vector<ValueType> trvals=compute_values(trials);
1420  if(trvals.size()!=trials.size()) {
1421  std::ostringstream oss;
1422  oss << "Error - requested computation of " << trials.size() << " values but got only " << trvals.size() << "!\n";
1423  throw std::runtime_error(oss.str());
1424  }
1425 
1426  // Compute mogs
1427  arma::vec mogs(trvals.size());
1428  for(size_t i=0;i<trvals.size();i++)
1429  mogs(i)=compute_mog(trvals[i],curval,0.0);
1430 
1431  // Decompose mogs per angular momentum
1432  std::vector<arma::uvec> amidx;
1433  std::vector<arma::vec> ammog;
1434  mog_decompose(tram,mogs,amidx,ammog);
1435  printf("\n\t%2s %12s\n","am","mog");
1436  for(size_t am=0;am<amidx.size();am++)
1437  if(amidx[am].n_elem) {
1438  printf("\t%-2c %e\n",shell_types[am],arma::max(ammog[am]));
1439  } else
1440  printf("\t%-2c %e\n",shell_types[am],0.0);
1441  fflush(stdout);
1442 
1443  // Figure out maximal mog
1444  arma::uword maxind;
1445  double maxmog=mogs.max(maxind);
1446 
1447  // Converged?
1448  if(maxmog < tau) {
1449  printf("Maximal mog is %e, converged.\n\n",maxmog);
1450  print_value(curval,"Final value");
1451  print_limits(cpl,"Final limits");
1452  fflush(stdout);
1453 
1454  if(forcesub) {
1455  bool upd=comply_subset(cpl,curval);
1456  if(upd)
1457  continue;
1458  }
1459 
1460  return maxmog;
1461 
1462  } else {
1463  // Still stuff to add
1464 
1465  // Update profile
1466  cpl=trials[maxind];
1467  // Update current value
1468  curval=trvals[maxind];
1469 
1470  // Print message
1471  printf("%s, range is now % .3f ... % .3f (%2i funcs), tol = %e, mog %e (%s).\n",descr[maxind].c_str(),cpl[tram[maxind]].start,cpl[tram[maxind]].end,(int) cpl[tram[maxind]].exps.size(),cpl[tram[maxind]].tol,maxmog,ttot.elapsed().c_str());
1472  fflush(stdout);
1473 
1474  // Update references
1475  update_reference(cpl);
1476  }
1477  }
1478  }
1479 
1480 
1482  double tighten_profile(std::vector<coprof_t> & cpl, ValueType & curval, double tau, int nxadd=1, bool forcesub=false) {
1483  printf("\n\n%s\n",print_bar("PROFILE TIGHTENING").c_str());
1484  printf("Final tolerance is %e.\n",tau);
1485  fflush(stdout);
1486 
1487  while(true) {
1488  Timer ttot;
1489 
1490  // Separator
1491  printf("\n");
1492  fflush(stdout);
1493 
1494  // Form trials
1495  std::vector< std::vector<coprof_t> > trials;
1496  std::vector< int > tram;
1497  for(int am=0;am<=maxam(cpl);am++) {
1498  Timer t;
1499 
1500  // Form trials
1501  std::vector<coprof_t> add(cpl);
1502 
1503  // Get exponents
1504  printf("Determining exponents for %c shell ... ",shell_types[am]);
1505  fflush(stdout);
1506 
1507  add[am].exps=optimize_completeness(am,add[am].start,add[am].end,add[am].exps.size()+nxadd,add[am].tol);
1508  printf("(%s)\n",t.elapsed().c_str());
1509  t.set();
1510 
1511  if (add[am].tol>=MINTAU) {
1512  trials.push_back(add);
1513  tram.push_back(am);
1514  }
1515  }
1516 
1517  // Sanity check
1518  if(!trials.size()) {
1519  printf("No shells can be tightened.\n");
1520  return 0.0;
1521  }
1522 
1523  // Compute values
1524  std::vector<ValueType> trvals=compute_values(trials);
1525  if(trvals.size()!=trials.size()) {
1526  std::ostringstream oss;
1527  oss << "Error - requested computation of " << trials.size() << " values but got only " << trvals.size() << "!\n";
1528  throw std::runtime_error(oss.str());
1529  }
1530 
1531  // and mogs
1532  arma::vec trmog(trvals.size());
1533  for(size_t i=0;i<trvals.size();i++)
1534  trmog(i)=compute_mog(trvals[i],curval,0.0);
1535 
1536  // Decompose mogs per angular momentum
1537  std::vector<arma::uvec> amidx;
1538  std::vector<arma::vec> ammog;
1539  mog_decompose(tram,trmog,amidx,ammog);
1540  printf("\n\t%2s %12s\n","am","mog");
1541  for(size_t am=0;am<amidx.size();am++)
1542  if(amidx[am].n_elem) {
1543  printf("\t%-2c %e\n",shell_types[am],arma::max(ammog[am]));
1544  } else
1545  printf("\t%-2c %e\n",shell_types[am],0.0);
1546  fflush(stdout);
1547 
1548  // Figure out maximal mog
1549  arma::uword maxind;
1550  double maxmog=trmog.max(maxind);
1551 
1552  // Converged?
1553  if(maxmog < tau) {
1554  printf("Maximal mog is %e, converged.\n\n",maxmog);
1555  print_value(curval,"Final value");
1556  print_limits(cpl,"Final tolerances");
1557  fflush(stdout);
1558 
1559  if(forcesub)
1560  comply_subset(cpl,curval);
1561 
1562  return maxmog;
1563 
1564  } else {
1565  // Still stuff to add
1566 
1567  // Update profile
1568  cpl=trials[maxind];
1569  // Update current value
1570  curval=trvals[maxind];
1571 
1572  printf("Added exponent to %c shell, range is now % .3f ... % .3f (%2i funcs), tol = %e, mog = %e (%s).\n",shell_types[maxind],cpl[maxind].start,cpl[maxind].end,(int) cpl[maxind].exps.size(),cpl[maxind].tol,maxmog,ttot.elapsed().c_str());
1573  fflush(stdout);
1574  }
1575 
1576  // Update references
1577  update_reference(cpl);
1578  }
1579  }
1580 
1582  virtual int atom_am() const=0;
1583 
1585  virtual void print_value(const ValueType & value, std::string msg)=0;
1586 
1588  void find_cbs_limit(std::vector<coprof_t> & cpl, ValueType & curval, double cotol, double minpol, double maxpol, double dpol, bool extend=true, int next=3, bool scan=true, int nscan=5, bool polinterp=true, int nxpol=1, bool doadd=true, int nxext=1, int am_max=max_am, bool cbsinterp=true, double cbsthr=0.0, double delta=0.9, bool allam=true, bool bothsd=true, bool forcesub=false) {
1589  // Amount of polarization shells
1590  int npol=maxam(cpl)-atom_am();
1591 
1592  // Check sanity
1593  if(!extend && !scan)
1594  throw std::runtime_error("You need to check for profile expansion either by scans or extension, or both!\n");
1595 
1596  // Compute initial value
1597  {
1598  std::vector< std::vector<coprof_t> > hlp(1,cpl);
1599  std::vector<ValueType> hlpvals(compute_values(hlp));
1600  if(hlpvals.size()!=1) {
1601  std::ostringstream oss;
1602  oss << "Error - requested computation of " << 1 << " values but got only " << hlpvals.size() << "!\n";
1603  throw std::runtime_error(oss.str());
1604  }
1605  curval=hlpvals[0];
1606  }
1607  print_value(curval,"Starting value");
1608 
1609  // Compute mog of next polarization shell
1610  double tau;
1611  {
1612  std::vector<coprof_t> polcpl(cpl);
1613  ValueType polval(curval);
1614  tau=check_polarization(polcpl,polval,am_max,minpol,maxpol,dpol,polinterp,cotol,nxpol,forcesub);
1615  }
1616 
1617  // Mogs of added polarization shells
1618  std::vector<double> polmogs;
1619 
1620  while(true) {
1621  // Extend existing shells
1622  double extmog=0.0;
1623 
1624  if(extend)
1625  extmog=extend_profile(cpl,curval,tau,next,nxext,forcesub);
1626  // Tighten existing shells
1627  if(doadd) {
1628  double amog=tighten_profile(cpl,curval,tau,1,forcesub);
1629  extmog=std::max(extmog,amog);
1630  }
1631 
1632  // Compute mog of further polarization shell
1633  std::vector<coprof_t> polcpl(cpl);
1634  ValueType polval(curval);
1635  double polmog=check_polarization(polcpl,polval,am_max,minpol,maxpol,dpol,polinterp,cotol,nxpol,forcesub);
1636 
1637  // What is the maximum mog?
1638  if(extmog >= std::max(polmog,cbsthr)) {
1639  // Need to expand existing shells.
1640  tau=exp((1.0-delta)*log(extmog) + delta*log(std::max(polmog,cbsthr)));
1641 
1642  } else {
1643  // Here extmog < max(polmog,cbsthr).
1644  // Going to add new polarization shell in next iteration; we are converged here.
1645 
1646  if(scan) {
1647  // Before adding new polarization shell, scan the stability of the existing shells.
1648  std::vector<coprof_t> scancpl(cpl);
1649  ValueType scanval(curval);
1650  double scanmog=scan_profile(scancpl,scanval,nscan,dpol,std::max(polmog,cbsthr),allam,bothsd,forcesub);
1651 
1652  print_value(scanval,"Compound value");
1653  print_limits(scancpl,"Compound limits");
1654  fflush(stdout);
1655 
1656  // Did the scan fail?
1657  bool scanfail;
1658 
1659  // Update extension mog
1660  extmog=std::max(extmog,scanmog);
1661 
1662  if(scanmog>=std::max(polmog,cbsthr)) {
1663  // Instability detected, real mog is
1664  double mog=compute_mog(scanval,curval,0.0);
1665 
1666  // Sanity check
1667  if(mog>0.0) {
1668  scanfail=false;
1669  cpl=scancpl;
1670  curval=scanval;
1671  printf("\n\nInstability detected with real mog = %e, restarting extension.\n",mog);
1672 
1673  } else {
1674  // Calculation didn't converge :(
1675 
1676  printf("\n\nInstability detected, but compound calculation failed to converge.\n");
1677  printf("Forming new, individual trials.\n");
1678 
1679  // Form new trials for the changed ams
1680  std::vector< std::vector<coprof_t> > trials;
1681  std::vector<int> tram;
1682  for(int am=0;am<=maxam(cpl);am++)
1683  if(scancpl[am].exps.size() != cpl[am].exps.size()) {
1684  // Form trial
1685  std::vector<coprof_t> tr(cpl);
1686  tr[am]=scancpl[am];
1687  // and add it to the stack
1688  trials.push_back(tr);
1689  tram.push_back(am);
1690  }
1691 
1692  // Compute new values
1693  std::vector<ValueType> trvals=compute_values(trials);
1694  // and the corresponding mogs
1695  arma::vec trmog(trvals.size());
1696  for(size_t i=0;i<trvals.size();i++)
1697  trmog[i]=compute_mog(trvals[i],curval,0.0);
1698 
1699  printf("\n\t%2s %12s\n","am","mog");
1700  for(size_t i=0;i<tram.size();i++)
1701  printf("\t%-2c %e\n",shell_types[tram[i]],trmog[i]);
1702  fflush(stdout);
1703 
1704  // Pick out maximum
1705  arma::uword maxind;
1706  double maxmog=trmog.max(maxind);
1707  if(maxmog>0.0) {
1708  // Switch values
1709  cpl=trials[maxind];
1710  curval=trvals[maxind];
1711  scanfail=false;
1712  } else
1713  scanfail=true;
1714  }
1715 
1716  print_value(curval,"Current value");
1717  print_limits(cpl,"Current limits");
1718 
1719  if(!scanfail)
1720  continue;
1721  }
1722  }
1723 
1724  // Save polarization mog
1725  if(maxam(polcpl)>maxam(cpl))
1726  polmogs.push_back(polmog);
1727 
1728  // Switch to polarized basis.
1729  cpl=polcpl;
1730  curval=polval;
1731  npol++;
1732 
1733  printf("\n\n%s\n",print_bar("CONVERGENCE ACHIEVED",'#').c_str());
1734  printf("\nConverged: extension mog was %e, polarization mog %e.\n",extmog,polmog);
1735 
1736  printf("\nFinal composition for %i polarization shells (tau = %e):\n",npol,polmog);
1737  print_value(curval,"Current value");
1738  print_limits(cpl,"Current limits");
1739  fflush(stdout);
1740 
1741  // Save basis set
1742  {
1743  BasisSetLibrary baslib=form_basis(cpl);
1744  char fname[180];
1745  sprintf(fname,"un-co-ref-%i.gbs",npol);
1746  baslib.save_gaussian94(fname);
1747 
1748  // Save completeness range
1749  sprintf(fname,"completeness-%i.dat",npol);
1750  save_limits(cpl,fname);
1751  }
1752 
1753  // Converged?
1754  if(cbsthr>0.0) {
1755  // CBS threshold overrides am convergence check
1756  if(std::max(extmog,polmog)<cbsthr) {
1757  printf("\nCBS threshold reached: extension mog was %e, polarization mog %e.\n",extmog,polmog);
1758  break;
1759  }
1760  } else if(maxam(cpl)==am_max) {
1761  // Use current polarization mog to extend profile
1762  tau=polmog;
1763  // Interpolate to next polarization mog?
1764  if(cbsinterp && polmogs.size()>=2) {
1765  // Mog of highest shell is
1766  double hmog=polmogs[polmogs.size()-1];
1767  // and the one before that
1768  double lmog=polmogs[polmogs.size()-2];
1769 
1770  // Fit slope: log(mog) = a (delta L) + b,
1771  // so the next shell occurs at
1772  tau=hmog*hmog/lmog;
1773  }
1774 
1775  while(true) {
1776  // Check shell extension
1777  extend_profile(cpl,curval,tau,next,nxext,forcesub);
1778  if(doadd)
1779  tighten_profile(cpl,curval,tau,1,forcesub);
1780 
1781  double scanmog=0.0;
1782  if(scan) {
1783  // Before adding new polarization shell, scan the stability of the existing shells.
1784  std::vector<coprof_t> scancpl(cpl);
1785  ValueType scanval(curval);
1786  scanmog=scan_profile(scancpl,scanval,nscan,dpol,tau,allam,bothsd,forcesub);
1787 
1788  if(scanmog>=tau) {
1789  // Instability detected, real mog is
1790  double mog=compute_mog(scanval,curval,0.0);
1791 
1792  cpl=scancpl;
1793  curval=scanval;
1794  printf("\n\nInstability detected with real mog = %e, restarting extension.\n",mog);
1795 
1796  print_value(curval,"Current value");
1797  print_limits(cpl,"Current limits");
1798 
1799  // Restart extension
1800  continue;
1801  }
1802  }
1803 
1804  // Break last shell extension loop
1805  break;
1806  }
1807 
1808  // Break polarization loop
1809  break;
1810  }
1811 
1812  // Reduce tau. Sanity check for no functions
1813  if(extmog>0.0)
1814  tau=std::min(tau,std::max(extmog,cbsthr));
1815  }
1816  }
1817 
1818  print_value(curval,"\nFinal value");
1819  print_limits(cpl,"Final composition:");
1820  fflush(stdout);
1821 
1822  // Save basis set
1823  {
1824  BasisSetLibrary baslib=form_basis(cpl);
1825  char fname[180];
1826  sprintf(fname,"un-co-ref.gbs");
1827  baslib.save_gaussian94(fname);
1828  }
1829 
1830  // Save completeness range
1831  save_limits(cpl,"completeness.dat");
1832  }
1833 
1835  double reduce_profile(std::vector<coprof_t> & cpl, ValueType & curval, const ValueType & refval, double tol=0.0, int nred=3, bool saveall=false, bool allelectron=true, bool forcesub=false) {
1836  printf("\n\n%s\n",print_bar("PROFILE REDUCTION").c_str());
1837  if(tol==0.0)
1838  printf("Running until would drop last function of %c shell.\n",shell_types[maxam(cpl)]);
1839  else
1840  printf("Running until mog >= %e.\n",tol);
1841  fflush(stdout);
1842 
1843  if(nred<2)
1844  throw std::runtime_error("The NRed parameter must be at least two.\n");
1845 
1846  // Do the reduction.
1847  double minmog=0.0;
1848 
1849  while(true) {
1850  Timer ttot;
1851 
1852  // Separator
1853  printf("\n");
1854  fflush(stdout);
1855 
1856  // Form trial profiles
1857  std::vector< std::vector<coprof_t> > trials;
1858  std::vector< std::string > descr;
1859  std::vector< int > tram;
1860 
1861  for(int am=0;am<=maxam(cpl);am++) {
1862  if(!cpl[am].exps.size())
1863  continue;
1864 
1865  // Sanity check
1866  bool tryred=true;
1867  if(allelectron)
1868  // Only do composition check for all-electron calculations
1869  for(int sam=am+1;sam<=maxam(cpl);sam++) {
1870  // Check that the reduced shell would have more functions than the ones above.
1871  if(cpl[am].exps.size()-1 <= cpl[sam].exps.size()) {
1872  printf("%c shell limited due to %c shell.\n",shell_types[am],shell_types[am+1]);
1873  fflush(stdout);
1874  tryred=false;
1875  break;
1876  }
1877  }
1878  if(!tryred)
1879  continue;
1880 
1881  // Form trials.
1882  if(cpl[am].exps.size()==1) {
1883  std::vector<coprof_t> trcpl(cpl);
1884  trcpl[am].exps.clear();
1885  trcpl[am].start=0.0;
1886  trcpl[am].end=0.0;
1887 
1888  // Add the trial to the stack
1889  trials.push_back(trcpl);
1890  tram.push_back(am);
1891 
1892  char msg[200];
1893  sprintf(msg,"Dropped last exponent of %c shell",shell_types[am]);
1894  descr.push_back(msg);
1895 
1896  } else {
1897 
1898  Timer t;
1899 
1900  // Get exponents
1901  printf("Determining exponents for %c shell ... ",shell_types[am]);
1902  fflush(stdout);
1903 
1904  // Keep limits but drop flatness
1905  std::vector<coprof_t> delcpl(cpl);
1906  delcpl[am].exps=optimize_completeness(am,delcpl[am].start,delcpl[am].end,delcpl[am].exps.size()-1,delcpl[am].tol);
1907 
1908  // Keep flatness but change limits
1909  double step;
1910  double width;
1911  arma::vec exps=maxwidth_exps_table(am,cpl[am].tol,cpl[am].exps.size()-1,width);
1912  step=(cpl[am].end-cpl[am].start)-width;
1913  printf("step size is %7.5f (%s).\n",step,t.elapsed().c_str());
1914  fflush(stdout);
1915  t.set();
1916 
1917  // Form trials
1918  for(int itr=0;itr<nred;itr++) {
1919  // Displacement fraction is
1920  double d=itr*1.0/(nred-1);
1921 
1922  std::vector<coprof_t> trcpl(cpl);
1923  trcpl[am].start=cpl[am].start+(1.0-d)*step;
1924  trcpl[am].end=cpl[am].end-d*step;
1925  trcpl[am].exps=move_exps(exps,trcpl[am].start);
1926 
1927  bool dotrial=true;
1928  if(forcesub) {
1929  // Sanity check trial
1930  for(int ham=am+1;ham<=maxam(cpl);ham++) {
1931  if(cpl[ham].start<trcpl[am].start)
1932  dotrial=false;
1933  if(cpl[ham].end>trcpl[am].end)
1934  dotrial=false;
1935  }
1936  }
1937 
1938  // Add the trial to the stack
1939  if(dotrial) {
1940  trials.push_back(trcpl);
1941  tram.push_back(am);
1942 
1943  char msg[200];
1944  if(itr==0)
1945  sprintf(msg,"Moved starting point of %c shell by %.3f",shell_types[am],step);
1946  else if(itr==nred-1)
1947  sprintf(msg,"Moved ending point of %c shell by %.3f",shell_types[am],step);
1948  else
1949  sprintf(msg,"Moved starting point of %c shell by %.3f and ending point by %.3f",shell_types[am],(1.0-d)*step,d*step);
1950  descr.push_back(msg);
1951  }
1952  }
1953 
1954  // Keep limits but drop a function. The deviation from
1955  // completeness should increase monotonically with increasing am.
1956  bool dodrop=true;
1957  for(int ham=am+1;ham<=maxam(cpl);ham++)
1958  // Need to have at least two exponents for the restriction
1959  // for the deviation from completeness to make sense
1960  if(delcpl[ham].exps.size()>1 && delcpl[am].tol > delcpl[ham].tol)
1961  dodrop=false;
1962  if(dodrop) {
1963  trials.push_back(delcpl);
1964 
1965  char msg[200];
1966  sprintf(msg,"Dropped exponent from %c shell",shell_types[am]);
1967  descr.push_back(msg);
1968  tram.push_back(am);
1969  }
1970  }
1971  }
1972 
1973  // Empty basis set - nothing to remove
1974  if(!trials.size())
1975  return DBL_MAX;
1976 
1977  // Compute values
1978  std::vector<ValueType> trvals=compute_values(trials);
1979  if(trvals.size()!=trials.size()) {
1980  std::ostringstream oss;
1981  oss << "Error - requested computation of " << trials.size() << " values but got only " << trvals.size() << "!\n";
1982  throw std::runtime_error(oss.str());
1983  }
1984 
1985  // and mogs
1986  arma::vec trmog(trvals.size());
1987  for(size_t i=0;i<trvals.size();i++)
1988  trmog[i]=compute_mog(trvals[i],refval,DBL_MAX);
1989 
1990  // Decompose mogs per angular momentum
1991  std::vector<arma::uvec> amidx;
1992  std::vector<arma::vec> ammog;
1993  mog_decompose(tram,trmog,amidx,ammog);
1994  printf("\n\t%2s %12s\n","am","mog");
1995  for(size_t am=0;am<amidx.size();am++)
1996  if(amidx[am].n_elem) {
1997  printf("\t%-2c %e\n",shell_types[am],arma::min(ammog[am]));
1998  } else
1999  printf("\t%-2c %e\n",shell_types[am],0.0);
2000  fflush(stdout);
2001 
2002  // Figure out minimal mog.
2003  arma::uword minind;
2004  minmog=trmog.min(minind);
2005 
2006  // Check if a higher am can be achieved with the same mog
2007  for(size_t i=0;i<trmog.n_elem;i++)
2008  if(trmog[i] == minmog && tram[i] > tram[minind])
2009  minind=i;
2010 
2011  // Converged?
2012  if(tol==0.0 && ((int) tram[minind] == maxam(cpl)) && (cpl[tram[minind]].exps.size()==1)) {
2013  printf("Would remove last exponent of %c shell with mog %e, converged.\n\n",shell_types[tram[minind]],minmog);
2014  print_value(curval,"Final value");
2015  print_limits(cpl,"Final limits");
2016  fflush(stdout);
2017 
2018  return minmog;
2019 
2020  } else if(tol>0.0 && minmog>tol) {
2021 
2022  printf("Minimal mog is %e for %c shell, converged.\n\n",minmog,shell_types[tram[minind]]);
2023  print_value(curval,"Final value");
2024  print_limits(cpl,"Final limits");
2025  fflush(stdout);
2026  return minmog;
2027 
2028  } else {
2029  // Still stuff to remove.
2030 
2031  // Update profile
2032  cpl=trials[minind];
2033  // Update current value
2034  curval=trvals[minind];
2035  // Print message
2036  printf("%s, range is now % .3f ... % .3f (%2i funcs), tol = %e, mog %e (%s).\n",descr[minind].c_str(),cpl[tram[minind]].start,cpl[tram[minind]].end,(int) cpl[tram[minind]].exps.size(),cpl[tram[minind]].tol,minmog,ttot.elapsed().c_str());
2037  fflush(stdout);
2038  }
2039 
2040  // Save out interim basis set
2041  if(saveall) {
2042  static size_t isave=0;
2043  std::ostringstream oss;
2044  oss << "reduce_" << isave << ".gbs";
2045  form_basis(cpl).save_gaussian94(oss.str());
2046  isave++;
2047  }
2048 
2049  // Update references
2050  update_reference(cpl);
2051  }
2052 
2053  return minmog;
2054  }
2055 
2056  void reduce_basis(const std::vector<coprof_t> & cbscpl, std::vector<coprof_t> & cpl, int nred=3, bool docontr=true, bool restr=true, double nelcutoff=0.01, double Porth=true, double saveall=false, double tol=0.0, bool allelectron=true, bool forcesub=false) {
2057  // Reference value
2058  ValueType curval;
2059  {
2060  std::vector< std::vector<coprof_t> > hlp(1,cbscpl);
2061  std::vector<ValueType> hlpvals(compute_values(hlp));
2062  if(hlpvals.size()!=1) {
2063  std::ostringstream oss;
2064  oss << "Error - requested computation of " << 1 << " values but got only " << hlpvals.size() << "!\n";
2065  throw std::runtime_error(oss.str());
2066  }
2067  curval=hlpvals[0];
2068  }
2069  const ValueType cbsval(curval);
2070 
2071  print_value(curval,"CBS limit value");
2072  print_limits(cbscpl,"CBS limit basis");
2073  fflush(stdout);
2074 
2075  // Amount of polarization functions
2076  int npol=maxam(cpl)-atom_am();
2077 
2078  // Current value
2079  {
2080  std::vector< std::vector<coprof_t> > hlp(1,cpl);
2081  std::vector<ValueType> hlpvals(compute_values(hlp));
2082  if(hlpvals.size()!=1) {
2083  std::ostringstream oss;
2084  oss << "Error - requested computation of " << 1 << " values but got only " << hlpvals.size() << "!\n";
2085  throw std::runtime_error(oss.str());
2086  }
2087  curval=hlpvals[0];
2088  }
2089 
2090  print_value(curval,"Starting point value");
2091  print_limits(cpl,"Starting point basis");
2092  fflush(stdout);
2093 
2094  double tau=0.0;
2095  while((tol>0.0 && tau<tol) || (tol==0.0 && npol>=1)) {
2096  // Reduce the profile
2097  tau=reduce_profile(cpl,curval,cbsval,tol,nred,saveall,allelectron,forcesub);
2098 
2099  if(tol==0.0) {
2100  printf("Final composition for %i polarization shells (mog = %e):\n",npol,tau);
2101  print_value(curval,"Current value");
2102  print_limits(cpl,"Current limits");
2103  fflush(stdout);
2104 
2105  // Save basis set
2106  {
2107  char fname[180];
2108  sprintf(fname,"un-co-%i.gbs",npol);
2109  form_basis(cpl).save_gaussian94(fname);
2110 
2111  sprintf(fname,"reduced-%i.dat",npol);
2112  save_limits(cpl,fname);
2113  }
2114  } else {
2115  printf("Final composition for tol = %e, mog = %e:\n",tol,tau);
2116  print_value(curval,"Current value");
2117  print_limits(cpl,"Current limits");
2118  fflush(stdout);
2119 
2120  // Save basis set
2121  {
2122  char fname[180];
2123  sprintf(fname,"un-co-%e.gbs",tol);
2124  form_basis(cpl).save_gaussian94(fname);
2125 
2126  sprintf(fname,"reduced-%e.dat",tol);
2127  save_limits(cpl,fname);
2128  }
2129  }
2130 
2131  // Contract the basis
2132  if(docontr) {
2133  // Use CBS value as reference for contraction
2134  // ValueType contrref(cbsval);
2135 
2136  // Use current value as reference for contraction
2137  ValueType contrref(curval);
2138 
2139  // Threshold to use
2140  double thr;
2141  if(tol==0.0)
2142  thr=tau;
2143  else
2144  thr=tol;
2145 
2146  // Contract the basis, compute mog possibly using P-orthogonalization
2147  BasisSetLibrary contrbas=contract_basis(cpl,contrref,thr,nelcutoff,Porth,restr);
2148 
2149  char fname[180];
2150  if(tol==0.0) {
2151  // General contractions
2152  sprintf(fname,"co-general-%i.gbs",npol);
2153  contrbas.save_gaussian94(fname);
2154 
2155  // Segmented contractions
2156  contrbas.P_orthogonalize();
2157  sprintf(fname,"co-segmented-%i.gbs",npol);
2158  contrbas.save_gaussian94(fname);
2159  } else {
2160  // General contractions
2161  sprintf(fname,"co-general-%e.gbs",tol);
2162  contrbas.save_gaussian94(fname);
2163 
2164  // Segmented contractions
2165  contrbas.P_orthogonalize();
2166  sprintf(fname,"co-segmented-%e.gbs",tol);
2167  contrbas.save_gaussian94(fname);
2168  }
2169  }
2170 
2171  if(tol==0.0) {
2172  // Erase polarization shell
2173  int delam=maxam(cpl);
2174  cpl[delam].start=0.0;
2175  cpl[delam].end=0.0;
2176  cpl[delam].exps.clear();
2177  npol--;
2178  }
2179  }
2180  }
2181 
2183  virtual std::vector<size_t> update_contraction(const std::vector<coprof_t> & cpl, double cutoff)=0;
2184 
2186  BasisSetLibrary contract_basis(const std::vector<coprof_t> & cpl, const ValueType & refval, double tol, double nelcutoff, bool Porth=true, bool restr=true) {
2187 
2188  printf("\n\n%s\n",print_bar("BASIS CONTRACTION").c_str());
2189  fflush(stdout);
2190 
2191  // Update the coefficients, get amount of electrons in each am shell
2192  std::vector<size_t> nel=update_contraction(cpl,nelcutoff);
2193 
2194  Timer ttot;
2195 
2196  // How many functions to contract, counting from the tightest functions.
2197  // Initialize with amount of contractions on each shell; this
2198  // doesn't yet affect the variational freedom of the basis set.
2199  std::vector<size_t> contract(nel);
2200 
2201  // Do the contraction.
2202  while(true) {
2203  Timer tc;
2204 
2205  // Trial contractions and basis sets
2206  std::vector< std::vector<size_t> > trials;
2207  std::vector<BasisSetLibrary> trbas;
2208  std::vector<int> tram;
2209  for(size_t am=0;am<contract.size();am++) {
2210  // Check if there are free functions left
2211  bool free=true;
2212 
2213  // No free functions available.
2214  if(contract[am] >= cpl[am].exps.size())
2215  free=false;
2216 
2217  if(restr && am+1<cpl.size()) {
2218  // Check that the shell has more free functions than the one above.
2219  int Ncur=cpl[am].exps.size()-contract[am];
2220  int Nnext=cpl[am+1].exps.size();
2221  if(am+1<contract.size())
2222  Nnext-=contract[am+1];
2223 
2224  // If the amount of functions is the same or smaller, don't try to contract.
2225  if(Ncur<=Nnext) {
2226  printf("%c shell (%i free funcs) limited due to %c shell (%i free funcs).\n",shell_types[am],Ncur,shell_types[am+1],Nnext);
2227  fflush(stdout);
2228  free=false;
2229  }
2230  }
2231 
2232  if(free) {
2233  // We still have free functions left, form trial contraction
2234  std::vector<size_t> tr(contract);
2235  tr[am]++;
2236 
2237  trials.push_back(tr);
2238  trbas.push_back(form_basis(cpl,tr,Porth));
2239  tram.push_back(am);
2240  }
2241  }
2242 
2243  if(!trbas.size()) {
2244  printf("No functions left to contract.\n");
2245  fflush(stdout);
2246  break;
2247  }
2248 
2249  // Compute values
2250  std::vector<ValueType> trvals=compute_values(trbas);
2251  if(trvals.size()!=trbas.size()) {
2252  std::ostringstream oss;
2253  oss << "Error - requested computation of " << trbas.size() << " values but got only " << trvals.size() << "!\n";
2254  throw std::runtime_error(oss.str());
2255  }
2256 
2257  // and mogs
2258  arma::vec trmogs(trvals.size());
2259  for(size_t i=0;i<trvals.size();i++) {
2260  trmogs(i)=compute_mog(trvals[i],refval,DBL_MAX);
2261  print_scheme(trbas[i],40);
2262  printf("%e\n",trmogs(i));
2263  }
2264  fflush(stdout);
2265 
2266  // Determine minimal mog
2267  arma::uword minind;
2268  double minmog=trmogs.min(minind);
2269 
2270  // Accept move?
2271  if(minmog<=tol) {
2272  printf("Contracting a %c function, mog is %e (%s).\n\n",shell_types[tram[minind]],minmog,tc.elapsed().c_str());
2273  fflush(stdout);
2274  contract=trials[minind];
2275  } else {
2276  // Converged
2277  printf("Minimal mog is %e, converged (%s). Scheme is ",minmog,tc.elapsed().c_str());
2278  fflush(stdout);
2279  // Break while loop
2280  break;
2281  }
2282  }
2283 
2284  // Compile library
2285  BasisSetLibrary ret=form_basis(cpl,contract,false);
2286 
2287  print_scheme(ret);
2288  printf(".\nContraction took %s.\n\n",ttot.elapsed().c_str());
2289  fflush(stdout);
2290 
2291  return ret;
2292  }
2293 };
2294 
2295 #endif
double check_polarization(std::vector< coprof_t > &cpl, ValueType &curval, int am_max, double minpol, double maxpol, double dpol, bool polinterp, double cotol, int nx, bool forcesub=false)
Definition: co-opt.h:749
void save_gaussian94(const std::string &filename, bool append=false) const
Save basis set to file in Gaussian&#39;94 format.
Definition: basislibrary.cpp:1623
double scan_profile(std::vector< coprof_t > &cpl, ValueType &curval, int npoints, double dpol, double tol, bool allam=true, bool bothsd=true, bool forcesub=false)
Definition: co-opt.h:1075
A timer routine.
Definition: timer.h:43
bool comply_subset(std::vector< coprof_t > &cpl, ValueType &curval)
Reform to comply to subset requirement.
Definition: co-opt.h:981
arma::vec lga
Logarithms of scanning exponents .
Definition: completeness_profile.h:37
arma::vec span_width(int am, double tol, double &width, int nx)
Span necessary width.
Definition: co-opt.h:251
BasisSetLibrary contract_basis(const std::vector< coprof_t > &cpl, const ValueType &refval, double tol, double nelcutoff, bool Porth=true, bool restr=true)
Contract basis set. nel holds amount of states of each angular momentum. P toggles intermediate P-ort...
Definition: co-opt.h:2186
std::vector< coprof_t > initial_profile(const std::vector< size_t > &nfuncs, double cotol=1e-4)
Definition: co-opt.h:497
Completeness profiles for a given element.
Definition: completeness_profile.h:35
virtual std::vector< size_t > update_contraction(const std::vector< coprof_t > &cpl, double cutoff)=0
Update the contraction coefficients, return the amount of electrons in each am shell.
Helper for maxwidth_exps_table.
Definition: co-opt.h:102
Structure for completeness-optimized basis set.
Definition: co-opt.h:71
std::string elapsed() const
Get pretty-printed elapsed time.
Definition: timer.cpp:102
double tighten_profile(std::vector< coprof_t > &cpl, ValueType &curval, double tau, int nxadd=1, bool forcesub=false)
Tighten the shells until mog &lt; tau. Returns maximal mog.
Definition: co-opt.h:1482
int get_nfull() const
Get amount of fully optimized exponents.
Definition: co-opt.h:214
double end
End of completeness profile.
Definition: co-opt.h:75
virtual void update_reference(const std::vector< coprof_t > &cpl)
Update reference wave function.
Definition: co-opt.h:140
arma::vec optimize_completeness(int am, double min, double max, int Nf, double &mog)
Get exponents.
Definition: co-opt.h:224
void add_function(FunctionShell f)
Add a shell of functions to the basis.
Definition: basislibrary.cpp:317
int n_tau
Moment of profile to optimize.
Definition: co-opt.h:121
double start
Start of completeness profile.
Definition: co-opt.h:73
arma::vec maxwidth_exps_table(int am, double tol, size_t nexp, double &width)
Get completeness-optimized exponents.
Definition: co-opt.h:229
int get_ntau() const
Get moment of completeness profile to optimize.
Definition: co-opt.h:204
virtual int atom_am() const =0
Get SCF free-atom angular momentum.
std::vector< coprof_t > initial_profile(const std::vector< int > &nshells, int nexp, double cotol=1e-4)
Definition: co-opt.h:486
void find_cbs_limit(std::vector< coprof_t > &cpl, ValueType &curval, double cotol, double minpol, double maxpol, double dpol, bool extend=true, int next=3, bool scan=true, int nscan=5, bool polinterp=true, int nxpol=1, bool doadd=true, int nxext=1, int am_max=max_am, bool cbsinterp=true, double cbsthr=0.0, double delta=0.9, bool allam=true, bool bothsd=true, bool forcesub=false)
Extend the basis set till the CBS limit.
Definition: co-opt.h:1588
void mog_decompose(const std::vector< int > &tram, const arma::vec &mogs, std::vector< arma::uvec > &amidx, std::vector< arma::vec > &ammog) const
Decompose vectorized mogs per angular momentum.
Definition: co-opt.h:146
Basis set library class.
Definition: basislibrary.h:231
CompletenessOptimizer()
Constructor.
Definition: co-opt.h:191
System default location for basis sets.
Definition: basislibrary.h:55
size_t politer
Iteration number (for saving out polarization and stability scans)
Definition: co-opt.h:126
double tol
Tolerance.
Definition: co-opt.h:77
virtual double compute_energy(const std::vector< coprof_t > &cpl)=0
Compute the energy.
double extend_profile(std::vector< coprof_t > &cpl, ValueType &curval, double tau, int next=3, int nxadd=1, bool forcesub=false)
Extend the current shells until mog &lt; tau. Returns maximal mog.
Definition: co-opt.h:1347
virtual std::vector< ValueType > compute_values(const std::vector< BasisSetLibrary > &baslib)=0
Compute values with given basis sets (implementation can be parallellized)
virtual void print_value(const ValueType &value, std::string msg)=0
Print value.
Basis set for an element.
Definition: basislibrary.h:102
void decontract()
Decontract set.
Definition: basislibrary.cpp:467
virtual ~CompletenessOptimizer()
Destructor.
Definition: co-opt.h:200
int n_full
Amount of fully optimized exponents on each end.
Definition: co-opt.h:123
void set_nfull(int n)
Set amount of fully optimized exponents.
Definition: co-opt.h:219
virtual BasisSetLibrary form_basis(const std::vector< coprof_t > &cpl, std::vector< size_t > contract, bool Porth=true)
Form a contracted basis set. Porth toggles P-orthogonalization.
Definition: co-opt.h:132
void add_element(const ElementBasisSet &el)
Add element to basis set.
Definition: basislibrary.cpp:2084
void save_limits(const std::vector< coprof_t > &cpl, const std::string &fname) const
Save the limits.
Definition: co-opt.h:284
std::vector< coprof_t > initial_profile(const std::vector< size_t > &nfuncs, std::vector< double > offsets, double cotol=1e-4)
Definition: co-opt.h:507
std::vector< coprof_t > load_limits(const BasisSetLibrary &baslib, double tol, int maxam)
Generate limits from existing basis set.
Definition: co-opt.h:365
Completeness optimizer class.
Definition: co-opt.h:118
virtual BasisSetLibrary form_basis(const std::vector< coprof_t > &cpl)
Form a basis set.
Definition: co-opt.h:275
arma::vec exps
The completeness-optimized exponents.
Definition: co-opt.h:80
void P_orthogonalize(double cutoff=1e-8, double Cortho=1e-4)
Definition: basislibrary.cpp:2196
void set_ntau(int n)
Set moment of completeness profile to optimize.
Definition: co-opt.h:209
std::vector< ElementBasisSet > get_elements() const
Get elements.
Definition: basislibrary.cpp:2110
std::vector< ValueType > compute_values(const std::vector< std::vector< coprof_t > > &cpl)
Compute value with given basis set.
Definition: co-opt.h:726
std::vector< compprof_am_t > shells
Completeness profiles for all angular momenta.
Definition: completeness_profile.h:39
void set()
Zero timer.
Definition: timer.cpp:62
double reduce_profile(std::vector< coprof_t > &cpl, ValueType &curval, const ValueType &refval, double tol=0.0, int nred=3, bool saveall=false, bool allelectron=true, bool forcesub=false)
Reduce the profile until only a single function is left on the highest am shell (polarization / corre...
Definition: co-opt.h:1835
void add_exponent(double C, double z)
Add exponent into contraction.
Definition: basislibrary.cpp:214
std::vector< coprof_t > load_limits(const std::string &fname)
Load limits.
Definition: co-opt.h:296