1 ////////////////////////////////////////////////////////////////////////////////////// 2 // This file is distributed under the University of Illinois/NCSA Open Source License. 3 // See LICENSE file in top directory for details. 4 // 5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. 6 // 7 // File developed by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign 8 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory 9 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign 10 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 11 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory 12 // 13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 14 ////////////////////////////////////////////////////////////////////////////////////// 15 16 17 #ifndef QMCPLUSPLUS_COSTFUNCTIONBASE_H 18 #define QMCPLUSPLUS_COSTFUNCTIONBASE_H 19 20 #include <deque> 21 #include <set> 22 #include "Configuration.h" 23 #include "Optimize/OptimizeBase.h" 24 #include "QMCHamiltonians/QMCHamiltonian.h" 25 #include "QMCWaveFunctions/TrialWaveFunction.h" 26 #include "Message/MPIObjectBase.h" 27 28 #ifdef HAVE_LMY_ENGINE 29 //#include "Eigen/Dense" 30 #include "formic/utils/matrix.h" 31 #include "formic/utils/lmyengine/engine.h" 32 #endif 33 34 namespace qmcplusplus 35 { 36 class MCWalkerConfiguration; 37 class DescentEngine; 38 39 /** @ingroup QMCDrivers 40 * @brief Implements wave-function optimization 41 * 42 * Optimization by correlated sampling method with configurations 43 * generated from VMC. 44 */ 45 46 class QMCCostFunctionBase : public CostFunctionBase<QMCTraits::RealType>, public MPIObjectBase 47 { 48 public: 49 enum FieldIndex_OPT 50 { 51 LOGPSI_FIXED = 0, 52 LOGPSI_FREE = 1, 53 ENERGY_TOT = 2, 54 ENERGY_FIXED = 3, 55 ENERGY_NEW = 4, 56 REWEIGHT = 5 57 }; 58 enum SumIndex_OPT 59 { 60 SUM_E_BARE = 0, 61 SUM_ESQ_BARE, 62 SUM_ABSE_BARE, 63 SUM_E_WGT, 64 SUM_ESQ_WGT, 65 SUM_ABSE_WGT, 66 SUM_WGT, 67 SUM_WGTSQ, 68 SUM_INDEX_SIZE 69 }; 70 71 72 ///Constructor. 73 QMCCostFunctionBase(MCWalkerConfiguration& w, TrialWaveFunction& psi, QMCHamiltonian& h, Communicate* comm); 74 75 ///Destructor 76 virtual ~QMCCostFunctionBase(); 77 78 ///process xml node 79 bool put(xmlNodePtr cur); 80 void resetCostFunction(std::vector<xmlNodePtr>& cset); 81 ///Save opt parameters to HDF5 82 bool reportH5; 83 bool CI_Opt; 84 ///Path and name of the HDF5 prefix where CI coeffs are saved 85 std::string newh5; 86 ///assign optimization parameter i Params(int i)87 Return_t& Params(int i) { return OptVariables[i]; } 88 ///return optimization parameter i Params(int i)89 Return_t Params(int i) const { return OptVariables[i]; } getType(int i)90 int getType(int i) const { return OptVariables.getType(i); } 91 ///return the cost value for CGMinimization 92 Return_rt Cost(bool needGrad = true); 93 94 ///return the cost value for CGMinimization 95 Return_rt computedCost(); 96 void printEstimates(); 97 ///return the gradient of cost value for CGMinimization 98 virtual void GradCost(std::vector<Return_rt>& PGradient, 99 const std::vector<Return_rt>& PM, 100 Return_rt FiniteDiff = 0){}; 101 ///return the number of optimizable parameters getNumParams()102 inline int getNumParams() const { return OptVariables.size(); } 103 ///return the number of optimizable parameters getNumSamples()104 inline int getNumSamples() const { return NumSamples; } setNumSamples(int newNumSamples)105 inline void setNumSamples(int newNumSamples) { NumSamples = newNumSamples; } 106 ///reset the wavefunction 107 virtual void resetPsi(bool final_reset = false) = 0; 108 getParameterTypes(std::vector<int> & types)109 inline void getParameterTypes(std::vector<int>& types) { return OptVariablesForPsi.getParameterTypeList(types); } 110 111 ///dump the current parameters and other report 112 void Report(); 113 ///report parameters at the end 114 void reportParameters(); 115 116 ///report parameters in HDF5 at the end 117 void reportParametersH5(); 118 ///return the counter which keeps track of optimization steps getReportCounter()119 inline int getReportCounter() const { return ReportCounter; } 120 setWaveFunctionNode(xmlNodePtr cur)121 void setWaveFunctionNode(xmlNodePtr cur) { m_wfPtr = cur; } 122 123 void setTargetEnergy(Return_rt et); 124 setRootName(const std::string & aroot)125 void setRootName(const std::string& aroot) { RootName = aroot; } 126 setStream(std::ostream * os)127 void setStream(std::ostream* os) { msg_stream = os; } 128 129 void addCoefficients(xmlXPathContextPtr acontext, const char* cname); 130 131 void printCJParams(xmlNodePtr cur, std::string& rname); 132 133 void addCJParams(xmlXPathContextPtr acontext, const char* cname); 134 135 /** implement the virtual function 136 * @param x0 current parameters 137 * @param gr gradients or conjugate gradients 138 * @param dl return the displacelement to minimize the cost function 139 * @param val_proj projected cost 140 * 141 * If successful, any optimization object updates the parameters by x0 + dl*gr 142 * and proceeds with a new step. 143 */ 144 bool lineoptimization(const std::vector<Return_rt>& x0, 145 const std::vector<Return_rt>& gr, 146 Return_rt val0, 147 Return_rt& dl, 148 Return_rt& val_proj, 149 Return_rt& lambda_max); 150 151 virtual Return_rt fillOverlapHamiltonianMatrices(Matrix<Return_rt>& Left, Matrix<Return_rt>& Right) = 0; 152 153 #ifdef HAVE_LMY_ENGINE 154 Return_rt LMYEngineCost(const bool needDeriv, cqmc::engine::LMYEngine<Return_t>* EngineObj); 155 #endif 156 157 virtual void getConfigurations(const std::string& aroot) = 0; 158 159 virtual void checkConfigurations() = 0; 160 161 #ifdef HAVE_LMY_ENGINE 162 virtual void engine_checkConfigurations(cqmc::engine::LMYEngine<Return_t>* EngineObj, 163 DescentEngine& descentEngineObj, 164 const std::string& MinMethod) = 0; 165 166 #endif 167 168 void setRng(std::vector<RandomGenerator_t*>& r); 169 getneedGrads()170 inline bool getneedGrads() const { return needGrads; } 171 setneedGrads(bool tf)172 inline void setneedGrads(bool tf) { needGrads = tf; } setDMC()173 inline void setDMC() { vmc_or_dmc = 1.0; } 174 getParamName(int i)175 inline std::string getParamName(int i) const { return OptVariables.name(i); } 176 getOptVariables()177 inline const opt_variables_type& getOptVariables() const { return OptVariables; } 178 179 /// return variance after checkConfigurations getVariance()180 inline Return_rt getVariance() const 181 { 182 return SumValue[SUM_ESQ_WGT] / SumValue[SUM_WGT] - 183 (SumValue[SUM_E_WGT] / SumValue[SUM_WGT]) * (SumValue[SUM_E_WGT] / SumValue[SUM_WGT]); 184 } 185 186 protected: 187 ///walker ensemble 188 MCWalkerConfiguration& W; 189 190 ///trial function 191 TrialWaveFunction& Psi; 192 193 ///Hamiltonian 194 QMCHamiltonian& H; 195 196 ///if true, do not write the *.opt.#.xml 197 bool Write2OneXml; 198 ///if true, use analytic derivatives for the non-local potential component 199 bool useNLPPDeriv; 200 /** |E-E_T|^PowerE is used for the cost function 201 * 202 * default PowerE=1 203 */ 204 int PowerE; 205 ///number of times cost function evaluated 206 int NumCostCalls; 207 ///total number of samples to use in correlated sampling 208 int NumSamples; 209 ///total number of optimizable variables 210 int NumOptimizables; 211 ///counter for output 212 int ReportCounter; 213 ///weights for energy and variance in the cost function 214 Return_rt w_en, w_var, w_abs, w_w; 215 ///value of the cost function 216 Return_rt CostValue; 217 ///target energy 218 Return_rt Etarget; 219 ///real target energy with the Correlation Factor 220 Return_rt EtargetEff; 221 ///effective number of walkers 222 Return_rt NumWalkersEff; 223 ///fraction of the number of walkers below which the costfunction becomes invalid 224 Return_rt MinNumWalkers; 225 ///maximum weight beyond which the weight is set to 1 226 Return_rt MaxWeight; 227 ///current Average 228 Return_rt curAvg; 229 ///current Variance 230 Return_rt curVar; 231 ///current weighted average (correlated sampling) 232 Return_rt curAvg_w; 233 ///current weighted variance (correlated sampling) 234 Return_rt curVar_w; 235 ///current variance of SUM_ABSE_WGT/SUM_WGT 236 Return_rt curVar_abs; 237 238 Return_rt w_beta; 239 std::string GEVType; 240 Return_rt vmc_or_dmc; 241 bool needGrads; 242 ///whether we are targeting an excited state 243 std::string targetExcitedStr; 244 ///whether we are targeting an excited state 245 bool targetExcited; 246 ///the shift to use when targeting an excited state 247 double omega_shift; 248 249 ///list of optimizables 250 opt_variables_type OptVariables; 251 /** full list of optimizables 252 * 253 * The size of OptVariablesForPsi is equal to or larger than 254 * that of OptVariables due to the dependent variables. 255 * This is used for TrialWaveFunction::resetParameters and 256 * is normally the same as OptVariables. 257 */ 258 opt_variables_type OptVariablesForPsi; 259 // unchanged initial checked-in variables 260 opt_variables_type InitVariables; 261 /** index mapping for <equal> constraints 262 * 263 * - equalVarMap[i][0] : index in OptVariablesForPsi 264 * - equalVarMap[i][1] : index in OptVariables 265 */ 266 std::vector<TinyVector<int, 2>> equalVarMap; 267 /** index mapping for <negate> constraints 268 * 269 * - negateVarMap[i][0] : index in OptVariablesForPsi 270 * - negateVarMap[i][1] : index in OptVariables 271 */ 272 ///index mapping for <negative> constraints 273 std::vector<TinyVector<int, 2>> negateVarMap; 274 ///stream to which progress is sent 275 std::ostream* msg_stream; 276 ///xml node to be dumped 277 xmlNodePtr m_wfPtr; 278 ///document node to be dumped 279 xmlDocPtr m_doc_out; 280 ///parameters to be updated 281 std::map<std::string, xmlNodePtr> paramNodes; 282 ///coefficients to be updated 283 std::map<std::string, xmlNodePtr> coeffNodes; 284 ///attributes to be updated 285 std::map<std::string, std::pair<xmlNodePtr, std::string>> attribNodes; 286 ///string for the file root 287 std::string RootName; 288 ///Hamiltonians that depend on the optimization: KE 289 QMCHamiltonian H_KE; 290 291 ///Random number generators 292 std::vector<RandomGenerator_t*> RngSaved, MoverRng; 293 std::string includeNonlocalH; 294 295 296 /** Sum of energies and weights for averages 297 * 298 * SumValues[k] where k is one of SumIndex_opt 299 */ 300 std::vector<Return_rt> SumValue; 301 /** Saved properties of all the walkers 302 * 303 * Records(iw,field_id) returns the field_id value of the iw-th walker 304 * field_id is one of FieldIndex_opt 305 */ 306 Matrix<Return_rt> Records; 307 ///** Saved derivative properties and Hderivative properties of all the walkers 308 //*/ 309 //vector<std::vector<vector<Return_t> >* > DerivRecords; 310 //vector<std::vector<vector<Return_t> >* > HDerivRecords; 311 312 typedef ParticleSet::ParticleGradient_t ParticleGradient_t; 313 typedef ParticleSet::ParticleLaplacian_t ParticleLaplacian_t; 314 ///** Fixed Gradients , \f$\nabla\ln\Psi\f$, components */ 315 std::vector<ParticleGradient_t*> dLogPsi; 316 ///** Fixed Laplacian , \f$\nabla^2\ln\Psi\f$, components */ 317 std::vector<ParticleLaplacian_t*> d2LogPsi; 318 ///stream for debug 319 std::ostream* debug_stream; 320 321 bool checkParameters(); 322 void updateXmlNodes(); 323 324 325 virtual Return_rt correlatedSampling(bool needGrad = true) = 0; 326 327 #ifdef HAVE_LMY_ENGINE LMYEngineCost_detail(cqmc::engine::LMYEngine<Return_t> * EngineObj)328 virtual Return_rt LMYEngineCost_detail(cqmc::engine::LMYEngine<Return_t>* EngineObj) 329 { 330 APP_ABORT("NOT IMPLEMENTED"); 331 return 0; 332 } 333 #endif 334 }; 335 } // namespace qmcplusplus 336 #endif 337