[WIP] Anisotropic Ewald code

src/LongRange/LRCoulombSingleton.cpp

@@ -119,10 +119,10 @@ LRCoulombSingleton::LRHandlerType* LRCoulombSingleton::getHandler(ParticleSet& r
     }
     else //if(ref.LRBox.SuperCellEnum == SUPERCELL_BULK)
     {
-      app_log() << "\n  Creating CoulombHandler with the optimal breakup. " << std::endl;
-      CoulombHandler = new LRHandlerTemp<CoulombFunctor<mRealType>, LPQHIBasis>(ref);
-      //  app_log() << "\n  Creating CoulombHandler with the Ewald3D breakup. " << std::endl;
-      //  CoulombHandler= new EwaldHandler3D(ref);
+      //app_log() << "\n  Creating CoulombHandler with the optimal breakup. " << std::endl;
+      //CoulombHandler = new LRHandlerTemp<CoulombFunctor<mRealType>, LPQHIBasis>(ref);
+        app_log() << "\n  Creating CoulombHandler with the Ewald3D breakup. " << std::endl;
+        CoulombHandler= new EwaldHandler3D(ref);
       //  CoulombHandler = new LRHandlerSRCoulomb<CoulombFunctor<mRealType>, LPQHISRCoulombBasis>(ref);
     }
 //        else if(ref.LRBox.SuperCellEnum == SUPERCELL_SLAB)

src/LongRange/StructFact.cpp

@@ -76,6 +76,8 @@ void StructFact::UpdateAllPart(ParticleSet& P) { FillRhok(P); }
  */
 void StructFact::FillRhok(ParticleSet& P)
 {
+  ScopedTimer rhoktimer(TimerManager.createTimer("StructFact::FillRhok"));
+
   int npart = P.getTotalNum();
 #if defined(USE_REAL_STRUCT_FACTOR)
   rhok_r = 0.0;

src/QMCHamiltonians/CMakeLists.txt

@@ -66,6 +66,7 @@ IF(OHMMS_DIM MATCHES 3)
     ECPComponentBuilder.1.cpp
     ECPComponentBuilder.2.cpp
     ECPComponentBuilder_L2.cpp
+    EwaldRef.cpp
     )
 
 

src/QMCHamiltonians/CoulombPBCAA.cpp

@@ -14,14 +14,15 @@
 //////////////////////////////////////////////////////////////////////////////////////
 
 
-#include "EwaldRef.h"
 #include "QMCHamiltonians/CoulombPBCAA.h"
 #include "Particle/DistanceTableData.h"
 #include "Utilities/ProgressReportEngine.h"
 #include <numeric>
 
 namespace qmcplusplus
 {
+
+
 CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
     : ForceBase(ref, ref),
       AA(0),
@@ -48,26 +49,203 @@ CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
     ref.turnOnPerParticleSK();
     update_source(ref);
   }
+
+
+  // Setup anistropic ewald
+  ewaldref::RealMat A;
+  ewaldref::ChargeArray Q;
+
+  A = Ps.Lattice.R;
+
+  Q.resize(NumCenters);
+  for(int i=0;i<NumCenters;++i)
+    Q[i] = Zat[i];
+
+
+
+  //RealType qmcpack_kappa = AA->LR_kc;
+  RealType qmcpack_sigma = std::sqrt(AA->LR_kc / (2.0 * AA->LR_rc));
+  RealType qmcpack_kappa = 1./(std::sqrt(2.0)*qmcpack_sigma);
+  ewaldtools::AnisotropicEwald ewald_qp(A,Q,1e-10,qmcpack_kappa);
+
+  //ewald.initialize(A,Q);
+  ewald.initialize(A,Q,1e-10,qmcpack_kappa);
+
+
+  //ewaldref::IntVec nmax = 7;
+  //ewald.setupOpt(nmax);
+
   if (!is_active)
   {
     update_source(ref);
+
 
-    ewaldref::RealMat A;
     ewaldref::PosArray R;
-    ewaldref::ChargeArray Q;
-
-    A = Ps.Lattice.R;
-
     R.resize(NumCenters);
-    Q.resize(NumCenters);
-    for (int i = 0; i < NumCenters; ++i)
-    {
+    for(int i=0;i<NumCenters;++i)
       R[i] = Ps.R[i];
-      Q[i] = Zat[i];
-    }
 
     RealType Vii_ref        = ewaldref::ewaldEnergy(A, R, Q);
     RealType Vdiff_per_atom = std::abs(Value - Vii_ref) / NumCenters;
+
+
+
+
+    //app_log()<<std::setprecision(14);
+    app_log()<<std::setprecision(16);
+    app_log()<<std::endl;
+    app_log()<<std::endl;
+    app_log()<<"   ewald ref i-i   energy: "<<Vii_ref<<std::endl;
+    app_log()<<"   qmcpack   i-i   energy: "<<Value<<std::endl;
+    app_log()<<std::endl;
+    app_log()<<std::endl;
+
+    const DistanceTableData& dt(Ps.getDistTable(d_aa_ID));
+
+    //RealType E_ref = ewald.ewaldEnergy(R);
+    //RealType E_fix = ewald.fixedGridEwaldEnergy(R);
+    RealType c     = ewald.ewaldEnergyConst();
+    RealType sr    = ewald.ewaldEnergySR(Ps.R);
+    RealType lr    = ewald.ewaldEnergyLR(Ps.R);
+
+
+    //app_log()<<"   ewald adaptive  energy: "<<E_ref<<std::endl;
+    //app_log()<<"   ewald fixed     energy: "<<E_fix<<std::endl;
+    app_log()<<"   ewald component energy: "<<sr+lr+c<<std::endl;
+
+    app_log()<<std::endl;
+    app_log()<<std::endl;
+
+    app_log()<<"   sr ref : "<< sr <<std::endl;
+    app_log()<<"   sr opt : "<< ewald.ewaldEnergySROpt(dt) <<"  "<<ewald.getKappa() <<std::endl;
+
+    app_log()<<std::endl;
+    app_log()<<std::endl;
+
+    // make print out of SR accuracy vs kappa
+    //ewald.printSRErrorTable(R,dt);
+
+    // LR accuracy vs cutoff
+    //app_log()<<"   lr ref : "<< lr <<std::endl;
+
+    //ewald.findOptKGrid(R);
+
+
+
+
+    //auto g = 2*ewald.nkmax+1;
+    //size_t nkpoints_ref = g[0]*g[1]*g[2];
+
+    //ewald.setupOpt(R,dt);
+    //
+    //app_log()<<"   sr ref : "<< sr_ref <<std::endl;
+    //app_log()<<"   sr opt : "<< ewald.ewaldEnergySROpt(dt) <<"  "<<ewald.getKappa() <<std::endl;
+    //
+    //
+    //app_log()<<std::endl;
+    //app_log()<<std::endl;
+    //
+    //app_log()<<"   lr ref : "<< lr_ref<<" "<< nkpoints_ref <<std::endl;
+    //app_log()<<"   lr opt : "<< ewald.ewaldEnergyLROpt(R) <<"  "<< ewald.getNkpoints() <<std::endl;
+    //
+    //app_log()<<"   ewald optimized energy: "<<ewald.ewaldEnergyOpt(R,dt)<<std::endl;
+
+
+
+    app_log()<<std::endl;
+    app_log()<<std::endl;
+
+    ewald.printErrorAndTimeVSTolerance(R,dt);
+
+
+
+    //RealType eE   = ewald.ewaldEnergy(R);
+    //RealType eEf  = ewald.fixedGridEwaldEnergy(R);
+    ////RealType eEt  = ewald.ewaldEnergy(Ps.R);
+    ////RealType eEft = ewald.fixedGridEwaldEnergy(Ps.R);
+    //
+    //RealType c  = ewald.ewaldEnergyConst();
+    //RealType sr = ewald.ewaldEnergySR(Ps.R);
+    //RealType lr = ewald.ewaldEnergyLR(Ps.R);
+    //
+    //const DistanceTableData& dt(Ps.getDistTable(d_aa_ID));
+    ////RealType srdt = ewald.ewaldEnergySRDT(dt);
+    ////RealType lrdt = ewald.ewaldEnergyLRDT(dt);
+    //
+    //
+    ////ewald_qp.madelungEnergy0();
+    //RealType eE_qp = ewald_qp.ewaldEnergy(Ps.R);
+    ////RealType c_qp   = ewald_qp.ewaldEnergyConst();
+    ////RealType sr_qp  = ewald_qp.ewaldEnergySRDT(dt);
+    ////RealType lr_qp  = ewald_qp.ewaldEnergyLRDT(dt);
+    ////RealType sr0_qp = ewald_qp.ewaldEnergySR0DT(dt);
+    //
+    //RealType c_qmc = myConst;
+    //RealType sr_qmc = evalSR(Ps);
+    //RealType lr_qmc = evalLR(Ps);
+    //RealType E_qmc  = c_qmc + sr_qmc + lr_qmc;
+    //
+    //
+    //app_log()<<std::setprecision(14);
+    //app_log()<<std::endl;
+    //app_log()<<"   adaptive aniso ewald energy: "<<eE<<std::endl;
+    //app_log()<<"   fixed aniso ewald energy   : "<<eEf<<std::endl;
+    //
+    ////app_log()<<"   adaptive aniso ewald energy: "<<eEt<<std::endl;
+    ////app_log()<<"   fixed aniso ewald energy   : "<<eEft<<std::endl;
+    //app_log()<<"   ewald energy c+sr+lr       : "<<c+sr+lr<<std::endl;
+    ////app_log()<<"   ewald energy c+sr+lr dt    : "<<c+srdt+lrdt<<std::endl;
+    //app_log()<<"   aniso ewald nmax           : "<<ewald.getNmax()<<std::endl;
+    //app_log()<<std::setprecision(14);
+    ////app_log()<<"   qmcpack constant : "<<c_qmc<<std::endl;
+    ////app_log()<<"   ewald   constant : "<<c_qp<<std::endl;
+    ////app_log()<<"   qmcpack SR  : "<<sr_qmc<<std::endl;
+    ////app_log()<<"   ewald   SR  : "<<sr_qp<<std::endl;
+    ////app_log()<<"   ewald   SR0 : "<<sr0_qp<<std::endl;
+    ////app_log()<<"   qmcpack LR  : "<<lr_qmc<<std::endl;
+    ////app_log()<<"   ewald   LR  : "<<lr_qp<<std::endl;
+    ////app_log()<<"   qmcpack tot : "<<c_qmc+sr_qmc+lr_qmc<<std::endl;
+    ////app_log()<<"   ewald   tot : "<< c_qp+sr_qp+lr_qp <<std::endl;
+    ////app_log()<<"   ewald   tot0: "<< c_qp+sr0_qp+lr_qp <<std::endl;
+    ////app_log()<<std::endl;
+    ////app_log()<<"   SR error    : "<<sr_qmc-sr_qp<<std::endl;
+    ////app_log()<<"   tot error   : "<<c_qmc+sr_qmc+lr_qmc-c_qp-sr_qp-lr_qp<<std::endl;
+    ////app_log()<<std::endl;
+    //app_log()<<"   ewald_qp aniso energy      : "<<eE<<std::endl;
+    //app_log()<<"   aniso ewald nmax           : "<<ewald_qp.getNmax()<<std::endl;
+    //app_log()<<std::setprecision(14);
+    //
+    //
+    ////ewaldtools::IntVec nmax = 20;
+    ////ewald.setupOpt(nmax);
+    //
+    ////auto& ewald_opt = ewald_qp;
+    //auto& ewald_opt = ewald;
+    //
+    //RealType c_ref  = ewald_opt.ewaldEnergyConst();
+    //RealType sr_ref = ewald_opt.ewaldEnergySR(Ps.R);
+    //RealType lr_ref = ewald_opt.ewaldEnergyLR(Ps.R);
+    ////ewald_opt.setupOpt();
+    //ewald_opt.setupOpt(nmax);
+    //RealType c_opt  = ewald_opt.ewaldEnergyConst();
+    //RealType sr_opt = ewald_opt.ewaldEnergySROpt(dt);
+    //RealType lr_opt = ewald_opt.ewaldEnergyLROpt(Ps.R);
+    //app_log()<<"   ewald opt energy           : "<<c_opt+sr_opt+lr_opt<<std::endl;
+    //app_log()<<"   ewald opt energy2          : "<<ewald.ewaldEnergyOpt(Ps.R,dt)<<std::endl;
+    //app_log()<<"   qmcpack   energy           : "<<E_qmc<<std::endl;
+    //app_log()<<"   ewald     LR  : "<<lr_ref<<std::endl;
+    //app_log()<<"   ewald opt LR  : "<<lr_opt<<std::endl;
+    //app_log()<<"   ewald     SR  : "<<sr_ref<<std::endl;
+    //app_log()<<"   ewald opt SR  : "<<sr_opt<<std::endl;
+    //app_log()<<"   qmcpack   SR  : "<<sr_qmc<<std::endl;
+    //
+    //app_log()<<std::endl;
+    //app_log()<<std::endl;
+
+
+    APP_ABORT("explore")
+
+
     app_log() << "Checking ion-ion Ewald energy against reference..." << std::endl;
     if (Vdiff_per_atom > Ps.Lattice.LR_tol)
     {
@@ -95,6 +273,7 @@ CoulombPBCAA::CoulombPBCAA(ParticleSet& ref, bool active, bool computeForces)
     {
       app_log() << "  Check passed." << std::endl;
     }
+
   }
   prefix = "F_AA";
   app_log() << "  Maximum K shell " << AA->MaxKshell << std::endl;
@@ -163,14 +342,96 @@ void CoulombPBCAA::delete_particle_quantities()
 
 CoulombPBCAA::Return_t CoulombPBCAA::evaluate(ParticleSet& P)
 {
+  ScopedTimer evaltimer(TimerManager.createTimer("CoulombPBCAA_eval_total"));
+
+  const DistanceTableData& dt(P.getDistTable(d_aa_ID));
+
   if (is_active)
   {
-#if !defined(REMOVE_TRACEMANAGER)
-    if (streaming_particles)
-      Value = evaluate_sp(P);
-    else
-#endif
-      Value = evalLR(P) + evalSR(P) + myConst;
+//    {
+//    ScopedTimer evaltimer(TimerManager.createTimer("QMCPACK eval"));
+//
+//#if !defined(REMOVE_TRACEMANAGER)
+//    if (streaming_particles)
+//      Value = evaluate_sp(P);
+//    else
+//#endif
+//      Value = evalLR(P) + evalSR(P) + myConst;
+//    }
+//
+//    RealType eval;
+//    {
+//      ScopedTimer evaltimer(TimerManager.createTimer("Aniso eval"));
+//      //eval = ewald.fixedGridEwaldEnergy(P.R);
+//      eval = ewald.ewaldEnergyOpt(P.R,dt);
+//    }
+
+
+    {
+      ScopedTimer evaltimer(TimerManager.createTimer("QMCPACK eval"));
+      RealType sr_qmc;
+      RealType lr_qmc;
+      {
+        ScopedTimer evaltimer(TimerManager.createTimer("SR qmc"));
+        sr_qmc = evalSR(P);
+      }
+      {
+        ScopedTimer evaltimer(TimerManager.createTimer("LR qmc"));
+        lr_qmc = evalLR(P);
+      }
+      Value = lr_qmc + sr_qmc + myConst;
+    }
+
+    RealType E_opt;
+    RealType c_opt;
+    RealType sr_opt;
+    RealType lr_opt;
+    {
+      ScopedTimer evaltimer(TimerManager.createTimer("Aniso eval"));
+      c_opt = ewald.ewaldEnergyConst();
+      {
+        ScopedTimer evaltimer(TimerManager.createTimer("SR opt"));
+        sr_opt = ewald.ewaldEnergySROpt(dt);
+      }
+      {
+        ScopedTimer evaltimer(TimerManager.createTimer("LR opt"));
+        lr_opt = ewald.ewaldEnergyLROpt(P.R);
+      }
+    }
+    E_opt = lr_opt + sr_opt + c_opt;
+
+
+    RealType sr_opt_orig = ewald.ewaldEnergySROpt_ref(dt);
+    RealType lr_opt_orig = ewald.ewaldEnergyLROpt_ref(P.R);
+    RealType lr_opt_qmc  = ewald.ewaldEnergyLROpt_qmcpack(P.R);
+
+    RealType E_opt2 = c_opt + sr_opt_orig + lr_opt_orig;
+    RealType E_opt3 = c_opt + sr_opt_orig + lr_opt_qmc;
+
+    //RealType E_ref  = ewald.ewaldEnergy(P.R);
+
+    //RealType c_ref  = ewald.ewaldEnergyConst();
+    //RealType sr_ref = ewald.ewaldEnergySR(P.R);
+    //RealType lr_ref = ewald.ewaldEnergyLR(P.R);
+    //
+    //RealType sr_opt = ewald.ewaldEnergySROpt(dt);
+    //RealType lr_opt = ewald.ewaldEnergyLROpt(P.R);
+
+    app_log()<<std::endl;
+    app_log()<<"  QMCPACK    value: "<<Value<<std::endl;
+    //app_log()<<"  Aniso ref  value: "<<E_ref<<std::endl;
+    //app_log()<<"  Aniso ref  value: "<<c_ref+sr_ref+lr_ref<<std::endl;
+    app_log()<<"  Aniso opt  value: "<<E_opt<<std::endl;
+    app_log()<<"  Aniso opt2 value: "<<E_opt2<<std::endl;
+    app_log()<<"  Aniso opt3 value: "<<E_opt3<<std::endl;
+    //app_log()<<"  QMCPACK       SR: "<<evalSR(P)<<std::endl;
+    //app_log()<<"  Aniso     ref SR: "<<sr_ref<<std::endl;
+    //app_log()<<"  Aniso     opt SR: "<<sr_opt<<std::endl;
+    //app_log()<<"  QMCPACK       LR: "<<evalLR(P)+myConst<<std::endl;
+    //app_log()<<"  Aniso     ref LR: "<<lr_ref+c_ref<<std::endl;
+    //app_log()<<"  Aniso     opt LR: "<<lr_opt+c_ref<<std::endl;
+    app_log()<<std::endl;
+
   }
   return Value;
 }

src/QMCHamiltonians/CoulombPBCAA.h

@@ -20,6 +20,8 @@
 #include "QMCHamiltonians/ForceBase.h"
 #include "LongRange/LRCoulombSingleton.h"
 #include "Particle/DistanceTableData.h"
+//#include "QMCHamiltonians/EwaldRef.h"
+#include "QMCHamiltonians/EwaldTools.h"
 
 namespace qmcplusplus
 {
@@ -36,6 +38,8 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
   typedef LRCoulombSingleton::RadFunctorType RadFunctorType;
   typedef LRHandlerType::mRealType mRealType;
 
+  ewaldtools::AnisotropicEwald ewald;
+
   // energy-optimized
   LRHandlerType* AA;
   GridType* myGrid;