1 /* 2 * This source code is part of 3 * 4 * E R K A L E 5 * - 6 * HF/DFT from Hel 7 * 8 * Written by Susi Lehtola, 2010-2011 9 * Copyright (c) 2010-2011, Susi Lehtola 10 * 11 * This program is free software; you can redistribute it and/or 12 * modify it under the terms of the GNU General Public License 13 * as published by the Free Software Foundation; either version 2 14 * of the License, or (at your option) any later version. 15 */ 16 17 #ifndef ERKALE_OPTCOMP 18 #define ERKALE_OPTCOMP 19 20 /// Minimum allowed value of deviation from completeness (for numerical stability) 21 #define MINTAU pow(10.0,-5.5) 22 23 #include "../global.h" 24 #include <armadillo> 25 #include <vector> 26 27 extern "C" { 28 #include <gsl/gsl_vector.h> 29 } 30 31 /// Parameters for completeness scan 32 typedef struct { 33 /// Angular momentum of shell to optimize 34 int am; 35 /// Which moment to optimize 36 int n; 37 38 /// Scanning exponents to optimize against 39 arma::vec scanexp; 40 41 /// Add in exponent at the center? 42 bool odd; 43 /// Amount of even-tempered exponents near the center 44 size_t neven; 45 /// Amount of fully optimized exponents at the edges 46 size_t nfull; 47 } completeness_scan_t; 48 49 /// Get exponents. x contains the natural logarithms 50 arma::vec get_exponents(const gsl_vector *x, const completeness_scan_t & p); 51 52 /// Compute self-overlap \f$ S_{ij} \f$ 53 arma::mat self_overlap(const arma::vec & z, int am); 54 55 /// Compute completeness profile 56 arma::vec completeness_profile(const gsl_vector * x, void * params); 57 58 /// Evaluate measure of goodness 59 double compl_mog(const gsl_vector * x, void * params); 60 61 /// Evaluate gradient of measure of goodness 62 void compl_mog_df(const gsl_vector * x, void * params, gsl_vector * g); 63 64 /// Evaluate gradient and measure of goodness 65 void compl_mog_fdf(const gsl_vector * x, void * params, double * f, gsl_vector * df); 66 67 /** 68 * Optimize completeness profile for angular momentum am in exponent 69 * range from 10^{min} to 10^{max}. n gives the moment to optimize for 70 * (1 for maximal area, 2 for minimal rms deviation from unity). 71 * 72 * By default, only the 4 outermost exponents on each end are fully 73 * optimized, while the inner exponents are well described by an 74 * even-tempered formula. If nfull>=Nf, a full optimization is performed. 75 */ 76 arma::vec optimize_completeness(int am, double min, double max, int Nf, int n=1, bool verbose=true, double *mog=NULL, int nfull=4); 77 78 /** 79 * Optimize completeness profile for angular momentum am in exponent 80 * range from 10^{min} to 10^{max}. n gives the moment to optimize for 81 * (1 for maximal area, 2 for minimal rms deviation from unity). 82 * 83 * This routine uses Fletcher-Reeves conjugate gradients. 84 * 85 * By default, only the 4 outermost exponents on each end are fully 86 * optimized, while the inner exponents are well described by an 87 * even-tempered formula. If nfull>=Nf, a full optimization is performed. 88 */ 89 arma::vec optimize_completeness_cg(int am, double min, double max, int Nf, int n=1, bool verbose=true, double *mog=NULL, int nfull=4); 90 91 /** 92 * Optimize completeness profile for angular momentum am in exponent 93 * range from 10^{min} to 10^{max}. n gives the moment to optimize for 94 * (1 for maximal area, 2 for minimal rms deviation from unity). 95 * 96 * This is the old version of the routine, which uses the Nead-Miller algorithm. 97 * 98 * By default, only the 4 outermost exponents on each end are fully 99 * optimized, while the inner exponents are well described by an 100 * even-tempered formula. If nfull>=Nf, a full optimization is 101 * performed. 102 */ 103 arma::vec optimize_completeness_simplex(int am, double min, double max, int Nf, int n=1, bool verbose=true, double *mog=NULL, int nfull=4); 104 105 /// Calculate maximum width to obtain tolerance with given amount of exponents 106 double maxwidth(int am, double tol, int nexp, int n=1, int nfull=4); 107 108 /// Calculate exponents corresponding to maximum width to obtain tolerance with given amount of exponents 109 arma::vec maxwidth_exps(int am, double tol, int nexp, double & width, int n=1, int nfull=4); 110 111 /// Move exponents in logarithmic scale to start at x instead of 0.0 112 arma::vec move_exps(const arma::vec & exps, double start); 113 114 /// Get exponents that span [start,end] with tau <= tol 115 arma::vec get_exponents(int am, double start, double end, double tol, int n=1, bool verbose=false, int nfull=4); 116 117 #endif 118