Changeset ffc368


Ignore:
Timestamp:
Dec 19, 2012, 3:25:53 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:
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.
Location:
src/Potentials/Specifics
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

    re7579e rffc368  
    4343
    4444#include "CodePatterns/Assert.hpp"
     45//#include "CodePatterns/Info.hpp"
    4546
    4647#include "Potentials/helpers.hpp"
    … …  
    5455
    5556ManyBodyPotential_Tersoff::ManyBodyPotential_Tersoff(
    56     const double &_cutoff_offset,
    57     const double &_cutoff_halfwidth,
     57    const double &_R,
     58    const double &_S,
    5859    const double &_A,
    5960    const double &_B,
    60     const double &_lambda1,
    61     const double &_lambda2,
     61    const double &_lambda,
     62    const double &_mu,
    6263    const double &_lambda3,
    6364    const double &_alpha,
    6465    const double &_beta,
     66    const double &_chi,
     67    const double &_omega,
    6568    const double &_n,
    6669    const double &_c,
    … …  
    7174  triplefunction(_triplefunction)
    7275{
    73   params[cutoff_offset] = _cutoff_offset;
    74   params[cutoff_halfwidth] = _cutoff_halfwidth;
     76//  Info info(__func__);
     77  params[R] = _R;
     78  params[S] = _S;
    7579  params[A] = _A;
    7680  params[B] = _B;
    77   params[lambda1] = _lambda1;
    78   params[lambda2] = _lambda2;
     81  params[lambda] = _lambda;
     82  params[mu] = _mu;
    7983  params[lambda3] = _lambda3;
    8084  params[alpha] = _alpha;
    8185  params[beta] = _beta;
     86  params[chi] = _chi;
     87  params[omega] = _omega;
    8288  params[n] = _n;
    8389  params[c] = _c;
    … …  
    9197    ) const
    9298{
    93   const double &distance = arguments[0].distance;
    94   const double cutoff = function_cutoff(distance);
    95   const double result = (cutoff == 0.) ? 0. : cutoff * (
    96       function_prefactor(
    97           alpha,
    98           boost::bind(&ManyBodyPotential_Tersoff::function_eta,
    99               boost::cref(*this),
    100               boost::cref(arguments[0])))
    101       * function_smoother(
    102           distance,
    103           A,
    104           lambda1)
    105       +
    106       function_prefactor(
    107           beta,
    108           boost::bind(&ManyBodyPotential_Tersoff::function_zeta,
    109               boost::cref(*this),
    110               boost::cref(arguments[0])))
    111       * function_smoother(
    112           distance,
    113           -B,
    114           lambda2)
    115   );
     99//  Info info(__func__);
     100  const argument_t &r_ij = arguments[0];
     101  const double cutoff = function_cutoff(r_ij.distance);
     102  const double result = (cutoff == 0.) ?
     103      0. :
     104      cutoff * (
     105          function_prefactor(
     106              params[alpha],
     107              function_eta(r_ij))
     108          * function_smoother(
     109              params[A],
     110              params[lambda],
     111              r_ij.distance)
     112          +
     113          function_prefactor(
     114              params[beta],
     115              function_zeta(r_ij))
     116          * function_smoother(
     117              -params[B],
     118              params[mu],
     119              r_ij.distance)
     120      );
     121//  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
    116122  return std::vector<result_t>(1, result);
     123}
     124
     125ManyBodyPotential_Tersoff::derivative_components_t
     126ManyBodyPotential_Tersoff::derivative(
     127    const arguments_t &arguments
     128    ) const
     129{
     130//  Info info(__func__);
     131  return ManyBodyPotential_Tersoff::derivative_components_t();
     132}
     133
     134ManyBodyPotential_Tersoff::results_t
     135ManyBodyPotential_Tersoff::parameter_derivative(
     136    const arguments_t &arguments,
     137    const size_t index
     138    ) const
     139{
     140//  Info info(__func__);
     141//  ASSERT( arguments.size() == 1,
     142//      "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
     143  const argument_t &r_ij = arguments[0];
     144  switch (index) {
     145    case R:
     146    {
     147      const double result = 0.;
     148      return results_t(1, result);
     149      break;
     150    }
     151    case S:
     152    {
     153      const double result = 0.;
     154      return results_t(1, result);
     155      break;
     156    }
     157    case A:
     158    {
     159      const double result = 0.;
     160      return results_t(1, result);
     161      break;
     162    }
     163    case B:
     164    {
     165      const double result = 0.;
     166      return results_t(1, result);
     167      break;
     168    }
     169    case lambda:
     170    {
     171      const double cutoff = function_cutoff(r_ij.distance);
     172      const double result = (cutoff == 0.) ?
     173          0. :
     174          -r_ij.distance * cutoff * params[lambda] * (
     175          function_prefactor(
     176              params[alpha],
     177              function_eta(r_ij))
     178          * function_smoother(
     179              params[A],
     180              params[lambda],
     181              r_ij.distance)
     182          );
     183      return results_t(1, result);
     184      break;
     185    }
     186    case mu:
     187    {
     188      const double cutoff = function_cutoff(r_ij.distance);
     189      const double result = (cutoff == 0.) ?
     190          0. :
     191          -r_ij.distance * cutoff * params[mu] *(
     192          function_prefactor(
     193              params[beta],
     194              function_zeta(r_ij))
     195          * function_smoother(
     196              -params[B],
     197              params[mu],
     198              r_ij.distance)
     199      );
     200      return results_t(1, result);
     201      break;
     202    }
     203    case lambda3:
     204    {
     205      const double result = 0.;
     206      return results_t(1, result);
     207      break;
     208    }
     209    case alpha:
     210    {
     211      const double temp =
     212          pow(params[alpha]*function_eta(r_ij), params[n]);
     213      const double cutoff = function_cutoff(r_ij.distance);
     214      const double result = (cutoff == 0.) || (params[alpha] == 0. )?
     215          0. :
     216          function_smoother(
     217              -params[A],
     218              params[lambda],
     219              r_ij.distance)
     220          * (-.5) * params[alpha] * (temp/params[alpha])
     221          / (1. + temp)
     222          ;
     223      return results_t(1, result);
     224      break;
     225    }
     226    case beta:
     227    {
     228      const double temp =
     229          pow(params[beta]*function_zeta(r_ij), params[n]);
     230      const double cutoff = function_cutoff(r_ij.distance);
     231      const double result = (cutoff == 0.) || (params[beta] == 0. )?
     232          0. : cutoff *
     233          function_smoother(
     234              -params[B],
     235              params[mu],
     236              r_ij.distance)
     237          * (-.5) * params[beta] * (temp/params[beta])
     238          / (1. + temp)
     239          ;
     240      return results_t(1, result);
     241      break;
     242    }
     243    case n:
     244    {
     245      const double temp =
     246          pow(params[beta]*function_zeta(r_ij), params[n]);
     247      const double cutoff = function_cutoff(r_ij.distance);
     248      const double result = (cutoff == 0.) ?
     249          0. : cutoff *
     250          function_smoother(
     251              -params[B],
     252              params[mu],
     253              r_ij.distance)
     254          * params[B]
     255          * ( log(1.+temp)/(params[n]*params[n]) - temp
     256              * (log(function_zeta(r_ij)) + log(params[beta]))
     257              /(params[n]*(1.+temp)))
     258          ;
     259      return results_t(1, result);
     260      break;
     261    }
     262    case c:
     263    {
     264      const double result = 0.;
     265      return results_t(1, result);
     266      break;
     267    }
     268    case d:
     269    {
     270      const double result = 0.;
     271      return results_t(1, result);
     272      break;
     273    }
     274    case h:
     275    {
     276      const double result = 0.;
     277      return results_t(1, result);
     278      break;
     279    }
     280    default:
     281      break;
     282  }
     283  return results_t(1, 0.);
    117284}
    118285
    … …  
    122289  ) const
    123290{
    124   const double offset = (distance - cutoff_offset);
    125   if (offset < - cutoff_halfwidth)
    126     return 1.;
    127   else if (offset > cutoff_halfwidth)
    128     return 0.;
     291//  Info info(__func__);
     292  double result = 0.;
     293  if (distance < params[R])
     294    result = 1.;
     295  else if (distance > params[S])
     296    result = 0.;
    129297  else {
    130     return (0.5 - 0.5 * sin( .5 * M_PI * offset/cutoff_halfwidth));
     298    result = (0.5 + 0.5 * cos( M_PI * (distance - params[R])/(params[S]-params[R])));
    131299  }
     300//  LOG(2, "DEBUG: function_cutoff(" << distance << ") = " << result);
     301  return result;
    132302}
    133303
    … …  
    135305ManyBodyPotential_Tersoff::function_prefactor(
    136306    const double &alpha,
    137     boost::function<result_t()> etafunction
     307    const double &eta
    138308  ) const
    139309{
    140   return pow(
    141       (1. + Helpers::pow(alpha * etafunction(), n)),
    142       -1./(2.*n));
     310//  Info info(__func__);
     311  const double result = params[chi] * pow(
     312      (1. + pow(alpha * eta, params[n])),
     313      -1./(2.*params[n]));
     314//  LOG(2, "DEBUG: function_prefactor(" << alpha << "," << eta << ") = " << result);
     315  return result;
     316}
     317
     318ManyBodyPotential_Tersoff::result_t
     319ManyBodyPotential_Tersoff::function_smoother(
     320    const double &prefactor,
     321    const double &lambda,
     322    const double &distance
     323    ) const
     324{
     325//  Info info(__func__);
     326  const double result = prefactor * exp(-lambda * distance);
     327//  LOG(2, "DEBUG: function_smoother(" << prefactor << "," << lambda << "," << distance << ") = " << result);
     328  return result;
    143329}
    144330
    … …  
    148334  ) const
    149335{
     336//  Info info(__func__);
    150337  result_t result = 0.;
    151338
    152339  // get all triples within the cutoff
    153   std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
     340  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
    154341  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
    155342      iter != triples.end(); ++iter) {
    … …  
    158345    const argument_t &r_ik = (*iter)[0];
    159346    result += function_cutoff(r_ik.distance)
    160         * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
     347        * exp( Helpers::pow(params[lambda3] * (r_ij.distance - r_ik.distance) ,3));
    161348  }
    162349
     350//  LOG(2, "DEBUG: function_eta(" << r_ij.distance << ") = " << result);
    163351  return result;
    164352}
    … …  
    169357  ) const
    170358{
     359//  Info info(__func__);
    171360  result_t result = 0.;
    172361
    173362  // get all triples within the cutoff
    174   std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
     363  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
    175364  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
    176365      iter != triples.end(); ++iter) {
    … …  
    181370    result +=
    182371        function_cutoff(r_ik.distance)
     372        * params[omega]
    183373        * function_angle(r_ij.distance, r_ik.distance, r_jk.distance)
    184         * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
     374        * exp( Helpers::pow(params[lambda3] * (r_ij.distance - r_ik.distance) ,3));
    185375  }
    186376
     377//  LOG(2, "DEBUG: function_zeta(" << r_ij.distance << ") = " << result);
    187378  return result;
    188379}
    … …  
    195386  ) const
    196387{
     388//  Info info(__func__);
    197389  const double angle = Helpers::pow(r_ij,2) + Helpers::pow(r_ik,2) - Helpers::pow(r_jk,2);
    198   const double divisor = r_ij * r_ik;
    199   const double result =
    200       1.
    201       + (Helpers::pow(c, 2)/Helpers::pow(d, 2))
    202       - Helpers::pow(c, 2)/(Helpers::pow(d, 2) +
    203           Helpers::pow(h - cos(angle/divisor),2));
    204   return result;
    205 }
     390  const double divisor = 2.* r_ij * r_ik;
     391//  LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
     392  if (divisor != 0.) {
     393    const double result =
     394        1.
     395        + (Helpers::pow(params[c]/params[d], 2))
     396        - Helpers::pow(params[c], 2)/(Helpers::pow(params[d], 2) +
     397            Helpers::pow(params[h] - angle/divisor,2));
     398
     399//    LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
     400    return result;
     401  } else
     402    return 0.;
     403}
     404

  • src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp

    re7579e rffc368  
    1717#include <cmath>
    1818
     19#include "CodePatterns/Assert.hpp"
     20
    1921#include "Potentials/EmpiricalPotential.hpp"
    2022#include "FunctionApproximation/FunctionModel.hpp"
    21 
    22 #include "CodePatterns/Assert.hpp"
    2323
    2424/** This class is the implementation of the Tersoff potential function.
    … …  
    3131class ManyBodyPotential_Tersoff : virtual public EmpiricalPotential, virtual public FunctionModel
    3232{
     33  //!> grant unit test access to internal parts
     34  friend class ManyBodyPotential_TersoffTest;
    3335  // some repeated typedefs to avoid ambiguities
    3436  typedef FunctionModel::arguments_t arguments_t;
    … …  
    5052  /** Constructor for class ManyBodyPotential_Tersoff.
    5153   *
    52    * @param _cutoff_offset offset for cutoff
    53    * @param _cutoff_halfwidth halfwidth for cutoff relative to \a _cutoff_offset
     54   * @param _R offset for cutoff
     55   * @param _S halfwidth for cutoff relative to \a _R
    5456   * @param A
    5557   * @param B
    56    * @param lambda1
    57    * @param lambda2
     58   * @param lambda
     59   * @param mu
    5860   * @param lambda3
    5961   * @param alpha
    6062   * @param beta
     63   * @param chi
     64   * @param omega
    6165   * @param n
    6266   * @param c
    … …  
    6872   */
    6973  ManyBodyPotential_Tersoff(
    70       const double &_cutoff_offset,
    71       const double &_cutoff_halfwidth,
    72       const double &A,
    73       const double &B,
    74       const double &lambda1,
    75       const double &lambda2,
    76       const double &lambda3,
    77       const double &alpha,
    78       const double &beta,
    79       const double &n,
    80       const double &c,
    81       const double &d,
    82       const double &h,
     74      const double &_R,
     75      const double &_S,
     76      const double &_A,
     77      const double &_B,
     78      const double &_lambda,
     79      const double &_mu,
     80      const double &_lambda3,
     81      const double &_alpha,
     82      const double &_beta,
     83      const double &_chi,
     84      const double &_omega,
     85      const double &_n,
     86      const double &_c,
     87      const double &_d,
     88      const double &_h,
    8389      boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction);
    8490
    … …  
    132138  /** This  function has the exponential feature from the Morse potential.
    133139   *
    134    * @param distance variable of the function
    135140   * @param prefactor prefactor parameter to exp function
    136141   * @param lambda scale parameter of exp function's argument
     142   * @param distance variable of the function
    137143   * @return
    138144   */
    139145  result_t function_smoother(
    140       const double &distance,
    141146      const double &prefactor,
    142       const double &lambda
    143       ) const
    144   {
    145     return prefactor * exp(-lambda * distance);
    146   }
     147      const double &lambda,
     148      const double &distance
     149      ) const;
    147150
    148151  /** This function represents \f$ (1 + \alpha^n \eta^n)^{-1/2n} \f$.
    … …  
    150153   * @param alpha prefactor to eta function
    151154   * @param r_ij distance argument
    152    * @param etafunction eta or zeta
     155   * @param eta result value of eta or zeta
    153156   * @return \f$ (1 + \alpha^n \eta^n)^{-1/2n} \f$
    154157   */
    155158  result_t function_prefactor(
    156159      const double &alpha,
    157       boost::function<result_t()> etafunction
     160      const double &eta
    158161        ) const;
    159162
    … …  
    177180private:
    178181  enum parameter_enum_t {
    179     cutoff_offset=0,
    180     equilibrium_distance=1,
    181     cutoff_halfwidth=2,
    182     A=3,
    183     B=4,
    184     lambda1=5,
    185     lambda2=6,
    186     lambda3=7,
    187     alpha=8,
    188     beta=9,
    189     n=10,
    190     c=11,
    191     d=12,
    192     h=13,
     182    R=0,
     183    S=1,
     184    A=2,
     185    B=3,
     186    lambda=4,
     187    mu=5,
     188    lambda3=6,
     189    alpha=7,
     190    beta=8,
     191    chi=9,
     192    omega=10,
     193    n=11,
     194    c=12,
     195    d=13,
     196    h=14,
    193197    MAXPARAMS
    194198  };

  • src/Potentials/Specifics/unittests/Makefile.am

    re7579e rffc368  
    33
    44POTENTIALSSPECIFICSTESTSSOURCES = \
     5        ../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp \
    56        ../Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.cpp \
    67        ../Potentials/Specifics/unittests/PairPotential_MorseUnitTest.cpp
    78
    89POTENTIALSSPECIFICSTESTSHEADERS = \
     10        ../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.hpp \
    911        ../Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.hpp \
    1012        ../Potentials/Specifics/unittests/PairPotential_MorseUnitTest.hpp
    1113
    1214POTENTIALSSPECIFICSTESTS = \
     15        ManyBodyPotential_TersoffUnitTest \
    1316        PairPotential_HarmonicUnitTest \
    1417        PairPotential_MorseUnitTest
    … …  
    2023POTENTIALSSPECIFICSLIBS = \
    2124        ../libMolecuilderPotentials.la \
     25        $(top_builddir)/LinearAlgebra/src/LinearAlgebra/libLinearAlgebra.la \
    2226        ${CodePatterns_LIBS} \
    2327  $(BOOST_SERIALIZATION_LDFLAGS) $(BOOST_SERIALIZATION_LIBS) \
    2428        $(BOOST_LIB)
     29
     30ManyBodyPotential_TersoffUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
     31        ../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp \
     32        ../Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.hpp
     33ManyBodyPotential_TersoffUnitTest_LDADD = ${POTENTIALSSPECIFICSLIBS}
    2534
    2635PairPotential_HarmonicUnitTest_SOURCES = $(top_srcdir)/src/unittests/UnitTestMain.cpp \
Note: See TracChangeset for help on using the changeset viewer.