Changeset 2ba2ed for src/Potentials/Specifics/SaturationPotential.cpp

Timestamp:

Feb 25, 2013, 5:28:57 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:

dd8094

Parents:

b17e0f

git-author:

Frederik Heber <heber@…> (11/13/12 13:03:13)

git-committer:

Frederik Heber <heber@…> (02/25/13 17:28:57)

Message:

FIX: SaturationPotential's parameter derivatives seem to be working now.

• the problem seems to be the fitting itself is no good, maybe fitting Morse and Angle to CH4 at the same time just does not work ...

File:

• src/Potentials/Specifics/SaturationPotential.cpp (modified) (5 diffs)

src/Potentials/Specifics/SaturationPotential.cpp

@@ -45 +45 @@
 
 SaturationPotential::SaturationPotential(
+    const double _saturation_cutoff,
     boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
   energy_offset(0.),
   triplefunction(_triplefunction),
-  saturation_cutoff(2.5)
+  saturation_cutoff(_saturation_cutoff)
 {}
 
@@ -58 +59 @@
     const double _angle_equilibrium_distance,
     const double _all_energy_offset,
+    const double _saturation_cutoff,
     boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
   energy_offset(_all_energy_offset),
   triplefunction(_triplefunction),
-  saturation_cutoff(2.5)
+  saturation_cutoff(_saturation_cutoff)
 {
   parameters_t morse_params(morse.getParameterDimension());
@@ -148 +150 @@
       ++argiter) {
     const argument_t &r_ij = *argiter;
-    arguments_t args(1, r_ij);
-    if ((r_ij.indices.first == 0)) { // || (r_ij.indices.second == 0)) {
+    if ((r_ij.indices.first == 0)) { // first item must be the non-hydrogen
+      arguments_t args(1, r_ij);
+
+      // Morse contribution
       result += morse(args)[0];
       if (result != result)
         ELOG(1, "result is NAN.");
+
+      // Angle contribution
       std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
+      args.resize(3, r_ij);
       for (std::vector<arguments_t>::const_iterator iter = triples.begin();
           iter != triples.end(); ++iter) {
         ASSERT( iter->size() == 2,
-            "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
+            "SaturationPotential::function_derivative_c() - the triples result must contain exactly two distances.");
         const argument_t &r_ik = (*iter)[0];
         const argument_t &r_jk = (*iter)[1];
-        arguments_t args;
-        args.push_back(r_ij);
-        args.push_back(r_ik);
-        args.push_back(r_jk);
-        result += .5 * angle(args)[0];
+        args[1] = r_ik;
+        args[2] = r_jk;
+        result += .5*angle(args)[0];  // as we have all distances we get both jk and kj
         if (result != result)
           ELOG(1, "result is NAN.");
@@ -191 +196 @@
 {
   double result = 0.;
-  for(arguments_t::const_iterator argiter = arguments.begin();
-      argiter != arguments.end();
-      ++argiter) {
-    const argument_t &r_ij = *argiter;
-    arguments_t args(1, r_ij);
-    if ((r_ij.indices.first == 0)) { // || (r_ij.indices.second == 0)) {
-      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:
-        {
-          std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
-          for (std::vector<arguments_t>::const_iterator iter = triples.begin();
-              iter != triples.end(); ++iter) {
-            ASSERT( iter->size() == 2,
-                "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
-            const argument_t &r_ik = (*iter)[0];
-            const argument_t &r_jk = (*iter)[1];
-            arguments_t args;
-            args.push_back(r_ij);
-            args.push_back(r_ik);
-            args.push_back(r_jk);
-            switch (index) {
-              case angle_spring_constant:
-              {
-                result += .5*angle.parameter_derivative(args, PairPotential_Angle::spring_constant)[0];
-                break;
+  if (index == all_energy_offset) {
+    result = 1.;
+  } else {
+    for(arguments_t::const_iterator argiter = arguments.begin();
+        argiter != arguments.end();
+        ++argiter) {
+      const argument_t &r_ij = *argiter;
+      if ((r_ij.indices.first == 0)) { // first item must be the non-hydrogen
+        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:
+          {
+            args.resize(3, r_ij);
+            std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
+            for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+                iter != triples.end(); ++iter) {
+              ASSERT( iter->size() == 2,
+                  "SaturationPotential::parameter_derivative() - the triples result must contain exactly two distances.");
+              const argument_t &r_ik = (*iter)[0];
+              ASSERT( r_ik.indices.first == r_ij.indices.first,
+                  "SaturationPotential::parameter_derivative() - i not same in ij, ik.");
+              const argument_t &r_jk = (*iter)[1];
+              ASSERT( r_jk.indices.first == r_ij.indices.second,
+                  "SaturationPotential::parameter_derivative() - j not same in ij, jk.");
+              ASSERT( r_ik.indices.second == r_jk.indices.second,
+                  "SaturationPotential::parameter_derivative() - k not same in ik, jk.");
+              args[1] = r_ik;
+              args[2] = r_jk;
+              switch (index) {   // .5 due to we have all distances we get both jk and kj
+                case angle_spring_constant:
+                {
+                  result += .5*angle.parameter_derivative(args, PairPotential_Angle::spring_constant)[0];
+                  break;
+                }
+                case angle_equilibrium_distance:
+                {
+                  result += .5*angle.parameter_derivative(args, PairPotential_Angle::equilibrium_distance)[0];
+                  break;
+                }
+                default:
+                  break;
               }
-              case angle_equilibrium_distance:
-              {
-                result += .5*angle.parameter_derivative(args, PairPotential_Angle::equilibrium_distance)[0];
-                break;
-              }
-              default:
-                break;
             }
+            break;
           }
         }
@@ -245 +260 @@
     }
   }
-  if (index == all_energy_offset)
-    result = 1.;
   return SaturationPotential::results_t(1, result);
 }