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"
58 #include "../stringutil.h"
59 #include "../unitary.h"
60 #include "../elements.h"
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);
91 std::vector<coprof_t> augment_basis(
const std::vector<coprof_t> & cplo,
int ndiffuse,
int ntight);
94 ElementBasisSet get_element_library(
const std::string & el,
const std::vector<coprof_t> & cpl);
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);
117 template<
typename ValueType>
129 virtual double compute_energy(
const std::vector<coprof_t> & cpl)=0;
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());
161 int maxam=arma::conv_to<arma::ivec>::from(tram).max();
164 amidx.resize(maxam+1);
166 ammog.resize(maxam+1);
169 std::vector< std::vector<double> > amm(maxam+1);
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));
178 for(
int am=0;am<=maxam;am++) {
180 amidx[am]=arma::conv_to<arma::uvec>::from(ami[am]);
181 ammog[am]=arma::conv_to<arma::vec>::from(amm[am]);
225 return ::optimize_completeness(am, min, max, Nf,
n_tau,
false, &mog,
n_full);
231 static std::vector<co_exps_t> opt(max_am+1);
234 if(opt[am].tol!=tol || opt[am].exps.size()!=nexp) {
236 opt[am].exps=maxwidth_exps(am,tol,nexp,opt[am].w,
n_tau,
n_full);
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());
251 arma::vec
span_width(
int am,
double tol,
double & width,
int nx) {
260 for(nx++;nx<=NFMAX;nx++) {
277 for(
int Z=1;Z<=maxZ;Z++)
278 baslib.
add_element(get_element_library(element_symbols[Z],cpl));
284 void save_limits(
const std::vector<coprof_t> & cpl,
const std::string & fname)
const {
285 FILE *out=fopen(fname.c_str(),
"w");
287 throw std::runtime_error(
"Error opening completeness range output file.\n");
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());
297 FILE *in=fopen(fname.c_str(),
"r");
299 throw std::runtime_error(
"Error opening completeness range file.\n");
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());
311 std::vector<coprof_t> cpl(max+1);
314 for(
int am=0;am<=max;am++) {
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());
328 std::ostringstream oss;
329 oss <<
"Read in am " << amt <<
" does not match supposed am " << am <<
"!\n";
330 throw std::runtime_error(oss.str());
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";
341 printf(
"Warning: %s",oss.str().c_str());
345 double width=cpl[am].end-cpl[am].start;
353 cpl[am].start-=dw/2.0;
357 cpl[am].exps=move_exps(exps,cpl[am].start);
367 double cplthr=1.0-tol;
370 std::vector<coprof_t> cpl(max_am+1);
371 for(
size_t i=0;i<cpl.size();i++) {
381 for(
size_t iel=0;iel<els.size();iel++) {
391 for(
size_t am=0;am<prof.
shells.size();am++) {
394 for(
size_t ia=0;ia<prof.
lga.size();ia++)
395 if(prof.
shells[am].Y[ia]>=cplthr) {
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) {
408 if(cpl[am].start>low)
416 for(
int am=maxam;am<max_am;am++) {
419 cpl[am].exps.clear();
422 for(
int am=0;am<std::min(maxam,max_am);am++)
423 if(cpl[am].start==DBL_MAX && cpl[am].end==-DBL_MAX) {
428 printf(
"Initial composition:\n");
431 for(
int am=0;am<max_am;am++) {
432 if(cpl[am].start == 0.0 && cpl[am].end == 0.0)
436 double width=cpl[am].end-cpl[am].start;
439 std::vector<double> widths;
440 std::vector< std::vector<double> > explist;
442 widths.push_back(0.0);
445 for(
int nf=1;nf<=NFMAX;nf++) {
450 explist.push_back(exps);
460 for(
size_t i=0;i<widths.size();i++)
461 if(fabs(widths[i]-width)<mind) {
463 mind=fabs(widths[i]-width);
467 cpl[am].start-=(widths[nf]-width)/2.0;
468 cpl[am].end+=(widths[nf]-width)/2.0;
470 cpl[am].exps=move_exps(explist[nf],cpl[am].start);
472 printf(
"%c % .2f % .2f %2i\n",shell_types[am],cpl[am].start,cpl[am].end,(
int) cpl[am].exps.size());
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++)
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);
507 std::vector<coprof_t>
initial_profile(
const std::vector<size_t> & nfuncs, std::vector<double> offsets,
double cotol=1e-4) {
511 std::vector<coprof_t> cpl(nfuncs.size());
512 for(
size_t i=0;i<cpl.size();i++) {
519 std::vector<arma::vec> exps(cpl.size());
520 arma::vec widths(cpl.size());
523 printf(
"\nStarting composition:\n");
524 for(
size_t am=0;am<cpl.size();am++) {
526 cpl[am].start=offsets[am];
527 cpl[am].end=cpl[am].start+widths[am];
528 cpl[am].exps=exps[am];
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());
535 arma::vec relw(widths);
536 for(
size_t am=0;am<relw.size();am++)
537 relw(am)/=exps[am].size();
540 const double h0=1e-3;
550 bool cg=iiter%cpl.size()!=0;
552 printf(
"\nConjugate gradient iteration %i\n",(
int) iiter+1);
554 printf(
"\nSteepest descent iteration %i\n",(
int) iiter+1);
558 arma::vec g(cpl.size());
559 for(
size_t am=0;am<cpl.size();am++) {
561 double h=h0*relw(am);
563 std::vector<coprof_t> lcpl(cpl);
566 lcpl[am].exps=move_exps(exps[am],lcpl[am].start);
568 std::vector<coprof_t> rcpl(cpl);
571 rcpl[am].exps=move_exps(exps[am],rcpl[am].start);
577 bool Elok=(std::isfinite(El) && fabs(El)!=DBL_MAX);
578 bool Erok=(std::isfinite(Er) && fabs(Er)!=DBL_MAX);
581 g(am)=(Er-El)/(2*h) * relw(am);
585 printf(
"\t-grad(%c) = % e\n", shell_types[am], -g(am));
593 memset(minus,
'-',Nm);
595 printf(
"\t%s\n",minus);
597 printf(
"\tgrad norm = % e\n",arma::norm(g,2));
605 g-=arma::dot(g,gold)*gold;
624 for(
int i=1;i<Ntr;i++)
625 h(i)=hs*std::pow(2.0,i-1);
627 printf(
"\n\tLine search:\n");
628 printf(
"\t%5s %12s %13s %13s\n",
"trial",
"step",
"E",
"dE");
632 for(itr=0;itr<Ntr;itr++) {
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);
645 printf(
"\t%2i/%-2i %e % e\n",(
int) (itr+1), (
int) Ntr, h(itr), E(itr));
647 printf(
"\t%2i/%-2i %e % e % e\n",(
int) (itr+1), (
int) Ntr, h(itr), E(itr), E(itr)-E(0));
651 if(itr>0 && E(itr)>E(itr-1))
661 printf(
"\tMinimum with step %e, decrease by % e.\n",h(ind),E(ind)-E(0));
667 arma::vec x=h.subvec(ind-1,ind+1);
668 arma::vec y=E.subvec(ind-1,ind+1);
671 arma::vec c=fit_polynomial(x,y);
673 arma::vec dc=derivative_coefficients(c);
675 double st=smallest_positive(solve_roots(dc));
677 printf(
"\tInterpolation gives step %e, ",st);
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);
693 printf(
"further decrease by % e, accepted.\n",Etr-E.min());
695 printf(
"further increase by % e, rejected.\n",Etr-E.min());
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);
706 printf(
"Iteration took %s.\n",tstep.
elapsed().c_str());
707 print_limits(cpl,
"\nCurrent limits");
716 print_limits(cpl,
"Initial composition");
717 printf(
"Generation of initial profile took %s.\n",t.
elapsed().c_str());
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++)
734 virtual double compute_mog(
const ValueType & val,
const ValueType & ref,
double failval)=0;
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) {
751 printf(
"\n\n%s\n",print_bar(
"POLARIZATION SHELLS").c_str());
756 const int max=maxam(cpl);
758 printf(
"Won't add another polarization shell, MaxAM=%i reached.\n",am_max);
769 const int addam=max+1;
773 std::ostringstream oss;
774 oss <<
"Invalid value " << dpol <<
" for DPol.\n";
775 throw std::runtime_error(oss.str());
786 sp=(wnext-wcur)/dpol;
789 double ww=maxpol-minpol;
791 size_t Np=1+round(ww/sp);
797 double maxp=maxpol+(width-ww);
800 arma::vec startp=arma::linspace(minp,maxp,Np);
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);
807 if(cpl.size() <= (size_t) addam) {
810 cpl[addam].tol=cotol;
814 std::vector< std::vector<coprof_t> > cpls(startp.n_elem,cpl);
815 for(
size_t iexp=0;iexp<startp.n_elem;iexp++) {
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);
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());
831 printf(
"%11s %8s %12s\n",
"trial ",
"exponent",
"mog");
833 arma::vec mogs(startp.n_elem);
834 for(
size_t i=0;i<vals.size();i++) {
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));
848 std::ostringstream oss;
849 oss <<
"polmog_" <<
politer <<
"_" << shell_types[addam] <<
".dat";
853 arma::mat savemog(mogs.n_rows,2);
854 savemog.col(0)=startp;
856 savemog.save(fname,arma::raw_ascii);
859 for(
int am=0;am<=maxam(cpl);am++) {
868 lim(0,1)=cpl[am].start;
869 lim(1,1)=cpl[am].start;
872 lim(0,2)=cpl[am].end;
873 lim(1,2)=cpl[am].end;
876 std::ostringstream oss;
877 oss <<
"limits_" <<
politer <<
"_" << shell_types[am] <<
".dat";
878 lim.save(oss.str(),arma::raw_ascii);
883 std::ostringstream oss;
884 oss <<
"polbas_" <<
politer <<
".gbs";
891 double maxmog=mogs.max(imax);
894 printf(
"%11s % 7.5f %e\n",
"max mog",startp(imax),maxmog);
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);
903 arma::uword iintmean;
905 intmean=moginterp.max(iintmean);
908 printf(
"%11s % 7.5f %e\n",
"interp mog",finergrid(iintmean),intmean);
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);
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());
932 intmog=compute_mog(intval,curval,0.0);
938 printf(
"%11s %8s %e\n",
"real value",
"",intmog);
941 if(intmog > maxmog) {
942 printf(
"Interpolation accepted.\n");
948 printf(
"Interpolation rejected.\n");
957 printf(
"Interpolation wouldn't result in a better value.\n");
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());
988 for(
int am=maxam(cpl);am>=0;am--) {
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;
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;
1008 double tmove=smove+emove;
1009 if(tmove==0.0 && Yam.
tol==cpl[am].tol)
1014 printf(
"\n\n%s\n",print_bar(
"Subset check",
'-').c_str());
1020 double curw=cpl[am].end-cpl[am].start;
1022 double nw=curw+tmove;
1029 int nadd=(int) (exps.size()-cpl[am].exps.size());
1030 double sreal=(realw-curw)*smove/tmove;
1031 double ereal=(realw-curw)*emove/tmove;
1034 printf(
"%c lower limit should be moved by % .3f; moved by % .3f.\n",shell_types[am],smove,sreal);
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());
1051 std::vector< std::vector<coprof_t> > cplv(1,cpl);
1055 print_limits(cpl,
"Updated limits");
1057 printf(
"%s\n\n",print_bar(
"End check",
'-').c_str());
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) {
1077 printf(
"\n\n%s\n",print_bar(
"SHELL STABILITY").c_str());
1086 printf(
"Scanning for shell stability.\n");
1087 printf(
"Using %.3f points per exponent interval.\n\n",dpol);
1091 std::vector<double> trexp;
1092 std::vector<int> tram;
1095 arma::vec spacing(maxam(cpl)+1);
1096 for(
int scanam=0;scanam<=maxam(cpl);scanam++) {
1107 arma::uword nexp=cpl[scanam].exps.n_elem+1;
1116 spacing(scanam)=step;
1119 printf(
"Spacing between points on %c shell is %.5f\n",shell_types[scanam],step);
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);
1129 for(
int n=nmin;n<=nmax;n++) {
1130 trexp.push_back(cpl[scanam].end+n*step/dpol);
1131 tram.push_back(scanam);
1136 std::vector< BasisSetLibrary > trbas;
1137 for(
size_t iexp=0;iexp<trexp.size();iexp++) {
1141 for(
size_t iel=0;iel<els.size();iel++) {
1151 trbas.push_back(baslib);
1156 printf(
"No shells to scan stability of.\n");
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());
1169 printf(
"%11s %2s %8s %12s\n",
"trial ",
"am",
"exponent",
"mog");
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));
1179 std::vector<arma::uvec> amidx;
1180 std::vector<arma::vec> 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) {
1187 double maxmog=ammog[am].max(maxind);
1188 printf(
"\t%-2c % 7.5f %e\n",shell_types[am],trexp[amidx[am](maxind)],maxmog);
1190 printf(
"\t%-2c %8s %e\n",shell_types[am],
"",0.0);
1195 for(
int am=0;am<=maxam(cpl);am++) {
1196 std::ostringstream fname;
1197 fname <<
"scanmog_" <<
politer <<
"_" << shell_types[am] <<
".dat";
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]];
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);
1213 double maxmog=mogs.max(imax);
1216 printf(
"%11s %-2c % 7.5f %e\n",
"max mog",shell_types[tram[imax]],trexp[imax],maxmog);
1225 for(
int scanam=0;scanam<=maxam(cpl);scanam++) {
1228 if(!allam && ((arma::uword) scanam != imax))
1232 if(!ammog[scanam].n_elem)
1236 std::vector<arma::uword> idxlist;
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]);
1248 stidx.push_back(itr);
1249 stmog.push_back(mogs[itr]);
1254 arma::vec stm(arma::conv_to<arma::vec>::from(stmog));
1255 arma::vec dfm(arma::conv_to<arma::vec>::from(dfmog));
1258 arma::uword stind, dfind;
1259 double stmax=stm.max(stind);
1260 double dfmax=dfm.max(dfind);
1262 idxlist.push_back(stidx[stind]);
1264 idxlist.push_back(dfidx[dfind]);
1268 double ammax=ammog[scanam].max(amind);
1272 idxlist.push_back(amidx[scanam](amind));
1276 for(
size_t i=0;i<idxlist.size();i++) {
1284 for(
int am=minam;am<=maxam;am++) {
1287 if(trexp[imax] < cpl[am].start) {
1289 double curw=cpl[am].end-cpl[am].start;
1291 double nw=cpl[am].end-trexp[imax];
1294 arma::vec exps=
span_width(am,cpl[am].tol,realw,cpl[am].exps.size());
1297 cpl[am].start-=realw-curw;
1298 cpl[am].exps=move_exps(exps,cpl[am].start);
1299 moved=-(realw-curw);
1301 }
else if(trexp[imax] > cpl[am].end) {
1303 double curw=cpl[am].end-cpl[am].start;
1305 double nw=trexp[imax]-cpl[am].start;
1308 arma::vec exps=
span_width(am,cpl[am].tol,realw,cpl[am].exps.size());
1311 cpl[am].end+=realw-curw;
1312 cpl[am].exps=move_exps(exps,cpl[am].start);
1313 moved=+(realw-curw);
1315 }
else if(am==scanam) {
1316 throw std::runtime_error(
"Possible bug in scan_limits - maximum inside profile!\n");
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());
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());
1330 std::vector< std::vector<coprof_t> > hlp(1,cpl);
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());
1340 printf(
"Stability scan done in %s.\n\n",t.
elapsed().c_str());
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);
1353 throw std::runtime_error(
"Tolerance must be positive!\n");
1355 throw std::runtime_error(
"The NExt parameter must be at least two.\n");
1365 std::vector< std::vector<coprof_t> > trials;
1367 std::vector< std::string > descr;
1369 std::vector< int > tram;
1371 for(
int am=0;am<=maxam(cpl);am++) {
1375 printf(
"Determining exponents for %c shell ... ",shell_types[am]);
1382 step=width-(cpl[am].end-cpl[am].start);
1383 printf(
"step size is %7.5f (%s).\n",step,t.
elapsed().c_str());
1388 for(
int itr=0;itr<next;itr++) {
1390 double d=itr*1.0/(next-1);
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);
1398 trials.push_back(trcpl);
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);
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);
1413 if(!trials.size()) {
1414 printf(
"No shells can be extended.\n");
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());
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);
1432 std::vector<arma::uvec> amidx;
1433 std::vector<arma::vec> 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]));
1440 printf(
"\t%-2c %e\n",shell_types[am],0.0);
1445 double maxmog=mogs.max(maxind);
1449 printf(
"Maximal mog is %e, converged.\n\n",maxmog);
1451 print_limits(cpl,
"Final limits");
1468 curval=trvals[maxind];
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());
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);
1495 std::vector< std::vector<coprof_t> > trials;
1496 std::vector< int > tram;
1497 for(
int am=0;am<=maxam(cpl);am++) {
1501 std::vector<coprof_t> add(cpl);
1504 printf(
"Determining exponents for %c shell ... ",shell_types[am]);
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());
1511 if (add[am].tol>=MINTAU) {
1512 trials.push_back(add);
1518 if(!trials.size()) {
1519 printf(
"No shells can be tightened.\n");
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());
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);
1537 std::vector<arma::uvec> amidx;
1538 std::vector<arma::vec> 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]));
1545 printf(
"\t%-2c %e\n",shell_types[am],0.0);
1550 double maxmog=trmog.max(maxind);
1554 printf(
"Maximal mog is %e, converged.\n\n",maxmog);
1556 print_limits(cpl,
"Final tolerances");
1570 curval=trvals[maxind];
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());
1582 virtual int atom_am()
const=0;
1585 virtual void print_value(
const ValueType & value, std::string msg)=0;
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) {
1590 int npol=maxam(cpl)-
atom_am();
1593 if(!extend && !scan)
1594 throw std::runtime_error(
"You need to check for profile expansion either by scans or extension, or both!\n");
1598 std::vector< std::vector<coprof_t> > hlp(1,cpl);
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());
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);
1618 std::vector<double> polmogs;
1629 extmog=std::max(extmog,amog);
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);
1638 if(extmog >= std::max(polmog,cbsthr)) {
1640 tau=exp((1.0-delta)*log(extmog) + delta*log(std::max(polmog,cbsthr)));
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);
1653 print_limits(scancpl,
"Compound limits");
1660 extmog=std::max(extmog,scanmog);
1662 if(scanmog>=std::max(polmog,cbsthr)) {
1664 double mog=compute_mog(scanval,curval,0.0);
1671 printf(
"\n\nInstability detected with real mog = %e, restarting extension.\n",mog);
1676 printf(
"\n\nInstability detected, but compound calculation failed to converge.\n");
1677 printf(
"Forming new, individual trials.\n");
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()) {
1685 std::vector<coprof_t> tr(cpl);
1688 trials.push_back(tr);
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);
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]);
1706 double maxmog=trmog.max(maxind);
1710 curval=trvals[maxind];
1717 print_limits(cpl,
"Current limits");
1725 if(maxam(polcpl)>maxam(cpl))
1726 polmogs.push_back(polmog);
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);
1736 printf(
"\nFinal composition for %i polarization shells (tau = %e):\n",npol,polmog);
1738 print_limits(cpl,
"Current limits");
1745 sprintf(fname,
"un-co-ref-%i.gbs",npol);
1749 sprintf(fname,
"completeness-%i.dat",npol);
1756 if(std::max(extmog,polmog)<cbsthr) {
1757 printf(
"\nCBS threshold reached: extension mog was %e, polarization mog %e.\n",extmog,polmog);
1760 }
else if(maxam(cpl)==am_max) {
1764 if(cbsinterp && polmogs.size()>=2) {
1766 double hmog=polmogs[polmogs.size()-1];
1768 double lmog=polmogs[polmogs.size()-2];
1784 std::vector<coprof_t> scancpl(cpl);
1785 ValueType scanval(curval);
1786 scanmog=
scan_profile(scancpl,scanval,nscan,dpol,tau,allam,bothsd,forcesub);
1790 double mog=compute_mog(scanval,curval,0.0);
1794 printf(
"\n\nInstability detected with real mog = %e, restarting extension.\n",mog);
1797 print_limits(cpl,
"Current limits");
1814 tau=std::min(tau,std::max(extmog,cbsthr));
1819 print_limits(cpl,
"Final composition:");
1826 sprintf(fname,
"un-co-ref.gbs");
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());
1838 printf(
"Running until would drop last function of %c shell.\n",shell_types[maxam(cpl)]);
1840 printf(
"Running until mog >= %e.\n",tol);
1844 throw std::runtime_error(
"The NRed parameter must be at least two.\n");
1857 std::vector< std::vector<coprof_t> > trials;
1858 std::vector< std::string > descr;
1859 std::vector< int > tram;
1861 for(
int am=0;am<=maxam(cpl);am++) {
1862 if(!cpl[am].exps.size())
1869 for(
int sam=am+1;sam<=maxam(cpl);sam++) {
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]);
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;
1889 trials.push_back(trcpl);
1893 sprintf(msg,
"Dropped last exponent of %c shell",shell_types[am]);
1894 descr.push_back(msg);
1901 printf(
"Determining exponents for %c shell ... ",shell_types[am]);
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);
1912 step=(cpl[am].end-cpl[am].start)-width;
1913 printf(
"step size is %7.5f (%s).\n",step,t.
elapsed().c_str());
1918 for(
int itr=0;itr<nred;itr++) {
1920 double d=itr*1.0/(nred-1);
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);
1930 for(
int ham=am+1;ham<=maxam(cpl);ham++) {
1931 if(cpl[ham].start<trcpl[am].start)
1933 if(cpl[ham].end>trcpl[am].end)
1940 trials.push_back(trcpl);
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);
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);
1957 for(
int ham=am+1;ham<=maxam(cpl);ham++)
1960 if(delcpl[ham].exps.size()>1 && delcpl[am].tol > delcpl[ham].tol)
1963 trials.push_back(delcpl);
1966 sprintf(msg,
"Dropped exponent from %c shell",shell_types[am]);
1967 descr.push_back(msg);
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());
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);
1991 std::vector<arma::uvec> amidx;
1992 std::vector<arma::vec> 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]));
1999 printf(
"\t%-2c %e\n",shell_types[am],0.0);
2004 minmog=trmog.min(minind);
2007 for(
size_t i=0;i<trmog.n_elem;i++)
2008 if(trmog[i] == minmog && tram[i] > tram[minind])
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);
2015 print_limits(cpl,
"Final limits");
2020 }
else if(tol>0.0 && minmog>tol) {
2022 printf(
"Minimal mog is %e for %c shell, converged.\n\n",minmog,shell_types[tram[minind]]);
2024 print_limits(cpl,
"Final limits");
2034 curval=trvals[minind];
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());
2042 static size_t isave=0;
2043 std::ostringstream oss;
2044 oss <<
"reduce_" << isave <<
".gbs";
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) {
2060 std::vector< std::vector<coprof_t> > hlp(1,cbscpl);
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());
2069 const ValueType cbsval(curval);
2072 print_limits(cbscpl,
"CBS limit basis");
2076 int npol=maxam(cpl)-
atom_am();
2080 std::vector< std::vector<coprof_t> > hlp(1,cpl);
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());
2091 print_limits(cpl,
"Starting point basis");
2095 while((tol>0.0 && tau<tol) || (tol==0.0 && npol>=1)) {
2097 tau=
reduce_profile(cpl,curval,cbsval,tol,nred,saveall,allelectron,forcesub);
2100 printf(
"Final composition for %i polarization shells (mog = %e):\n",npol,tau);
2102 print_limits(cpl,
"Current limits");
2108 sprintf(fname,
"un-co-%i.gbs",npol);
2111 sprintf(fname,
"reduced-%i.dat",npol);
2115 printf(
"Final composition for tol = %e, mog = %e:\n",tol,tau);
2117 print_limits(cpl,
"Current limits");
2123 sprintf(fname,
"un-co-%e.gbs",tol);
2126 sprintf(fname,
"reduced-%e.dat",tol);
2137 ValueType contrref(curval);
2152 sprintf(fname,
"co-general-%i.gbs",npol);
2157 sprintf(fname,
"co-segmented-%i.gbs",npol);
2161 sprintf(fname,
"co-general-%e.gbs",tol);
2166 sprintf(fname,
"co-segmented-%e.gbs",tol);
2173 int delam=maxam(cpl);
2174 cpl[delam].start=0.0;
2176 cpl[delam].exps.clear();
2183 virtual std::vector<size_t>
update_contraction(
const std::vector<coprof_t> & cpl,
double cutoff)=0;
2188 printf(
"\n\n%s\n",print_bar(
"BASIS CONTRACTION").c_str());
2199 std::vector<size_t> contract(nel);
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++) {
2214 if(contract[am] >= cpl[am].exps.size())
2217 if(restr && am+1<cpl.size()) {
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];
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);
2234 std::vector<size_t> tr(contract);
2237 trials.push_back(tr);
2244 printf(
"No functions left to contract.\n");
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());
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));
2268 double minmog=trmogs.min(minind);
2272 printf(
"Contracting a %c function, mog is %e (%s).\n\n",shell_types[tram[minind]],minmog,tc.
elapsed().c_str());
2274 contract=trials[minind];
2277 printf(
"Minimal mog is %e, converged (%s). Scheme is ",minmog,tc.
elapsed().c_str());
2288 printf(
".\nContraction took %s.\n\n",ttot.
elapsed().c_str());
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'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 < 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 < 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