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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign 8 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign 9 // Ye Luo, yeluo@anl.gov, Argonne National Laboratory 10 // Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory 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_COULOMBPOTENTIAL_H 18 #define QMCPLUSPLUS_COULOMBPOTENTIAL_H 19 #include "Particle/ParticleSet.h" 20 #include "Particle/WalkerSetRef.h" 21 #include "Particle/DistanceTableData.h" 22 #include "QMCHamiltonians/ForceBase.h" 23 #include "QMCHamiltonians/OperatorBase.h" 24 #include <numeric> 25 26 namespace qmcplusplus 27 { 28 /** CoulombPotential 29 * @tparam T type of the elementary data 30 * 31 * Hamiltonian operator for the Coulomb interaction for both AA and AB type for open systems. 32 */ 33 template<typename T> 34 struct CoulombPotential : public OperatorBase, public ForceBase 35 { 36 ///source particle set 37 ParticleSet& Pa; 38 ///target particle set 39 ParticleSet& Pb; 40 ///distance table index 41 const int myTableIndex; 42 ///true if the table is AA 43 const bool is_AA; 44 ///true, if CoulombAA for quantum particleset 45 bool is_active; 46 ///number of centers 47 int nCenters; 48 #if !defined(REMOVE_TRACEMANAGER) 49 ///single particle trace samples 50 Array<TraceReal, 1>* Va_sample; 51 Array<TraceReal, 1>* Vb_sample; 52 #endif 53 54 /// Flag for whether to compute forces or not 55 bool ComputeForces; 56 57 /** constructor for AA 58 * @param s source particleset 59 * @param active if true, new Value is computed whenver evaluate is used. 60 * @param computeForces if true, computes forces between inactive species 61 */ 62 inline CoulombPotential(ParticleSet& s, bool active, bool computeForces, bool copy = false) ForceBaseCoulombPotential63 : ForceBase(s, s), 64 Pa(s), 65 Pb(s), 66 myTableIndex(s.addTable(s)), 67 is_AA(true), 68 is_active(active), 69 ComputeForces(computeForces) 70 { 71 set_energy_domain(potential); 72 two_body_quantum_domain(s, s); 73 nCenters = s.getTotalNum(); 74 prefix = "F_AA"; 75 76 if (!is_active) //precompute the value 77 { 78 if (!copy) 79 s.update(); 80 Value = evaluateAA(s.getDistTable(myTableIndex), s.Z.first_address()); 81 if (ComputeForces) 82 evaluateAAForces(s.getDistTable(myTableIndex), s.Z.first_address()); 83 } 84 } 85 86 /** constructor for AB 87 * @param s source particleset 88 * @param t target particleset 89 * @param active if true, new Value is computed whenver evaluate is used. 90 * @param ComputeForces is not implemented for AB 91 */ 92 inline CoulombPotential(ParticleSet& s, ParticleSet& t, bool active, bool copy = false) ForceBaseCoulombPotential93 : ForceBase(s, t), 94 Pa(s), 95 Pb(t), 96 myTableIndex(t.addTable(s)), 97 is_AA(false), 98 is_active(active), 99 ComputeForces(false) 100 { 101 set_energy_domain(potential); 102 two_body_quantum_domain(s, t); 103 nCenters = s.getTotalNum(); 104 } 105 106 #if !defined(REMOVE_TRACEMANAGER) contribute_particle_quantitiesCoulombPotential107 virtual void contribute_particle_quantities() { request.contribute_array(myName); } 108 checkout_particle_quantitiesCoulombPotential109 virtual void checkout_particle_quantities(TraceManager& tm) 110 { 111 streaming_particles = request.streaming_array(myName); 112 if (streaming_particles) 113 { 114 Va_sample = tm.checkout_real<1>(myName, Pa); 115 if (!is_AA) 116 { 117 Vb_sample = tm.checkout_real<1>(myName, Pb); 118 } 119 else if (!is_active) 120 evaluate_spAA(Pa.getDistTable(myTableIndex), Pa.Z.first_address()); 121 } 122 } 123 delete_particle_quantitiesCoulombPotential124 virtual void delete_particle_quantities() 125 { 126 if (streaming_particles) 127 { 128 delete Va_sample; 129 if (!is_AA) 130 delete Vb_sample; 131 } 132 } 133 #endif 134 addObservablesCoulombPotential135 inline void addObservables(PropertySetType& plist, BufferType& collectables) 136 { 137 addValue(plist); 138 if (ComputeForces) 139 addObservablesF(plist); 140 } 141 142 /** evaluate AA-type interactions */ evaluateAACoulombPotential143 inline T evaluateAA(const DistanceTableData& d, const ParticleScalar_t* restrict Z) 144 { 145 T res = 0.0; 146 #if !defined(REMOVE_TRACEMANAGER) 147 if (streaming_particles) 148 res = evaluate_spAA(d, Z); 149 else 150 #endif 151 for (size_t iat = 1; iat < nCenters; ++iat) 152 { 153 const auto& dist = d.getDistRow(iat); 154 T q = Z[iat]; 155 for (size_t j = 0; j < iat; ++j) 156 res += q * Z[j] / dist[j]; 157 } 158 return res; 159 } 160 161 162 /** evaluate AA-type forces */ evaluateAAForcesCoulombPotential163 inline void evaluateAAForces(const DistanceTableData& d, const ParticleScalar_t* restrict Z) 164 { 165 forces = 0.0; 166 for (size_t iat = 1; iat < nCenters; ++iat) 167 { 168 const auto& dist = d.getDistRow(iat); 169 const auto& displ = d.getDisplRow(iat); 170 T q = Z[iat]; 171 for (size_t j = 0; j < iat; ++j) 172 { 173 forces[iat] += -q * Z[j] * displ[j] / (dist[j] * dist[j] * dist[j]); 174 forces[j] -= -q * Z[j] * displ[j] / (dist[j] * dist[j] * dist[j]); 175 } 176 } 177 } 178 179 180 /** JNKIM: Need to check the precision */ evaluateABCoulombPotential181 inline T evaluateAB(const DistanceTableData& d, 182 const ParticleScalar_t* restrict Za, 183 const ParticleScalar_t* restrict Zb) 184 { 185 constexpr T czero(0); 186 T res = czero; 187 #if !defined(REMOVE_TRACEMANAGER) 188 if (streaming_particles) 189 res = evaluate_spAB(d, Za, Zb); 190 else 191 #endif 192 { 193 const size_t nTargets = d.targets(); 194 for (size_t b = 0; b < nTargets; ++b) 195 { 196 const auto& dist = d.getDistRow(b); 197 T e = czero; 198 for (size_t a = 0; a < nCenters; ++a) 199 e += Za[a] / dist[a]; 200 res += e * Zb[b]; 201 } 202 } 203 return res; 204 } 205 206 207 #if !defined(REMOVE_TRACEMANAGER) 208 /** evaluate AA-type interactions */ evaluate_spAACoulombPotential209 inline T evaluate_spAA(const DistanceTableData& d, const ParticleScalar_t* restrict Z) 210 { 211 T res = 0.0; 212 T pairpot; 213 Array<RealType, 1>& Va_samp = *Va_sample; 214 Va_samp = 0.0; 215 for (size_t iat = 1; iat < nCenters; ++iat) 216 { 217 const auto& dist = d.getDistRow(iat); 218 T q = Z[iat]; 219 for (size_t j = 0; j < iat; ++j) 220 { 221 pairpot = 0.5 * q * Z[j] / dist[j]; 222 Va_samp(iat) += pairpot; 223 Va_samp(j) += pairpot; 224 res += pairpot; 225 } 226 } 227 res *= 2.0; 228 #if defined(TRACE_CHECK) 229 auto sptmp = streaming_particles; 230 streaming_particles = false; 231 T Vnow = res; 232 T Vsum = Va_samp.sum(); 233 T Vorig = evaluateAA(d, Z); 234 if (std::abs(Vorig - Vnow) > TraceManager::trace_tol) 235 { 236 app_log() << "versiontest: CoulombPotential::evaluateAA()" << std::endl; 237 app_log() << "versiontest: orig:" << Vorig << std::endl; 238 app_log() << "versiontest: mod:" << Vnow << std::endl; 239 APP_ABORT("Trace check failed"); 240 } 241 if (std::abs(Vsum - Vnow) > TraceManager::trace_tol) 242 { 243 app_log() << "accumtest: CoulombPotential::evaluateAA()" << std::endl; 244 app_log() << "accumtest: tot:" << Vnow << std::endl; 245 app_log() << "accumtest: sum:" << Vsum << std::endl; 246 APP_ABORT("Trace check failed"); 247 } 248 streaming_particles = sptmp; 249 #endif 250 return res; 251 } 252 253 evaluate_spABCoulombPotential254 inline T evaluate_spAB(const DistanceTableData& d, 255 const ParticleScalar_t* restrict Za, 256 const ParticleScalar_t* restrict Zb) 257 { 258 T res = 0.0; 259 T pairpot; 260 Array<RealType, 1>& Va_samp = *Va_sample; 261 Array<RealType, 1>& Vb_samp = *Vb_sample; 262 Va_samp = 0.0; 263 Vb_samp = 0.0; 264 const size_t nTargets = d.targets(); 265 for (size_t b = 0; b < nTargets; ++b) 266 { 267 const auto& dist = d.getDistRow(b); 268 T z = 0.5 * Zb[b]; 269 for (size_t a = 0; a < nCenters; ++a) 270 { 271 pairpot = z * Za[a] / dist[a]; 272 Va_samp(a) += pairpot; 273 Vb_samp(b) += pairpot; 274 res += pairpot; 275 } 276 } 277 res *= 2.0; 278 279 #if defined(TRACE_CHECK) 280 auto sptmp = streaming_particles; 281 streaming_particles = false; 282 T Vnow = res; 283 T Vasum = Va_samp.sum(); 284 T Vbsum = Vb_samp.sum(); 285 T Vsum = Vasum + Vbsum; 286 T Vorig = evaluateAB(d, Za, Zb); 287 if (std::abs(Vorig - Vnow) > TraceManager::trace_tol) 288 { 289 app_log() << "versiontest: CoulombPotential::evaluateAB()" << std::endl; 290 app_log() << "versiontest: orig:" << Vorig << std::endl; 291 app_log() << "versiontest: mod:" << Vnow << std::endl; 292 APP_ABORT("Trace check failed"); 293 } 294 if (std::abs(Vsum - Vnow) > TraceManager::trace_tol) 295 { 296 app_log() << "accumtest: CoulombPotential::evaluateAB()" << std::endl; 297 app_log() << "accumtest: tot:" << Vnow << std::endl; 298 app_log() << "accumtest: sum:" << Vsum << std::endl; 299 APP_ABORT("Trace check failed"); 300 } 301 if (std::abs(Vasum - Vbsum) > TraceManager::trace_tol) 302 { 303 app_log() << "sharetest: CoulombPotential::evaluateAB()" << std::endl; 304 app_log() << "sharetest: a share:" << Vasum << std::endl; 305 app_log() << "sharetest: b share:" << Vbsum << std::endl; 306 APP_ABORT("Trace check failed"); 307 } 308 streaming_particles = sptmp; 309 #endif 310 return res; 311 } 312 #endif 313 314 resetTargetParticleSetCoulombPotential315 void resetTargetParticleSet(ParticleSet& P) 316 { 317 //myTableIndex is the same 318 } 319 ~CoulombPotentialCoulombPotential320 ~CoulombPotential() {} 321 update_sourceCoulombPotential322 void update_source(ParticleSet& s) 323 { 324 if (is_AA) 325 { 326 Value = evaluateAA(s.getDistTable(myTableIndex), s.Z.first_address()); 327 } 328 } 329 evaluateCoulombPotential330 inline Return_t evaluate(ParticleSet& P) 331 { 332 if (is_active) 333 { 334 if (is_AA) 335 Value = evaluateAA(P.getDistTable(myTableIndex), P.Z.first_address()); 336 else 337 Value = evaluateAB(P.getDistTable(myTableIndex), Pa.Z.first_address(), P.Z.first_address()); 338 } 339 return Value; 340 } 341 evaluateWithIonDerivsCoulombPotential342 inline Return_t evaluateWithIonDerivs(ParticleSet& P, 343 ParticleSet& ions, 344 TrialWaveFunction& psi, 345 ParticleSet::ParticlePos_t& hf_terms, 346 ParticleSet::ParticlePos_t& pulay_terms) 347 { 348 if (is_active) 349 Value = evaluate(P); // No forces for the active 350 else 351 hf_terms -= forces; // No Pulay here 352 return Value; 353 } 354 putCoulombPotential355 bool put(xmlNodePtr cur) { return true; } 356 getCoulombPotential357 bool get(std::ostream& os) const 358 { 359 if (myTableIndex) 360 os << "CoulombAB source=" << Pa.getName() << std::endl; 361 else 362 os << "CoulombAA source/target " << Pa.getName() << std::endl; 363 return true; 364 } 365 setObservablesCoulombPotential366 void setObservables(PropertySetType& plist) 367 { 368 OperatorBase::setObservables(plist); 369 if (ComputeForces) 370 setObservablesF(plist); 371 } 372 setParticlePropertyListCoulombPotential373 void setParticlePropertyList(PropertySetType& plist, int offset) 374 { 375 OperatorBase::setParticlePropertyList(plist, offset); 376 if (ComputeForces) 377 setParticleSetF(plist, offset); 378 } 379 makeCloneCoulombPotential380 OperatorBase* makeClone(ParticleSet& qp, TrialWaveFunction& psi) 381 { 382 if (is_AA) 383 { 384 if (is_active) 385 return new CoulombPotential(qp, true, ComputeForces); 386 else 387 // Ye Luo April 16th, 2015 388 // avoid recomputing ion-ion DistanceTable when reusing ParticleSet 389 return new CoulombPotential(Pa, false, ComputeForces, true); 390 } 391 else 392 return new CoulombPotential(Pa, qp, true); 393 } 394 }; 395 396 } // namespace qmcplusplus 397 #endif 398