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) 2019 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 // Luke Shulenburger, lshulen@sandia.gov, Sandia National Laboratories 13 // 14 // File created by: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign 15 ////////////////////////////////////////////////////////////////////////////////////// 16 17 18 /* 19 DO NOT MAKE PERMANENT EDITS IN THIS FILE 20 This file is generated from src/QMCWavefunctions/Jastrow/codegen/user_jastrow.py and UserFunctor.h.in 21 22 To make changes, edit UserFunctor.h.in and rerun user_jastrow.py. 23 */ 24 25 26 /** @file UserFunctor.h 27 * @brief User-defined functor 28 */ 29 #ifndef QMCPLUSPLUS_USERFUNCTOR_H 30 #define QMCPLUSPLUS_USERFUNCTOR_H 31 #include "Numerics/OptimizableFunctorBase.h" 32 #include "OhmmsData/AttributeSet.h" 33 #include <cmath> 34 // #include <vector> 35 #include "OhmmsPETE/TinyVector.h" 36 37 38 namespace qmcplusplus 39 { 40 /** Implements the function 41 * f = A*r/(B*r + 1) - A/B 42 43 * 44 */ 45 template<class T> 46 struct UserFunctor : public OptimizableFunctorBase 47 { 48 49 50 /// Is optimizable 51 bool Opt_A; 52 /// Value 53 real_type A; 54 /// Id in XML input 55 std::string ID_A; 56 57 58 /// Is optimizable 59 bool Opt_B; 60 /// Value 61 real_type B; 62 /// Id in XML input 63 std::string ID_B; 64 65 66 67 ///default constructor UserFunctorUserFunctor68 UserFunctor() { reset(); } 69 70 // void setCusp(real_type cusp) 71 setCuspUserFunctor72 void setCusp(real_type cusp) 73 { 74 A = cusp; 75 Opt_A = false; 76 reset(); 77 } 78 79 makeCloneUserFunctor80 OptimizableFunctorBase* makeClone() const { return new UserFunctor(*this); } 81 resetUserFunctor82 void reset() 83 { 84 } 85 86 87 88 // inline real_type evaluate(real_type r) const 89 evaluateUserFunctor90 inline real_type evaluate(real_type r) const { 91 return A*r/(B*r + 1) - A/B; 92 } 93 94 95 // const inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) const 96 evaluateUserFunctor97 inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2) const 98 { 99 dudr = -A*B*r/((B*r + 1)*(B*r + 1)) + A/(B*r + 1); 100 d2udr2 = 2*A*B*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1)); 101 return A*r/(B*r + 1) - A/B; 102 } 103 104 105 106 // inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type& d3udr3) const 107 evaluateUserFunctor108 inline real_type evaluate(real_type r, real_type& dudr, real_type& d2udr2, real_type& d3udr3) const 109 { 110 dudr = -A*B*r/((B*r + 1)*(B*r + 1)) + A/(B*r + 1); 111 d2udr2 = 2*A*B*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1)); 112 d3udr3 = 6*A*((B)*(B))*(-B*r/(B*r + 1) + 1)/std::pow(B*r + 1, 3); 113 return A*r/(B*r + 1) - A/B; 114 } 115 116 117 118 evaluateVUserFunctor119 inline real_type evaluateV(const int iat, 120 const int iStart, 121 const int iEnd, 122 const T* restrict _distArray, 123 T* restrict distArrayCompressed) const 124 { 125 // specialized evaluation loop? 126 real_type sum(0); 127 for (int idx = iStart; idx < iEnd; idx++) 128 if (idx != iat) 129 sum += evaluate(_distArray[idx]); 130 return sum; 131 } 132 evaluateVGLUserFunctor133 inline void evaluateVGL(const int iat, 134 const int iStart, 135 const int iEnd, 136 const T* distArray, 137 T* restrict valArray, 138 T* restrict gradArray, 139 T* restrict laplArray, 140 T* restrict distArrayCompressed, 141 int* restrict distIndices) const 142 { 143 // specialized evaluation loop? 144 for (int idx = iStart; idx < iEnd; idx++) 145 { 146 valArray[idx] = evaluate(distArray[idx], gradArray[idx], laplArray[idx]); 147 gradArray[idx] /= distArray[idx]; 148 } 149 if (iat >= iStart && iat < iEnd) 150 valArray[iat] = gradArray[iat] = laplArray[iat] = T(0); 151 } 152 fUserFunctor153 inline real_type f(real_type r) { return evaluate(r); } 154 dfUserFunctor155 inline real_type df(real_type r) 156 { 157 real_type dudr, d2udr2; 158 real_type res = evaluate(r, dudr, d2udr2); 159 return dudr; 160 } 161 162 // inline bool evaluateDerivatives(real_type r, std::vector<TinyVector<real_type, 3>>& derivs) 163 evaluateDerivativesUserFunctor164 inline bool evaluateDerivatives(real_type r, std::vector<TinyVector<real_type, 3>>& derivs) 165 { 166 int i = 0; 167 168 if (Opt_A) 169 { 170 derivs[i][0] = r/(B*r + 1) - 1/B; 171 derivs[i][1] = -B*r/((B*r + 1)*(B*r + 1)) + 1.0/(B*r + 1); 172 derivs[i][2] = 2*B*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1)); 173 ++i; 174 } 175 176 if (Opt_B) 177 { 178 derivs[i][0] = -A*((r)*(r))/((B*r + 1)*(B*r + 1)) + A/((B)*(B)); 179 derivs[i][1] = 2*A*B*((r)*(r))/std::pow(B*r + 1, 3) - 2*A*r/((B*r + 1)*(B*r + 1)); 180 derivs[i][2] = -4*A*B*r*(B*r/(B*r + 1) - 1)/std::pow(B*r + 1, 3) + 2*A*B*(-B*((r)*(r))/((B*r + 1)*(B*r + 1)) + r/(B*r + 1))/((B*r + 1)*(B*r + 1)) + 2*A*(B*r/(B*r + 1) - 1)/((B*r + 1)*(B*r + 1)); 181 ++i; 182 } 183 184 return true; 185 } 186 187 188 // inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs) 189 190 /// compute derivatives with respect to variational parameters evaluateDerivativesUserFunctor191 inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs) 192 { 193 int i = 0; 194 195 if (Opt_A) 196 { 197 derivs[i] = r/(B*r + 1) - 1/B; 198 ++i; 199 } 200 201 if (Opt_B) 202 { 203 derivs[i] = -A*((r)*(r))/((B*r + 1)*(B*r + 1)) + A/((B)*(B)); 204 ++i; 205 } 206 207 208 return true; 209 } 210 211 212 // bool put(xmlNodePtr cur) 213 putUserFunctor214 bool put(xmlNodePtr cur) 215 { 216 cur = cur->xmlChildrenNode; 217 while (cur != NULL) 218 { 219 std::string cname((const char*)(cur->name)); 220 if (cname == "var") //only accept var 221 { 222 std::string id_in; 223 std::string p_name; 224 OhmmsAttributeSet rAttrib; 225 rAttrib.add(id_in, "id"); 226 rAttrib.add(p_name, "name"); 227 rAttrib.put(cur); 228 229 if (p_name == "A") 230 { 231 ID_A = id_in; 232 putContent(A, cur); 233 Opt_A = true; 234 } 235 236 if (p_name == "B") 237 { 238 ID_B = id_in; 239 putContent(B, cur); 240 Opt_B = true; 241 } 242 243 } 244 cur = cur->next; 245 } 246 reset(); 247 myVars.clear(); 248 249 if (Opt_A) 250 myVars.insert(ID_A, A, Opt_A, optimize::OTHER_P); 251 252 if (Opt_B) 253 myVars.insert(ID_B, B, Opt_B, optimize::OTHER_P); 254 255 return true; 256 } 257 258 checkInVariablesUserFunctor259 void checkInVariables(opt_variables_type& active) 260 { 261 active.insertFrom(myVars); 262 //myVars.print(std::cout); 263 } 264 checkOutVariablesUserFunctor265 void checkOutVariables(const opt_variables_type& active) 266 { 267 myVars.getIndex(active); 268 //myVars.print(std::cout); 269 } 270 271 //void resetParameters(const opt_variables_type& active) 272 resetParametersUserFunctor273 void resetParameters(const opt_variables_type& active) 274 { 275 if (myVars.size()) 276 { 277 int ia = myVars.where(0); 278 if (ia > -1) 279 { 280 int i = 0; 281 282 if (Opt_A) 283 A = std::real(myVars[i++] = active[ia++]); 284 285 if (Opt_B) 286 B = std::real(myVars[i++] = active[ia++]); 287 288 } 289 reset(); 290 } 291 } 292 293 }; 294 295 296 297 298 } // namespace qmcplusplus 299 #endif 300