Changeset 0dad5f for src

Timestamp:

Feb 27, 2013, 12:39:03 PM (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:

e36ba2

Parents:

775dd1a

git-author:

Frederik Heber <heber@…> (12/07/12 15:01:07)

git-committer:

Frederik Heber <heber@…> (02/27/13 12:39:03)

Message:

SaturationPotential now has an internal triplefunction.

• SaturationPotential simply to get all possible distances and from these we extract triplets of matching distances via an internal triplefunction as input to PairPotential_Angle.
• rewrote operator() and parameter_derivative(), gradients are remain ok.

Location:

src

Files:

• LevMartester.cpp (modified) (1 diff)
• Potentials/Specifics/SaturationPotential.cpp (modified) (5 diffs)
• Potentials/Specifics/SaturationPotential.hpp (modified) (3 diffs)

src/LevMartester.cpp

@@ -369 +369 @@
 
       SaturationPotential saturation(types);
-      saturation.setTriplefunction(2.5, triplefunction);
+      // saturation has internal triplefunction now
+//      saturation.setTriplefunction(2.5, triplefunction);
       FunctionModel &model = saturation;
       FunctionApproximation approximator(TersoffData, model); // CH4 has 5 atoms, hence 5*4/2 distances

src/Potentials/Specifics/SaturationPotential.cpp

@@ -176 +176 @@
 }
 
+std::vector<FunctionModel::arguments_t>
+triplefunction(
+    const argument_t &argument,
+    const FunctionModel::arguments_t& args)
+{
+  const size_t firstindex = argument.indices.first;
+  const size_t secondindex = argument.indices.second;
+//  LOG(2, "DEBUG: first index is " << firstindex << ", second index is " << secondindex << ".");
+
+  // place all arguments that share either index into a lookup map
+  typedef std::map< size_t, FunctionModel::arguments_t::const_iterator > IndexLookup_t;
+  IndexLookup_t LookuptoFirst;
+  IndexLookup_t LookuptoSecond;
+  for (FunctionModel::arguments_t::const_iterator iter = args.begin();
+      iter != args.end();
+      ++iter) {
+    if (((*iter).indices.first == argument.indices.first)
+        && ((*iter).indices.second == argument.indices.second))
+      continue;
+    if (firstindex == (*iter).indices.first) {
+      LookuptoFirst.insert( std::make_pair( (*iter).indices.second, iter) );
+    }
+    else if (firstindex == (*iter).indices.second) {
+      LookuptoFirst.insert( std::make_pair( (*iter).indices.first, iter) );
+    }
+    if (secondindex == (*iter).indices.first) {
+      LookuptoSecond.insert( std::make_pair( (*iter).indices.second, iter) );
+    }
+    else if (secondindex == (*iter).indices.second) {
+      LookuptoSecond.insert( std::make_pair((*iter).indices.first, iter) );
+    }
+  }
+//  {
+//    std::stringstream lookupstream;
+//    for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
+//        iter != LookuptoFirst.end();
+//        ++iter) {
+//      lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
+//    }
+//    LOG(2, "DEBUG: LookupToFirst is " << lookupstream.str() << ".");
+//  }
+//  {
+//    std::stringstream lookupstream;
+//    for (IndexLookup_t::const_iterator iter = LookuptoSecond.begin();
+//        iter != LookuptoSecond.end();
+//        ++iter) {
+//      lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
+//    }
+//    LOG(2, "DEBUG: LookuptoSecond is " << lookupstream.str() << ".");
+//  }
+
+  // now go through the first lookup as the second argument and pick the
+  // corresponding third argument by the matching index
+  std::vector<FunctionModel::arguments_t> results;
+  for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
+      iter != LookuptoFirst.end();
+      ++iter) {
+    IndexLookup_t::const_iterator otheriter = LookuptoSecond.find(iter->first);
+    ASSERT( otheriter != LookuptoSecond.end(),
+        "triplefunction() - cannot find index "+toString(iter->first)
+        +" in LookupToSecond");
+    FunctionModel::arguments_t result(1, argument);
+    result.reserve(3);
+    result.push_back(*(iter->second));
+    result.push_back(*(otheriter->second));
+    results.push_back(result);
+  }
+
+  return results;
+}
+
 SaturationPotential::results_t
 SaturationPotential::operator()(
@@ -197 +268 @@
       if (result != result)
         ELOG(1, "result is NAN.");
-    }
-  }
-  {
-    // Angle contribution
-    const double tmp = angle(arguments)[0];  // as we have all distances we get both jk and kj
-//    LOG(2, "DEBUG: angle yields " << tmp << " for << " << arguments << ".");
-    result += tmp;
-    if (result != result)
-      ELOG(1, "result is NAN.");
+
+      // Angle contribution
+      {
+        typedef std::vector<FunctionModel::arguments_t> tripleargs_t;
+        tripleargs_t tripleargs =
+            triplefunction(r_ij, arguments);
+        for (tripleargs_t::const_iterator iter = tripleargs.begin();
+            iter != tripleargs.end();
+            ++iter) {
+          FunctionModel::arguments_t tempargs =
+              Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
+          const double tmp = angle(tempargs)[0];  // as we have all distances we get both jk and kj
+//          LOG(2, "DEBUG: angle yields " << tmp << " for << " << arguments << ".");
+          result += tmp;
+          if (result != result)
+            ELOG(1, "result is NAN.");
+        }
+      }
+    }
   }
   return std::vector<result_t>(1, energy_offset + result);
@@ -228 +309 @@
 {
   double result = 0.;
-  switch (index) {
-    case all_energy_offset:
-    {
-      result = 1.;
-      break;
-    }
-    case morse_spring_constant:
-    case morse_equilibrium_distance:
-    case morse_dissociation_energy:
-    {
-      const ParticleTypes_t &morse_types = morse.getParticleTypes();
-      for(arguments_t::const_iterator argiter = arguments.begin();
-          argiter != arguments.end();
-          ++argiter) {
-        const argument_t &r_ij = *argiter;
-        if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
-            || ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
-          arguments_t args(1, r_ij);
-          switch (index) {
-            case morse_spring_constant:
-              result += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
-              break;
-            case morse_equilibrium_distance:
-              result += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
-              break;
-            case morse_dissociation_energy:
-              result += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
-              break;
-            default:
-              ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
-              break;
+  if (index == all_energy_offset) {
+    result = 1.;
+  } else {
+    const ParticleTypes_t &morse_types = morse.getParticleTypes();
+    for(arguments_t::const_iterator argiter = arguments.begin();
+        argiter != arguments.end();
+        ++argiter) {
+      const argument_t &r_ij = *argiter;
+      if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
+          || ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
+        arguments_t args(1, r_ij);
+
+        switch (index) {
+          case morse_spring_constant:
+            result += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
+            break;
+          case morse_equilibrium_distance:
+            result += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
+            break;
+          case morse_dissociation_energy:
+            result += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
+            break;
+          case angle_spring_constant:
+          case angle_equilibrium_distance:
+          {
+            typedef std::vector<FunctionModel::arguments_t> tripleargs_t;
+            tripleargs_t tripleargs =
+                triplefunction(r_ij, arguments);
+            for (tripleargs_t::const_iterator iter = tripleargs.begin();
+                iter != tripleargs.end();
+                ++iter) {
+              FunctionModel::arguments_t tempargs =
+                  Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
+              if (index == angle_spring_constant)
+                result += angle.parameter_derivative(tempargs, PairPotential_Angle::spring_constant)[0];
+              else if (index == angle_equilibrium_distance)
+                result += angle.parameter_derivative(tempargs, PairPotential_Angle::equilibrium_distance)[0];
+//              LOG(2, "DEBUG: angle yields " << tmp << " for << " << arguments << ".");
+//              result += tmp;
+              if (result != result)
+                ELOG(1, "result is NAN.");
+            }
           }
+            break;
+          default:
+            ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
+            break;
         }
       }
-      break;
-    }
-    case angle_spring_constant:
-    {
-      result = angle.parameter_derivative(arguments, PairPotential_Angle::spring_constant)[0];
-      break;
-    }
-    case angle_equilibrium_distance:
-    {
-      result = angle.parameter_derivative(arguments, PairPotential_Angle::equilibrium_distance)[0];
-      break;
-    }
-    default:
-      ELOG(1, "SaturationPotential::parameter_derivative() - index " << index << " invalid.");
-      break;
+    }
   }
   return SaturationPotential::results_t(1, result);
@@ -333 +415 @@
     // Finally, we also need to arrange them in correct order
     // (for PairPotentiale_Angle).
-    charges_t firstpair(2, boost::cref(charges[0]));
+    //charges_t firstpair(2, boost::cref(charges[0]));
+    // only that saturation potential never has its middle element twice!
+    // hence, we skip the firstpair but keep the code for later generalization
     charges_t secondpair(2, boost::cref(charges[1]));
     const charges_t &thirdpair = charges;
     returnfunction =
-        boost::bind(&Extractors::reorderArgumentsByParticleTypes,
+//        boost::bind(&Extractors::reorderArgumentsByParticleTypes,
           boost::bind(&Extractors::combineArguments,
-            boost::bind(&Extractors::combineArguments,
-              boost::bind(&Extractors::gatherAllDistancesFromFragment,
-                  boost::bind(&Fragment::getPositions, _1),
-                  boost::bind(&Fragment::getCharges, _1),
-                  firstpair,  // no crefs here as are temporaries!
-                  _2),
+//            boost::bind(&Extractors::combineArguments,
+//              boost::bind(&Extractors::gatherAllDistancesFromFragment,
+//                  boost::bind(&Fragment::getPositions, _1),
+//                  boost::bind(&Fragment::getCharges, _1),
+//                  firstpair,  // no crefs here as are temporaries!
+//                  _2),
               boost::bind(&Extractors::gatherAllDistancesFromFragment,
                   boost::bind(&Fragment::getPositions, _1),
                   boost::bind(&Fragment::getCharges, _1),
                   secondpair,  // no crefs here as are temporaries!
-                  _2)
-            ),
+                  _2),
+//            ),
             boost::bind(&Extractors::gatherAllDistancesFromFragment,
                 boost::bind(&Fragment::getPositions, _1),
@@ -356 +440 @@
                 boost::cref(thirdpair), // only the last one is no temporary
                 _2)
-          ),
-          boost::cref(angle.getParticleTypes())
-        );
+          );
+//          boost::cref(angle.getParticleTypes())
+//        );
 }
   return returnfunction;

src/Potentials/Specifics/SaturationPotential.hpp

@@ -90 +90 @@
   /** Sets the magic triple function that we use for getting angle distances.
    *
+   * \note does nothing as we have an internal triplefunction
+   *
    * \param _saturation_cutoff cutoff for the triplefunction
    * @param _triplefunction function that returns a list of triples (i.e. the
@@ -100 +102 @@
       )
   {
-    saturation_cutoff = _saturation_cutoff;
-    triplefunction = _triplefunction;
   }
 
@@ -207 +207 @@
   double energy_offset;
 
-  //!> bound function that obtains the triples for the internal coordinationb summation.
-  boost::function< std::vector< arguments_t >(const argument_t &, const double)> triplefunction;
-  double saturation_cutoff;
-
   //!> static definitions of the parameter name for this potential
   static const ParameterNames_t ParameterNames;