Changeset ffc368 for src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

Timestamp:

Dec 19, 2012, 3:25:53 PM (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, 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:

990a62

Parents:

e7579e

git-author:

Frederik Heber <heber@…> (10/04/12 22:01:47)

git-committer:

Frederik Heber <heber@…> (12/19/12 15:25:53)

Message:

ManyBodyPotential_Tersoff is now correctly implemented.

• changed implementation to match with [Tersoff, '89] with incorporates configurations consisting of different elements and which is also the one implemented in tremolo (for whose ptential parameters we aim for).
• we have unit test of operator() that checks against a set of values obtained from a brief run with tremolo of a configuration consisting of 5 carbon atoms: the results are the same by each function.

File:

• src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp (modified) (11 diffs)

src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

@@ -43 +43 @@
 
 #include "CodePatterns/Assert.hpp"
+//#include "CodePatterns/Info.hpp"
 
 #include "Potentials/helpers.hpp"
@@ -54 +55 @@
 
 ManyBodyPotential_Tersoff::ManyBodyPotential_Tersoff(
-    const double &_cutoff_offset,
-    const double &_cutoff_halfwidth,
+    const double &_R,
+    const double &_S,
     const double &_A,
     const double &_B,
-    const double &_lambda1,
-    const double &_lambda2,
+    const double &_lambda,
+    const double &_mu,
     const double &_lambda3,
     const double &_alpha,
     const double &_beta,
+    const double &_chi,
+    const double &_omega,
     const double &_n,
     const double &_c,
@@ -71 +74 @@
   triplefunction(_triplefunction)
 {
-  params[cutoff_offset] = _cutoff_offset;
-  params[cutoff_halfwidth] = _cutoff_halfwidth;
+//  Info info(__func__);
+  params[R] = _R;
+  params[S] = _S;
   params[A] = _A;
   params[B] = _B;
-  params[lambda1] = _lambda1;
-  params[lambda2] = _lambda2;
+  params[lambda] = _lambda;
+  params[mu] = _mu;
   params[lambda3] = _lambda3;
   params[alpha] = _alpha;
   params[beta] = _beta;
+  params[chi] = _chi;
+  params[omega] = _omega;
   params[n] = _n;
   params[c] = _c;
@@ -91 +97 @@
     ) const
 {
-  const double &distance = arguments[0].distance;
-  const double cutoff = function_cutoff(distance);
-  const double result = (cutoff == 0.) ? 0. : cutoff * (
-      function_prefactor(
-          alpha,
-          boost::bind(&ManyBodyPotential_Tersoff::function_eta,
-              boost::cref(*this),
-              boost::cref(arguments[0])))
-      * function_smoother(
-          distance,
-          A,
-          lambda1)
-      +
-      function_prefactor(
-          beta,
-          boost::bind(&ManyBodyPotential_Tersoff::function_zeta,
-              boost::cref(*this),
-              boost::cref(arguments[0])))
-      * function_smoother(
-          distance,
-          -B,
-          lambda2)
-  );
+//  Info info(__func__);
+  const argument_t &r_ij = arguments[0];
+  const double cutoff = function_cutoff(r_ij.distance);
+  const double result = (cutoff == 0.) ?
+      0. :
+      cutoff * (
+          function_prefactor(
+              params[alpha],
+              function_eta(r_ij))
+          * function_smoother(
+              params[A],
+              params[lambda],
+              r_ij.distance)
+          +
+          function_prefactor(
+              params[beta],
+              function_zeta(r_ij))
+          * function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+      );
+//  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
   return std::vector<result_t>(1, result);
+}
+
+ManyBodyPotential_Tersoff::derivative_components_t
+ManyBodyPotential_Tersoff::derivative(
+    const arguments_t &arguments
+    ) const
+{
+//  Info info(__func__);
+  return ManyBodyPotential_Tersoff::derivative_components_t();
+}
+
+ManyBodyPotential_Tersoff::results_t
+ManyBodyPotential_Tersoff::parameter_derivative(
+    const arguments_t &arguments,
+    const size_t index
+    ) const
+{
+//  Info info(__func__);
+//  ASSERT( arguments.size() == 1,
+//      "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
+  const argument_t &r_ij = arguments[0];
+  switch (index) {
+    case R:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case S:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case A:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case B:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case lambda:
+    {
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. :
+          -r_ij.distance * cutoff * params[lambda] * (
+          function_prefactor(
+              params[alpha],
+              function_eta(r_ij))
+          * function_smoother(
+              params[A],
+              params[lambda],
+              r_ij.distance)
+          );
+      return results_t(1, result);
+      break;
+    }
+    case mu:
+    {
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. :
+          -r_ij.distance * cutoff * params[mu] *(
+          function_prefactor(
+              params[beta],
+              function_zeta(r_ij))
+          * function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+      );
+      return results_t(1, result);
+      break;
+    }
+    case lambda3:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case alpha:
+    {
+      const double temp =
+          pow(params[alpha]*function_eta(r_ij), params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) || (params[alpha] == 0. )?
+          0. :
+          function_smoother(
+              -params[A],
+              params[lambda],
+              r_ij.distance)
+          * (-.5) * params[alpha] * (temp/params[alpha])
+          / (1. + temp)
+          ;
+      return results_t(1, result);
+      break;
+    }
+    case beta:
+    {
+      const double temp =
+          pow(params[beta]*function_zeta(r_ij), params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) || (params[beta] == 0. )?
+          0. : cutoff *
+          function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+          * (-.5) * params[beta] * (temp/params[beta])
+          / (1. + temp)
+          ;
+      return results_t(1, result);
+      break;
+    }
+    case n:
+    {
+      const double temp =
+          pow(params[beta]*function_zeta(r_ij), params[n]);
+      const double cutoff = function_cutoff(r_ij.distance);
+      const double result = (cutoff == 0.) ?
+          0. : cutoff *
+          function_smoother(
+              -params[B],
+              params[mu],
+              r_ij.distance)
+          * params[B]
+          * ( log(1.+temp)/(params[n]*params[n]) - temp
+              * (log(function_zeta(r_ij)) + log(params[beta]))
+              /(params[n]*(1.+temp)))
+          ;
+      return results_t(1, result);
+      break;
+    }
+    case c:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case d:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    case h:
+    {
+      const double result = 0.;
+      return results_t(1, result);
+      break;
+    }
+    default:
+      break;
+  }
+  return results_t(1, 0.);
 }
 
@@ -122 +289 @@
   ) const
 {
-  const double offset = (distance - cutoff_offset);
-  if (offset < - cutoff_halfwidth)
-    return 1.;
-  else if (offset > cutoff_halfwidth)
-    return 0.;
+//  Info info(__func__);
+  double result = 0.;
+  if (distance < params[R])
+    result = 1.;
+  else if (distance > params[S])
+    result = 0.;
   else {
-    return (0.5 - 0.5 * sin( .5 * M_PI * offset/cutoff_halfwidth));
+    result = (0.5 + 0.5 * cos( M_PI * (distance - params[R])/(params[S]-params[R])));
   }
+//  LOG(2, "DEBUG: function_cutoff(" << distance << ") = " << result);
+  return result;
 }
 
@@ -135 +305 @@
 ManyBodyPotential_Tersoff::function_prefactor(
     const double &alpha,
-    boost::function<result_t()> etafunction
+    const double &eta
   ) const
 {
-  return pow(
-      (1. + Helpers::pow(alpha * etafunction(), n)),
-      -1./(2.*n));
+//  Info info(__func__);
+  const double result = params[chi] * pow(
+      (1. + pow(alpha * eta, params[n])),
+      -1./(2.*params[n]));
+//  LOG(2, "DEBUG: function_prefactor(" << alpha << "," << eta << ") = " << result);
+  return result;
+}
+
+ManyBodyPotential_Tersoff::result_t
+ManyBodyPotential_Tersoff::function_smoother(
+    const double &prefactor,
+    const double &lambda,
+    const double &distance
+    ) const
+{
+//  Info info(__func__);
+  const double result = prefactor * exp(-lambda * distance);
+//  LOG(2, "DEBUG: function_smoother(" << prefactor << "," << lambda << "," << distance << ") = " << result);
+  return result;
 }
 
@@ -148 +334 @@
   ) const
 {
+//  Info info(__func__);
   result_t result = 0.;
 
   // get all triples within the cutoff
-  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
+  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
   for (std::vector<arguments_t>::const_iterator iter = triples.begin();
       iter != triples.end(); ++iter) {
@@ -158 +345 @@
     const argument_t &r_ik = (*iter)[0];
     result += function_cutoff(r_ik.distance)
-        * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
+        * exp( Helpers::pow(params[lambda3] * (r_ij.distance - r_ik.distance) ,3));
   }
 
+//  LOG(2, "DEBUG: function_eta(" << r_ij.distance << ") = " << result);
   return result;
 }
@@ -169 +357 @@
   ) const
 {
+//  Info info(__func__);
   result_t result = 0.;
 
   // get all triples within the cutoff
-  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
+  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
   for (std::vector<arguments_t>::const_iterator iter = triples.begin();
       iter != triples.end(); ++iter) {
@@ -181 +370 @@
     result +=
         function_cutoff(r_ik.distance)
+        * params[omega]
         * function_angle(r_ij.distance, r_ik.distance, r_jk.distance)
-        * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
+        * exp( Helpers::pow(params[lambda3] * (r_ij.distance - r_ik.distance) ,3));
   }
 
+//  LOG(2, "DEBUG: function_zeta(" << r_ij.distance << ") = " << result);
   return result;
 }
@@ -195 +386 @@
   ) const
 {
+//  Info info(__func__);
   const double angle = Helpers::pow(r_ij,2) + Helpers::pow(r_ik,2) - Helpers::pow(r_jk,2);
-  const double divisor = r_ij * r_ik;
-  const double result =
-      1.
-      + (Helpers::pow(c, 2)/Helpers::pow(d, 2))
-      - Helpers::pow(c, 2)/(Helpers::pow(d, 2) +
-          Helpers::pow(h - cos(angle/divisor),2));
-  return result;
-}
+  const double divisor = 2.* r_ij * r_ik;
+//  LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
+  if (divisor != 0.) {
+    const double result =
+        1.
+        + (Helpers::pow(params[c]/params[d], 2))
+        - Helpers::pow(params[c], 2)/(Helpers::pow(params[d], 2) +
+            Helpers::pow(params[h] - angle/divisor,2));
+
+//    LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
+    return result;
+  } else
+    return 0.;
+}
+