Changeset a82602

Timestamp:

Dec 10, 2012, 10:10:59 AM (12 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

82cbe4

Parents:

5a5196

git-author:

Frederik Heber <heber@…> (09/12/12 15:23:13)

git-committer:

Frederik Heber <heber@…> (12/10/12 10:10:59)

Message:

InterfaceVMGJob now also calculates nuclei interaction energy.

• necessiated to Particle::commMPI().
• CommParticles() call in ImportRightHandSide().
• in ExportSolution() we just copied stuff from interface_particles.cpp from VMG project.
• returndata.e_long now contains calculated energy from there.

Location:

src/Jobs

Files:

• InterfaceVMGJob.cpp (modified) (5 diffs)
• VMGJob.cpp (modified) (5 diffs)
• VMGJob.hpp (modified) (2 diffs)

src/Jobs/InterfaceVMGJob.cpp

@@ -41 +41 @@
 #include "grid/grid.hpp"
 #include "grid/multigrid.hpp"
+#include "units/particle/comm_mpi_particle.hpp"
+#include "units/particle/interpolation.hpp"
+#include "units/particle/linked_cell_list.hpp"
 #include "mg.hpp"
 
@@ -98 +101 @@
 
   Grid& grid = multigrid(multigrid.MaxLevel());
+  grid.Clear();
   //grid.ClearBoundary(); // we don't have a boundary under periodic boundary conditions
 
@@ -114 +118 @@
    * Charge assignment on the grid
    */
-  Comm& comm = *MG::GetComm();
+  Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
   Grid& particle_grid = comm.GetParticleGrid();
-
   particle_grid.Clear();
+
+  // distribute particles
+  particles.clear();
+  comm.CommParticles(grid, particles);
 
   assert(particle_grid.Global().LocalSize().IsComponentwiseGreater(
@@ -180 +187 @@
 void InterfaceVMGJob::ExportSolution(Grid& grid)
 {
+  /// sample the obtained potential to evaluate with the electron charge density
+
   // grid now contains the sough-for potential
-  Comm& comm = *MG::GetComm();
+  //Comm& comm = *MG::GetComm();
+  Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
 
   const Index begin_local = grid.Global().LocalBegin() - grid.Local().HaloSize1();
@@ -207 +217 @@
   comm.PrintStringOnce("Grid potential sum: %e", potential_sum);
 
-  //Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());
-  returndata.e_long = potential_sum;
+  {
+    Grid::iterator grid_iter = grid.Iterators().Local().Begin();
+    comm.PrintStringOnce("Grid potential at (0,0,0): %e", grid.GetVal(*grid_iter));
+  }
+
+  //Particle::CommMPI& comm = *dynamic_cast<Particle::CommMPI*>(MG::GetComm());  returndata.e_long = potential_sum;
+
+  /// Calculate potential energy of nuclei
+
+  vmg_float e = 0.0;
+  vmg_float e_long = 0.0;
+  vmg_float e_self = 0.0;
+  vmg_float e_short_peak = 0.0;
+  vmg_float e_short_spline = 0.0;
+
+  Factory& factory = MG::GetFactory();
+
+  /*
+   * Get parameters and arrays
+   */
+  const vmg_int& near_field_cells = factory.GetObjectStorageVal<int>("PARTICLE_NEAR_FIELD_CELLS");
+  const vmg_int& interpolation_degree = factory.GetObjectStorageVal<int>("PARTICLE_INTERPOLATION_DEGREE");
+
+  Particle::Interpolation ip(interpolation_degree);
+
+  const vmg_float r_cut = near_field_cells * grid.Extent().MeshWidth().Max();
+
+  /*
+   * Copy potential values to a grid with sufficiently large halo size.
+   * This may be optimized in future.
+   * The parameters of this grid have been set in the import step.
+   */
+  Grid& particle_grid = comm.GetParticleGrid();
+
+  for (i.X()=0; i.X()<grid.Local().Size().X(); ++i.X())
+    for (i.Y()=0; i.Y()<grid.Local().Size().Y(); ++i.Y())
+      for (i.Z()=0; i.Z()<grid.Local().Size().Z(); ++i.Z())
+        particle_grid(i + particle_grid.Local().Begin()) = grid.GetVal(i + grid.Local().Begin());
+
+  comm.CommToGhosts(particle_grid);
+
+  /*
+   * Compute potentials
+   */
+  Particle::LinkedCellList lc(particles, near_field_cells, grid);
+  Particle::LinkedCellList::iterator p1, p2;
+  Grid::iterator iter;
+
+  comm.CommLCListToGhosts(lc);
+
+  for (int i=lc.Local().Begin().X(); i<lc.Local().End().X(); ++i)
+    for (int j=lc.Local().Begin().Y(); j<lc.Local().End().Y(); ++j)
+      for (int k=lc.Local().Begin().Z(); k<lc.Local().End().Z(); ++k) {
+
+  if (lc(i,j,k).size() > 0)
+    ip.ComputeCoefficients(particle_grid, Index(i,j,k) - lc.Local().Begin() + particle_grid.Local().Begin());
+
+  for (p1=lc(i,j,k).begin(); p1!=lc(i,j,k).end(); ++p1) {
+
+    // Interpolate long-range part of potential and electric field
+    ip.Evaluate(**p1);
+
+    // Subtract self-induced potential
+    (*p1)->Pot() -= (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
+
+    e_long += 0.5 * (*p1)->Charge() * ip.EvaluatePotentialLR(**p1);
+    e_self += 0.5 * (*p1)->Charge() * (*p1)->Charge() * spl.GetAntiDerivativeAtZero();
+
+    for (int dx=-1*near_field_cells; dx<=near_field_cells; ++dx)
+      for (int dy=-1*near_field_cells; dy<=near_field_cells; ++dy)
+        for (int dz=-1*near_field_cells; dz<=near_field_cells; ++dz) {
+
+    for (p2=lc(i+dx,j+dy,k+dz).begin(); p2!=lc(i+dx,j+dy,k+dz).end(); ++p2)
+
+      if (*p1 != *p2) {
+
+        const Vector dir = (*p1)->Pos() - (*p2)->Pos();
+        const vmg_float length = dir.Length();
+
+        if (length < r_cut) {
+
+          (*p1)->Pot() += (*p2)->Charge() / length * (1.0 + spl.EvaluatePotential(length));
+          (*p1)->Field() += (*p2)->Charge() * dir * spl.EvaluateField(length);
+
+          e_short_peak += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length;
+          e_short_spline += 0.5 * (*p1)->Charge() * (*p2)->Charge() / length * spl.EvaluatePotential(length);
+        }
+      }
+        }
+  }
+      }
+
+  /* Remove average force term */
+  Vector average_force = 0.0;
+  for (std::list<Particle::Particle>::const_iterator iter=particles.begin(); iter!=particles.end(); ++iter)
+    average_force += iter->Charge() * iter->Field();
+  const vmg_int& npl = MG::GetFactory().GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
+  const vmg_int num_particles_global = comm.GlobalSum(npl);
+  average_force /= num_particles_global;
+  comm.GlobalSumArray(average_force.vec(), 3);
+  for (std::list<Particle::Particle>::iterator iter=particles.begin(); iter!=particles.end(); ++iter)
+    iter->Field() -= average_force / iter->Charge();
+
+  comm.CommParticlesBack(particles);
+
+  vmg_float* q = factory.GetObjectStorageArray<vmg_float>("PARTICLE_CHARGE_ARRAY");
+  const vmg_int& num_particles_local = factory.GetObjectStorageVal<vmg_int>("PARTICLE_NUM_LOCAL");
+  const vmg_float* p = factory.GetObjectStorageArray<vmg_float>("PARTICLE_POTENTIAL_ARRAY");
+//  const vmg_float* f = factory.GetObjectStorageArray<vmg_float>("PARTICLE_FIELD_ARRAY");
+
+
+  e_long = comm.GlobalSumRoot(e_long);
+  e_short_peak = comm.GlobalSumRoot(e_short_peak);
+  e_short_spline = comm.GlobalSumRoot(e_short_spline);
+  e_self = comm.GlobalSumRoot(e_self);
+
+  for (int j=0; j<num_particles_local; ++j)
+    e += 0.5 * p[j] * q[j];
+  e = comm.GlobalSumRoot(e);
+
+  comm.PrintStringOnce("E_long:         %e", e_long);
+  comm.PrintStringOnce("E_short_peak:   %e", e_short_peak);
+  comm.PrintStringOnce("E_short_spline: %e", e_short_spline);
+  comm.PrintStringOnce("E_self:         %e", e_self);
+  comm.PrintStringOnce("E_total:        %e", e);
+  comm.PrintStringOnce("E_total*:       %e", e_long + e_short_peak + e_short_spline - e_self);
+
+  returndata.e_long = e;
 }

src/Jobs/VMGJob.cpp

@@ -63 +63 @@
 //#include "solver/solver_regular.hpp"
 #include "solver/solver_singular.hpp"
-//#include "units/particle/comm_mpi_particle.hpp"
+#include "units/particle/comm_mpi_particle.hpp"
 
 #include "CodePatterns/MemDebug.hpp"
@@ -71 +71 @@
 #include "LinearAlgebra/defs.hpp"
 #include "Jobs/InterfaceVMGJob.hpp"
+
+#include "CodePatterns/Assert.hpp"
 
 using namespace VMG;
@@ -85 +87 @@
   particle_charges(_particle_charges),
   near_field_cells(_near_field_cells),
-  returndata(static_cast<const SamplingGridProperties &>(_density_grid))
+  returndata(static_cast<const SamplingGridProperties &>(_density_grid)),
+  num_particles(0),
+  f(NULL),
+  x(NULL),
+  p(NULL),
+  q(NULL)
 {}
 
 VMGJob::VMGJob() :
   FragmentJob(JobId::IllegalJob),
-  near_field_cells(5)
+  near_field_cells(5),
+  num_particles(0),
+  f(NULL),
+  x(NULL),
+  p(NULL),
+  q(NULL)
 {}
 
 VMGJob::~VMGJob()
-{}
+{
+  delete[] f;
+  delete[] x;
+  delete[] p;
+  delete[] q;
+}
+
+void VMGJob::InitVMGArrays()
+{
+  num_particles = particle_charges.size();
+  f = new double[num_particles*3];
+  x = new double[num_particles*3];
+  p = new double[num_particles];
+  q = new double[num_particles];
+
+  {
+    size_t index=0;
+    for (std::vector< std::vector< double > >::const_iterator iter = particle_positions.begin();
+        iter != particle_positions.end(); ++iter) {
+      for (std::vector< double >::const_iterator positer = (*iter).begin();
+          positer != (*iter).end(); ++positer) {
+        f[index] = 0.;
+        x[index++] = *positer;
+      }
+    }
+    ASSERT( index == num_particles*3,
+        "VMGJob::VMGJob() - too many particles or components in particle_positions.");
+  }
+  {
+    size_t index = 0;
+    for (std::vector< double >::const_iterator iter = particle_charges.begin();
+        iter != particle_charges.end(); ++iter) {
+      p[index] = 0.;
+      q[index++] = *iter;
+    }
+    ASSERT( index == num_particles,
+        "VMGJob::VMGJob() - too many charges in particle_charges.");
+  }
+}
+
 
 FragmentResult::ptr VMGJob::Work()
@@ -137 +188 @@
    */
 #ifdef HAVE_MPI
-  new CommMPI(boundary, new DomainDecompositionMPI());
+  new Particle::CommMPI(boundary, new DomainDecompositionMPI());
 #else
   new CommSerial(boundary);
@@ -169 +220 @@
   new ObjectStorage<int>("PARTICLE_INTERPOLATION_DEGREE", 3);
 
+  // first initialize f,x,p,q,... from STL vectors
+  InitVMGArrays();
+  // then initialize as objects with VMG's storage
+  new ObjectStorage<vmg_float*>("PARTICLE_POS_ARRAY", x);
+  new ObjectStorage<vmg_float*>("PARTICLE_CHARGE_ARRAY", q);
+  new ObjectStorage<vmg_float*>("PARTICLE_POTENTIAL_ARRAY", p);
+  new ObjectStorage<vmg_float*>("PARTICLE_FIELD_ARRAY", f);
+  new ObjectStorage<vmg_int>("PARTICLE_NUM_LOCAL", num_particles);
+
   /*
    * Post init

src/Jobs/VMGJob.hpp

@@ -53 +53 @@
   void InitVMG();
 
+  void InitVMGArrays();
+
 private:
   //!> sampled density required as input
@@ -66 +68 @@
   //!> temporary instance to hold return data
   VMGData returndata;
+
+  //!> number of particles
+  size_t num_particles;
+  //!> forces array
+  double *f;
+  //!> position array
+  double *x;
+  //!> potential array
+  double *p;
+  //!> charge array
+  double *q;
 
 private: