Changeset a82602 for src/Jobs/InterfaceVMGJob.cpp

Timestamp:

Dec 10, 2012, 10:10:59 AM (13 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, Candidate_v1.7.0, Candidate_v1.7.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.

File:

• src/Jobs/InterfaceVMGJob.cpp (modified) (5 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;
 }