Merge pull request #2605 from jtkrogel/ed_fix2

Fix energy density for periodic grids

docs/hamiltonianobservable.rst

@@ -1178,23 +1178,26 @@ In QMCPACK, the energy density can be accumulated on piecewise uniform 3D grids
 
 attributes:
 
-  +------------------------+--------------+----------------------+-------------+---------------------------+
-  | **Name**               | **Datatype** | **Values**           | **Default** | **Description**           |
-  +========================+==============+======================+=============+===========================+
-  | ``type``:math:`^r`     | text         | **EnergyDensity**    |             | Must be EnergyDensity     |
-  +------------------------+--------------+----------------------+-------------+---------------------------+
-  | ``name``:math:`^r`     | text         | *anything*           |             | Unique name for estimator |
-  +------------------------+--------------+----------------------+-------------+---------------------------+
-  | ``dynamic``:math:`^r`  | text         | ``particleset.name`` |             | Identify electrons        |
-  +------------------------+--------------+----------------------+-------------+---------------------------+
-  | ``static``:math:`^o`   | text         | ``particleset.name`` |             | Identify ions             |
-  +------------------------+--------------+----------------------+-------------+---------------------------+
+  +----------------------------+--------------+----------------------+-------------+----------------------------------------------+
+  | **Name**                   | **Datatype** | **Values**           | **Default** | **Description**                              |
+  +============================+==============+======================+=============+==============================================+
+  | ``type``:math:`^r`         | text         | **EnergyDensity**    |             | Must be EnergyDensity                        |
+  +----------------------------+--------------+----------------------+-------------+----------------------------------------------+
+  | ``name``:math:`^r`         | text         | *anything*           |             | Unique name for estimator                    |
+  +----------------------------+--------------+----------------------+-------------+----------------------------------------------+
+  | ``dynamic``:math:`^r`      | text         | ``particleset.name`` |             | Identify electrons                           |
+  +----------------------------+--------------+----------------------+-------------+----------------------------------------------+
+  | ``static``:math:`^o`       | text         | ``particleset.name`` |             | Identify ions                                |
+  +----------------------------+--------------+----------------------+-------------+----------------------------------------------+
+  | ``ion_points``:math:`^o`   | text         | yes/no               |  no         | Separate ion energy density onto point field |
+  +----------------------------+--------------+----------------------+-------------+----------------------------------------------+
 
 Additional information:
 
 -  ``name:`` Must be unique. A dataset with blocked statistical data for
    the energy density will appear in the ``stat.h5`` files labeled as
    ``name``.
+- **Important:** in order for the estimator to work, a traces XML input element (<traces array="yes" write="no"/>) must appear following the <qmcsystem/> element and prior to any <qmc/> element.
 
 .. code-block::
   :caption: Energy density estimator accumulated on a :math:`20 \times  10 \times 10` grid over the simulation cell.

nexus/lib/qmcpack_input.py

@@ -2158,9 +2158,10 @@ class localenergy(QIxml):
 
 class energydensity(QIxml):
     tag = 'estimator'
-    attributes = ['type','name','dynamic','static']
-    elements   = ['reference_points','spacegrid']
-    identifier = 'name'
+    attributes  = ['type','name','dynamic','static','ion_points']
+    elements    = ['reference_points','spacegrid']
+    identifier  = 'name'
+    write_types = obj(ion_points=yesno)
 #end class energydensity
 
 class reference_points(QIxml):
@@ -2802,7 +2803,7 @@ class gen(QIxml):
     type='spindensity',name='SpinDensity'
     )
 skall.defaults.set(
-    type='skall',name='skall'
+    type='skall',name='skall',source='ion0',target='e',hdf5=True
     )
 force.defaults.set(
     type='Force',name='force'
@@ -4922,6 +4923,8 @@ def generate_hamiltonian(name         = 'h0',
                     est = chiesa(name='KEcorr',type='chiesa',source=ename,psi=wfname)
                 elif estname=='localenergy':
                     est = localenergy(name='LocalEnergy')
+                elif estname=='skall':
+                    est = skall(name='SkAll',type='skall',source=iname,target=ename,hdf5=True)
                 elif estname=='energydensity':
                     est = energydensity(
                         type='EnergyDensity',name='EDvoronoi',dynamic=ename,static=iname,

nexus/lib/qmcpack_quantity_analyzers.py

@@ -301,7 +301,7 @@ def load_data_local(self,data=None):
         exclude = self.info.exclude
         self.data = QAHDFdata()
         for var in list(data.keys()):
-            if not var in exclude and not str(var)[0]=='_':
+            if not var in exclude and not str(var)[0]=='_' and not 'skall' in var.lower():
                 self.data[var] = data[var]
                 del data[var]
             #end if
@@ -522,6 +522,46 @@ def analyze_local(self):
             P = P
             )
 
+        # convert ion point data, if present
+        if 'ions' in self:
+            ions = QAobject()
+            ions.D  = QAobject()
+            ions.T  = QAobject()
+            ions.V  = QAobject()
+            ions.E  = QAobject()
+            ions.P  = QAobject()
+
+            value = self.ions.value.transpose()[...,nbe:]
+
+            mean,var,error,kappa = simstats(value)
+            ions.D.mean   = mean[iD]
+            ions.D.error  = error[iD]
+            ions.T.mean   = mean[iT]
+            ions.T.error  = error[iT]
+            ions.V.mean   = mean[iV]
+            ions.V.error  = error[iV]
+
+            E  = value[iT,:]+value[iV,:]
+            mean,var,error,kappa = simstats(E)
+            ions.E.mean  = mean
+            ions.E.error = error
+
+            P  = 2./3.*value[iT,:]+1./3.*value[iV,:]
+            mean,var,error,kappa = simstats(P)
+            ions.P.mean  = mean
+            ions.P.error = error
+
+            ions.data = obj(
+                D = value[iD,:],
+                T = value[iT,:],
+                V = value[iV,:],
+                E = E,
+                P = P
+                )
+
+            self.ions = ions
+        #end if
+
         return
     #end def analyze_local
 
@@ -2830,7 +2870,7 @@ def initialize(self): #like qmcpack SpaceGridBase.initialize
                     ndu_per_interval[iaxis][idom] = ndu_int[i]
                     idom+=1
                 #end 
-          #end       
+            #end       
         #end 
 
         axinv = inv(axes)

nexus/lib/simulation.py

@@ -1405,6 +1405,18 @@ def execute(self,run_command=None):
         #end if
     #end def execute
 
+
+    def show_input(self,exit=True):
+        print()
+        print(80*'=')
+        print('Input file for simulation "{}"\nDirectory: {}'.format(self.identifier,self.locdir))
+        print(80*'-')
+        print(self.input.write())
+        print(80*'=')
+        if exit:
+            exit_call()
+        #end if
+    #end def show_input
 #end class Simulation
 
 

src/QMCHamiltonians/EnergyDensityEstimator.cpp

@@ -93,7 +93,13 @@ bool EnergyDensityEstimator::put(xmlNodePtr cur)
   R.resize(nparticles);
   EDValues.resize(nparticles, nEDValues);
   if (ion_points)
+  {
     EDIonValues.resize(nions,nEDValues);
+    Rion.resize(nions,DIM);
+    for (int i=0; i<nions; i++)
+      for (int d=0; d<DIM; d++)
+        Rion(i,d) = Pstatic->R[i][d];
+  }
   particles_outside.resize(nparticles);
   fill(particles_outside.begin(), particles_outside.end(), true);
   //read xml element contents
@@ -127,6 +133,7 @@ bool EnergyDensityEstimator::put(xmlNodePtr cur)
     stop               = stop || !ref_succeeded;
   }
   //initialize grids or other cell partitions
+  bool periodic = Pdynamic->Lattice.SuperCellEnum != SUPERCELL_OPEN;
   bool grid_succeeded;
   element     = cur->children;
   int nvalues = (int)nEDValues;
@@ -141,11 +148,11 @@ bool EnergyDensityEstimator::put(xmlNodePtr cur)
       if (Pstatic)
       {
         set_ptcl();
-        grid_succeeded = sg->put(element, ref.points, Rptcl, Zptcl, Pdynamic->getTotalNum(), false);
+        grid_succeeded = sg->put(element, ref.points, Rptcl, Zptcl, Pdynamic->getTotalNum(), periodic, false);
         unset_ptcl();
       }
       else
-        grid_succeeded = sg->put(element, ref.points, false);
+        grid_succeeded = sg->put(element, ref.points, periodic, false);
       stop = stop || !grid_succeeded;
       ++i;
     }
@@ -265,7 +272,7 @@ EnergyDensityEstimator::Return_t EnergyDensityEstimator::evaluate(ParticleSet& P
         p++;
       }
     }
-    if (Pstatic)
+    if (Pstatic && !ion_points)
     {
       const ParticlePos_t& Rs = Pstatic->R;
       for (int i = 0; i < Rs.size(); i++)
@@ -297,13 +304,21 @@ EnergyDensityEstimator::Return_t EnergyDensityEstimator::evaluate(ParticleSet& P
       Vs_trace->combine();
       const ParticleSet& Ps            = *Pstatic;
       const std::vector<TraceReal>& Vs = Vs_trace->sample;
-      for (int i = 0; i < Ps.getTotalNum(); i++)
-      {
-        EDValues(p, W) = w;
-        EDValues(p, T) = 0.0;
-        EDValues(p, V) = w * Vs[i];
-        p++;
-      }
+      if (!ion_points)
+        for (int i = 0; i < Ps.getTotalNum(); i++)
+        {
+          EDValues(p, W) = w;
+          EDValues(p, T) = 0.0;
+          EDValues(p, V) = w * Vs[i];
+          p++;
+        }
+      else
+        for (int i = 0; i < Ps.getTotalNum(); i++)
+        {
+          EDIonValues(i, W) = w;
+          EDIonValues(i, T) = 0.0;
+          EDIonValues(i, V) = w * Vs[i];
+        }
     }
     //Accumulate energy density in spacegrids
     const DistanceTableData& dtab(P.getDistTable(dtable_index));
@@ -326,6 +341,16 @@ EnergyDensityEstimator::Return_t EnergyDensityEstimator::evaluate(ParticleSet& P
         }
       }
     }
+    if (ion_points)
+    {
+      // Accumulate energy density for ions at a point field
+      bi = ion_buffer_offset;
+      for (int i =0; i<nions; i++)
+        for (v=0; v<(int)nEDValues; v++, bi++)
+        {
+          P.Collectables[bi] += EDIonValues(i,v);
+        }
+    }
     nsamples++;
 #if defined(ENERGYDENSITY_CHECK)
     int thread = omp_get_thread_num();
@@ -345,6 +370,13 @@ EnergyDensityEstimator::Return_t EnergyDensityEstimator::evaluate(ParticleSet& P
         Tsum += EDValues(p, T);
         Vsum += EDValues(p, V);
       }
+      if (ion_points)
+        for (int i = 0; i < nions; i++)
+        {
+          Dsum += EDIonValues(i, W);
+          Tsum += EDIonValues(i, T);
+          Vsum += EDIonValues(i, V);
+        }
       Esum           = Tsum + Vsum;
       static int cnt = 0;
       //app_log()<<"eval ED Dsum"<<cnt<<" "<<Dsum<< std::endl;
@@ -363,6 +395,12 @@ EnergyDensityEstimator::Return_t EnergyDensityEstimator::evaluate(ParticleSet& P
       }
       for (int v = 0; v < nvals; v++)
         edvals[v] += P.Collectables[outside_buffer_offset + v];
+      if (ion_points)
+      {
+        for (int i = 0; i<nions; i++)
+          for (int v = 0; v < nvals; v++)
+            edvals[v] += P.Collectables[ion_buffer_offset + i*nvals + v];
+      }
       //app_log()<<"eval ES Dsum"<<cnt<<" "<<edvals[W]<< std::endl;
       app_log() << thread << " eval ES " << cnt << " " << edvals[T] << " " << edvals[V] << " " << edvals[T] + edvals[V]
                 << std::endl;
@@ -433,6 +471,13 @@ void EnergyDensityEstimator::addObservables(PropertySetType& plist, BufferType&
   {
     spacegrids[i]->allocate_buffer_space(collectables);
   }
+  if (ion_points)
+  {
+    ion_buffer_offset = collectables.size();
+    nvalues = nions*((int)nEDValues);
+    std::vector<RealType> tmp2(nvalues);
+    collectables.add(tmp2.begin(), tmp2.end());
+  }
 }
 
 
@@ -446,6 +491,11 @@ void EnergyDensityEstimator::registerCollectables(std::vector<observable_helper*
   int nspacegrids = spacegrids.size();
   oh->addProperty(const_cast<int&>(nspacegrids), "nspacegrids");
   oh->addProperty(const_cast<int&>(nsamples), "nsamples");
+  if (ion_points)
+  {
+    oh->addProperty(const_cast<int&>(nions), "nions");
+    oh->addProperty(const_cast<Matrix<RealType>&>(Rion), "ion_positions");
+  }
   h5desc.push_back(oh);
   ref.save(h5desc, g);
   oh = new observable_helper("outside");
@@ -459,6 +509,16 @@ void EnergyDensityEstimator::registerCollectables(std::vector<observable_helper*
     SpaceGrid& sg = *spacegrids[i];
     sg.registerCollectables(h5desc, g, i);
   }
+  if (ion_points)
+  {
+    oh = new observable_helper("ions");
+    std::vector<int> ng2(2);
+    ng2[0] = nions;
+    ng2[1] = (int)nEDValues;
+    oh->set_dimensions(ng2, ion_buffer_offset);
+    oh->open(g);
+    h5desc.push_back(oh);
+  }
 }
 
 void EnergyDensityEstimator::setObservables(PropertySetType& plist)

src/QMCHamiltonians/EnergyDensityEstimator.h

@@ -60,8 +60,10 @@ class EnergyDensityEstimator : public OperatorBase, public PtclOnLatticeTraits
   ParticleSet* get_particleset(std::string& psname);
   int dtable_index;
   int nparticles;
-  int nions;
   bool ion_points;
+  int nions;
+  int ion_buffer_offset;
+  Matrix<RealType> Rion;
   //collection of points from which to build spacegrid origin and axes
   ReferencePoints ref;
   //EnergyDenstity quantities