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) 2020 QMCPACK developers.
6 //
7 // File developed by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories
8 //
9 // File created by: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories
10 //////////////////////////////////////////////////////////////////////////////////////
11
12 #include "Particle/DistanceTableData.h"
13 #include "SOECPotential.h"
14 #include "Utilities/IteratorUtility.h"
15
16 namespace qmcplusplus
17 {
18 /** constructor
19 *\param ionic positions
20 *\param els electronic poitions
21 *\param psi Trial wave function
22 */
SOECPotential(ParticleSet & ions,ParticleSet & els,TrialWaveFunction & psi)23 SOECPotential::SOECPotential(ParticleSet& ions, ParticleSet& els, TrialWaveFunction& psi)
24 : myRNG(nullptr), IonConfig(ions), Psi(psi), Peln(els), ElecNeighborIons(els), IonNeighborElecs(ions)
25 {
26 set_energy_domain(potential);
27 two_body_quantum_domain(ions, els);
28 myTableIndex = els.addTable(ions);
29 NumIons = ions.getTotalNum();
30 PP.resize(NumIons, nullptr);
31 PPset.resize(IonConfig.getSpeciesSet().getTotalNum());
32 }
33
resetTargetParticleSet(ParticleSet & P)34 void SOECPotential::resetTargetParticleSet(ParticleSet& P) {}
35
evaluate(ParticleSet & P)36 SOECPotential::Return_t SOECPotential::evaluate(ParticleSet& P)
37 {
38 Value = 0.0;
39 for (int ipp = 0; ipp < PPset.size(); ipp++)
40 if (PPset[ipp])
41 PPset[ipp]->randomize_grid(*myRNG);
42 const auto& myTable = P.getDistTable(myTableIndex);
43 for (int iat = 0; iat < NumIons; iat++)
44 IonNeighborElecs.getNeighborList(iat).clear();
45 for (int jel = 0; jel < P.getTotalNum(); jel++)
46 ElecNeighborIons.getNeighborList(jel).clear();
47
48 for (int jel = 0; jel < P.getTotalNum(); jel++)
49 {
50 const auto& dist = myTable.getDistRow(jel);
51 const auto& displ = myTable.getDisplRow(jel);
52 std::vector<int>& NeighborIons = ElecNeighborIons.getNeighborList(jel);
53 for (int iat = 0; iat < NumIons; iat++)
54 if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax())
55 {
56 RealType pairpot = PP[iat]->evaluateOne(P, iat, Psi, jel, dist[iat], -displ[iat]);
57 Value += pairpot;
58 NeighborIons.push_back(iat);
59 IonNeighborElecs.getNeighborList(iat).push_back(jel);
60 }
61 }
62 return Value;
63 }
64
makeClone(ParticleSet & qp,TrialWaveFunction & psi)65 OperatorBase* SOECPotential::makeClone(ParticleSet& qp, TrialWaveFunction& psi)
66 {
67 SOECPotential* myclone = new SOECPotential(IonConfig, qp, psi);
68 for (int ig = 0; ig < PPset.size(); ++ig)
69 if (PPset[ig])
70 myclone->addComponent(ig, std::unique_ptr<SOECPComponent>(PPset[ig]->makeClone(qp)));
71 return myclone;
72 }
73
addComponent(int groupID,std::unique_ptr<SOECPComponent> && ppot)74 void SOECPotential::addComponent(int groupID, std::unique_ptr<SOECPComponent>&& ppot)
75 {
76 for (int iat = 0; iat < PP.size(); iat++)
77 if (IonConfig.GroupID[iat] == groupID)
78 PP[iat] = ppot.get();
79 PPset[groupID] = std::move(ppot);
80 }
81
82 } // namespace qmcplusplus
83