Merge pull request #2593 from ye-luo/LRCoulombSingleton-unique_ptr

Define LRCoulombSingleton handler ownership.

src/Particle/LongRange/LRCoulombSingleton.cpp

@@ -28,8 +28,8 @@
 namespace qmcplusplus
 {
 //initialization of the static data
-LRCoulombSingleton::LRHandlerType* LRCoulombSingleton::CoulombHandler      = 0;
-LRCoulombSingleton::LRHandlerType* LRCoulombSingleton::CoulombDerivHandler = 0;
+std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::CoulombHandler;
+std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::CoulombDerivHandler;
 LRCoulombSingleton::lr_type LRCoulombSingleton::this_lr_type = ESLER;
 /** CoulombFunctor
  *
@@ -108,37 +108,37 @@ struct PseudoCoulombFunctor
 };
 
 
-LRCoulombSingleton::LRHandlerType* LRCoulombSingleton::getHandler(ParticleSet& ref)
+std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::getHandler(ParticleSet& ref)
 {
   if (CoulombHandler == 0)
   {
 #if OHMMS_DIM == 3
     if (ref.SK->SuperCellEnum == SUPERCELL_SLAB)
     {
       app_log() << "\n   Creating CoulombHandler using quasi-2D Ewald method for the slab. " << std::endl;
-      CoulombHandler = new EwaldHandler(ref);
+      CoulombHandler = std::make_unique<EwaldHandler>(ref);
     }
     else //if(ref.LRBox.SuperCellEnum == SUPERCELL_BULK)
     {
-      if(this_lr_type == ESLER)
-      {  
+      if (this_lr_type == ESLER)
+      {
         app_log() << "\n  Creating CoulombHandler with the Esler Optimized Breakup. " << std::endl;
-        CoulombHandler = new LRHandlerTemp<CoulombFunctor<mRealType>, LPQHIBasis>(ref);
+        CoulombHandler = std::make_unique<LRHandlerTemp<CoulombFunctor<mRealType>, LPQHIBasis>>(ref);
       }
       else if (this_lr_type == EWALD)
       {
         app_log() << "\n  Creating CoulombHandler with the 3D Ewald Breakup. " << std::endl;
-        CoulombHandler= new EwaldHandler3D(ref); 
+        CoulombHandler = std::make_unique<EwaldHandler3D>(ref);
       }
       else if (this_lr_type == NATOLI)
       {
         app_log() << "\n  Creating CoulombHandler with the Natoli Optimized Breakup. " << std::endl;
-        CoulombHandler = new LRHandlerSRCoulomb<CoulombFunctor<mRealType>, LPQHISRCoulombBasis>(ref);
+        CoulombHandler = std::make_unique<LRHandlerSRCoulomb<CoulombFunctor<mRealType>, LPQHISRCoulombBasis>>(ref);
       }
       else
       {
         APP_ABORT("\n  Long range breakup method not recognized.\n");
-      } 
+      }
     }
 //        else if(ref.LRBox.SuperCellEnum == SUPERCELL_SLAB)
 //        {
@@ -147,18 +147,19 @@ LRCoulombSingleton::LRHandlerType* LRCoulombSingleton::getHandler(ParticleSet& r
 //        }
 #elif OHMMS_DIM == 2
     app_log() << "\n   Creating CoulombHandler using 2D Ewald method. " << std::endl;
-    CoulombHandler = new TwoDEwaldHandler(ref);
+    CoulombHandler = std::make_unique<TwoDEwaldHandler>(ref);
 #endif
     CoulombHandler->initBreakup(ref);
-    return CoulombHandler;
+    return std::unique_ptr<LRHandlerType>(CoulombHandler->makeClone(ref));
   }
   else
   {
     app_log() << "  Clone CoulombHandler. " << std::endl;
-    return CoulombHandler->makeClone(ref);
+    return std::unique_ptr<LRHandlerType>(CoulombHandler->makeClone(ref));
   }
 }
-LRCoulombSingleton::LRHandlerType* LRCoulombSingleton::getDerivHandler(ParticleSet& ref)
+
+std::unique_ptr<LRCoulombSingleton::LRHandlerType> LRCoulombSingleton::getDerivHandler(ParticleSet& ref)
 {
 #if OHMMS_DIM != 3
   APP_ABORT("energy derivative implemented for 3D only");
@@ -169,33 +170,28 @@ LRCoulombSingleton::LRHandlerType* LRCoulombSingleton::getDerivHandler(ParticleS
     if (this_lr_type == EWALD)
     {
       app_log() << "\n  Creating CoulombDerivHandler with the 3D Ewald Breakup. " << std::endl;
-      CoulombDerivHandler= new EwaldHandler3D(ref);
-    } 
+      CoulombDerivHandler = std::make_unique<EwaldHandler3D>(ref);
+    }
     else if (this_lr_type == NATOLI)
     {
       app_log() << "\n  Creating CoulombDerivHandler with the Natoli Optimized Breakup. " << std::endl;
-      CoulombDerivHandler = new LRHandlerSRCoulomb<CoulombFunctor<mRealType>, LPQHISRCoulombBasis>(ref);
-    } 
+      CoulombDerivHandler = std::make_unique<LRHandlerSRCoulomb<CoulombFunctor<mRealType>, LPQHISRCoulombBasis>>(ref);
+    }
     else if (this_lr_type == ESLER)
     {
-      APP_ABORT("\n  Derivatives are not supported with Esler Optimized Breakup.\n"); 
+      APP_ABORT("\n  Derivatives are not supported with Esler Optimized Breakup.\n");
     }
     else
     {
       APP_ABORT("\n  Long range breakup method for derivatives not recognized.\n");
     }
-    // CoulombDerivHandler = new LRDerivHandler<CoulombFunctor<mRealType>, LPQHIBasis> (ref);
-    //CoulombDerivHandler= new EwaldHandler(ref);
     CoulombDerivHandler->initBreakup(ref);
-    //return CoulombDerivHandler;
-
-    return CoulombDerivHandler;
+    return std::unique_ptr<LRHandlerType>(CoulombDerivHandler->makeClone(ref));
   }
   else
   {
     app_log() << "  Clone CoulombDerivHandler. " << std::endl;
-    return CoulombDerivHandler->makeClone(ref);
-    // return CoulombDerivHandler;
+    return std::unique_ptr<LRHandlerType>(CoulombDerivHandler->makeClone(ref));
   }
 }
 

src/Particle/LongRange/LRCoulombSingleton.h

@@ -19,6 +19,7 @@
 #ifndef QMCPLUSPLUS_LRCOULOMBSINGLETON_H
 #define QMCPLUSPLUS_LRCOULOMBSINGLETON_H
 
+#include <memory>
 #include <config.h>
 #include "LongRange/LRHandlerTemp.h"
 #include "LongRange/LRHandlerSRCoulomb.h"
@@ -37,16 +38,21 @@ struct LRCoulombSingleton
   typedef LinearGrid<pRealType> GridType;
   typedef OneDimCubicSpline<pRealType> RadFunctorType;
 
-  enum lr_type {ESLER=0, EWALD, NATOLI};
+  enum lr_type
+  {
+    ESLER = 0,
+    EWALD,
+    NATOLI
+  };
   static lr_type this_lr_type;
   ///Stores the energ optimized LR handler.
-  static LRHandlerType* CoulombHandler;
+  static std::unique_ptr<LRHandlerType> CoulombHandler;
   ///Stores the force/stress optimized LR handler.
-  static LRHandlerType* CoulombDerivHandler;
+  static std::unique_ptr<LRHandlerType> CoulombDerivHandler;
   ///This returns an energy optimized LR handler.  If non existent, it creates one.
-  static LRHandlerType* getHandler(ParticleSet& ref);
+  static std::unique_ptr<LRHandlerType> getHandler(ParticleSet& ref);
   ///This returns a force/stress optimized LR handler.  If non existent, it creates one.
-  static LRHandlerType* getDerivHandler(ParticleSet& ref);
+  static std::unique_ptr<LRHandlerType> getDerivHandler(ParticleSet& ref);
 
   //The following two helper functions are provided to spline the short-range component
   //of the coulomb potential and its derivative.  This is much faster than evaluating

src/QMCHamiltonians/CoulombPBCAA.cpp

@@ -24,10 +24,8 @@ namespace qmcplusplus
 {
 CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
     : ForceBase(ref, ref),
-      AA(0),
       myGrid(0),
       rVs(0),
-      dAA(0),
       myGridforce(0),
       rVsforce(0),
       is_active(active),
@@ -308,14 +306,14 @@ void CoulombPBCAA::initBreakup(ParticleSet& P)
   myRcut  = AA->get_rc(); //Basis.get_rc();
   if (rVs == 0)
   {
-    rVs = LRCoulombSingleton::createSpline4RbyVs(AA, myRcut, myGrid);
+    rVs = LRCoulombSingleton::createSpline4RbyVs(AA.get(), myRcut, myGrid);
   }
   if (ComputeForces)
   {
     dAA = LRCoulombSingleton::getDerivHandler(P);
     if (rVsforce == 0)
     {
-      rVsforce = LRCoulombSingleton::createSpline4RbyVs(dAA, myRcut, myGridforce);
+      rVsforce = LRCoulombSingleton::createSpline4RbyVs(dAA.get(), myRcut, myGridforce);
     }
   }
 

src/QMCHamiltonians/CoulombPBCAA.h

@@ -36,12 +36,18 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
   typedef LRCoulombSingleton::RadFunctorType RadFunctorType;
   typedef LRHandlerType::mRealType mRealType;
 
+  // using shared_ptr on AA, dAA is a compromise
+  // When ion-ion is_active = false, makeClone calls the copy constructor.
+  // AA, dAA are shared between clones.
+  // When elec-elec is_active = true, makeClone calls the constructor
+  // AA, dAA are not shared between clones and behave more like unique_ptr
+
   // energy-optimized
-  LRHandlerType* AA;
+  std::shared_ptr<LRHandlerType> AA;
   GridType* myGrid;
   RadFunctorType* rVs;
   // force-optimized
-  LRHandlerType* dAA;
+  std::shared_ptr<LRHandlerType> dAA;
   GridType* myGridforce;
   RadFunctorType* rVsforce;
 

src/QMCHamiltonians/CoulombPBCAB.cpp

@@ -470,7 +470,7 @@ void CoulombPBCAB::initBreakup(ParticleSet& P)
   // create the spline function for the short-range part assuming pure potential
   if (V0 == nullptr)
   {
-    V0 = LRCoulombSingleton::createSpline4RbyVs(AB, myRcut, myGrid);
+    V0 = LRCoulombSingleton::createSpline4RbyVs(AB.get(), myRcut, myGrid);
     if (Vat.size())
     {
       APP_ABORT("CoulombPBCAB::initBreakup.  Vat is not empty\n");
@@ -484,9 +484,9 @@ void CoulombPBCAB::initBreakup(ParticleSet& P)
   {
     dAB = LRCoulombSingleton::getDerivHandler(P);
     if (fV0 == nullptr)
-      fV0 = LRCoulombSingleton::createSpline4RbyVs(dAB, myRcut, myGrid);
+      fV0 = LRCoulombSingleton::createSpline4RbyVs(dAB.get(), myRcut, myGrid);
     if (dfV0 == nullptr)
-      dfV0 = LRCoulombSingleton::createSpline4RbyVsDeriv(dAB, myRcut, myGrid);
+      dfV0 = LRCoulombSingleton::createSpline4RbyVsDeriv(dAB.get(), myRcut, myGrid);
     if (fVat.size())
     {
       APP_ABORT("CoulombPBCAB::initBreakup.  Vat is not empty\n");

src/QMCHamiltonians/CoulombPBCAB.h

@@ -42,9 +42,9 @@ struct CoulombPBCAB : public OperatorBase, public ForceBase
   ///source particle set
   ParticleSet& PtclA;
   ///long-range Handler
-  LRHandlerType* AB;
+  std::unique_ptr<LRHandlerType> AB;
   ///long-range derivative handler
-  LRHandlerType* dAB;
+  std::unique_ptr<LRHandlerType> dAB;
   ///locator of the distance table
   const int myTableIndex;
   ///number of species of A particle set

src/QMCHamiltonians/ForceChiesaPBCAA.h

@@ -38,7 +38,7 @@ struct ForceChiesaPBCAA : public OperatorBase, public ForceBase
   ///source particle set
   ParticleSet& PtclA;
   ///long-range Handler
-  LRHandlerType* dAB;
+  std::unique_ptr<LRHandlerType> dAB;
   ///number of species of A particle set
   int NumSpeciesA;
   ///number of species of B particle set

src/QMCHamiltonians/tests/test_coulomb_CUDA.cpp

@@ -82,7 +82,7 @@ TEST_CASE("Coulomb PBC A-B CUDA", "[hamiltonian][CUDA]")
   ParticleSetPool ptcl = ParticleSetPool(c);
 
 
-  CoulombPBCAB_CUDA cab = CoulombPBCAB_CUDA(ions, elec);
+  CoulombPBCAB_CUDA cab(ions, elec);
 
   // Background charge term
   double consts = cab.evalConsts();
@@ -158,7 +158,7 @@ TEST_CASE("Coulomb PBC AB CUDA BCC H", "[hamiltonian][CUDA]")
   ParticleSetPool ptcl = ParticleSetPool(c);
 
 
-  CoulombPBCAB_CUDA cab = CoulombPBCAB_CUDA(ions, elec);
+  CoulombPBCAB_CUDA cab(ions, elec);
 
   // Background charge term
   double consts = cab.evalConsts();
@@ -227,7 +227,7 @@ TEST_CASE("Coulomb PBC A-A CUDA BCC H", "[hamiltonian][CUDA]")
   ions.updateLists_GPU();
 
 
-  CoulombPBCAA_CUDA caa = CoulombPBCAA_CUDA(ions, true);
+  CoulombPBCAA_CUDA caa(ions, true);
 
   // Background charge term
   double consts = caa.evalConsts();

src/QMCHamiltonians/tests/test_coulomb_pbcAA_ewald.cpp

@@ -53,10 +53,10 @@ TEST_CASE("Coulomb PBC A-A Ewald3D", "[hamiltonian]")
   int pMembersizeIdx                = ion_species.addAttribute("membersize");
   ion_species(pChargeIdx, pIdx)     = 1;
   ion_species(pMembersizeIdx, pIdx) = 1;
-  ions.Lattice = Lattice;
+  ions.Lattice                      = Lattice;
   ions.createSK();
 
-  LRCoulombSingleton::CoulombHandler = new EwaldHandler3D(ions);
+  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(ions);
   LRCoulombSingleton::CoulombHandler->initBreakup(ions);
 
 
@@ -68,8 +68,7 @@ TEST_CASE("Coulomb PBC A-A Ewald3D", "[hamiltonian]")
   double val = caa.evaluate(ions);
   REQUIRE(val == Approx(-1.418927)); // not validated
 
-  delete LRCoulombSingleton::CoulombHandler;
-  LRCoulombSingleton::CoulombHandler = 0;
+  LRCoulombSingleton::CoulombHandler.reset(nullptr);
 }
 
 TEST_CASE("Coulomb PBC A-A BCC H Ewald3D", "[hamiltonian]")
@@ -102,10 +101,10 @@ TEST_CASE("Coulomb PBC A-A BCC H Ewald3D", "[hamiltonian]")
   int pMembersizeIdx                = ion_species.addAttribute("membersize");
   ion_species(pChargeIdx, pIdx)     = 1;
   ion_species(pMembersizeIdx, pIdx) = 2;
-  ions.Lattice = Lattice;
+  ions.Lattice                      = Lattice;
   ions.createSK();
 
-  LRCoulombSingleton::CoulombHandler = new EwaldHandler3D(ions);
+  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(ions);
   LRCoulombSingleton::CoulombHandler->initBreakup(ions);
 
   CoulombPBCAA caa = CoulombPBCAA(ions, false);
@@ -117,8 +116,7 @@ TEST_CASE("Coulomb PBC A-A BCC H Ewald3D", "[hamiltonian]")
   double val = caa.evaluate(elec);
   REQUIRE(val == Approx(-0.963074)); // not validated
 
-  delete LRCoulombSingleton::CoulombHandler;
-  LRCoulombSingleton::CoulombHandler = 0;
+  LRCoulombSingleton::CoulombHandler.reset(nullptr);
 }
 
 TEST_CASE("Coulomb PBC A-A elec Ewald3D", "[hamiltonian]")
@@ -140,19 +138,19 @@ TEST_CASE("Coulomb PBC A-A elec Ewald3D", "[hamiltonian]")
   elec.R[0][1] = 0.5;
   elec.R[0][2] = 0.0;
 
-  SpeciesSet& tspecies         = elec.getSpeciesSet();
-  int upIdx                    = tspecies.addSpecies("u");
-  int chargeIdx                = tspecies.addAttribute("charge");
-  int massIdx                  = tspecies.addAttribute("mass");
-  int pMembersizeIdx           = tspecies.addAttribute("membersize");
-  tspecies(pMembersizeIdx, upIdx)   = 1;
-  tspecies(chargeIdx, upIdx)   = -1;
-  tspecies(massIdx, upIdx)     = 1.0;
+  SpeciesSet& tspecies            = elec.getSpeciesSet();
+  int upIdx                       = tspecies.addSpecies("u");
+  int chargeIdx                   = tspecies.addAttribute("charge");
+  int massIdx                     = tspecies.addAttribute("mass");
+  int pMembersizeIdx              = tspecies.addAttribute("membersize");
+  tspecies(pMembersizeIdx, upIdx) = 1;
+  tspecies(chargeIdx, upIdx)      = -1;
+  tspecies(massIdx, upIdx)        = 1.0;
 
   elec.createSK();
   elec.update();
 
-  LRCoulombSingleton::CoulombHandler = new EwaldHandler3D(elec);
+  LRCoulombSingleton::CoulombHandler = std::make_unique<EwaldHandler3D>(elec);
   LRCoulombSingleton::CoulombHandler->initBreakup(elec);
 
   CoulombPBCAA caa = CoulombPBCAA(elec, false);
@@ -164,8 +162,7 @@ TEST_CASE("Coulomb PBC A-A elec Ewald3D", "[hamiltonian]")
   double val = caa.evaluate(elec);
   REQUIRE(val == Approx(-1.418927)); // not validated
 
-  delete LRCoulombSingleton::CoulombHandler;
-  LRCoulombSingleton::CoulombHandler = 0;
+  LRCoulombSingleton::CoulombHandler.reset(nullptr);
 }
 
 

src/QMCHamiltonians/tests/test_coulomb_pbcAB.cpp

@@ -85,7 +85,7 @@ TEST_CASE("Coulomb PBC A-B", "[hamiltonian]")
   ParticleSetPool ptcl = ParticleSetPool(c);
 
 
-  CoulombPBCAB cab = CoulombPBCAB(ions, elec);
+  CoulombPBCAB cab(ions, elec);
 
   // Self energy plus Background charge term
   double consts = cab.evalConsts();
@@ -167,7 +167,7 @@ TEST_CASE("Coulomb PBC A-B BCC H", "[hamiltonian]")
   ParticleSetPool ptcl = ParticleSetPool(c);
 
 
-  CoulombPBCAB cab = CoulombPBCAB(ions, elec);
+  CoulombPBCAB cab(ions, elec);
 
   // Background charge term
   double consts = cab.evalConsts();