Changeset ca8d82 for src


Ignore:
Timestamp:
Dec 19, 2012, 3:25:54 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:
2e9486
Parents:
f904d5
git-author:
Frederik Heber <heber@…> (10/07/12 15:28:49)
git-committer:
Frederik Heber <heber@…> (12/19/12 15:25:54)
Message:

FIX: Corrected parameter derivatives of ManyBodyPotential_Tersoff.

  • setParameters also allows setting of less parameters than required.
Location:
src/Potentials/Specifics
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

    rf904d5 rca8d82  
    170170    case A:
    171171    {
    172       const double result = 0.;
     172      const double cutoff = function_cutoff(r_ij.distance);
     173      const double result = (cutoff == 0.) ?
     174          0. :
     175          cutoff *
     176              function_prefactor(
     177                  alpha,
     178                  function_eta(r_ij))
     179              * function_smoother(
     180                  1.,
     181                  params[lambda],
     182                  r_ij.distance);
     183//          cutoff * function_prefactor(
     184//              alpha,
     185//              function_eta(r_ij))
     186//          * function_smoother(
     187//              1.,
     188//              params[mu],
     189//              r_ij.distance);
    173190      return results_t(1, result);
    174191      break;
    … …  
    176193    case B:
    177194    {
    178       const double result = 0.;
     195      const double cutoff = function_cutoff(r_ij.distance);
     196      const double result = (cutoff == 0.) ?
     197          0. :
     198          cutoff * function_prefactor(
     199              params[beta],
     200              function_zeta(r_ij))
     201          * function_smoother(
     202              -1.,
     203              params[mu],
     204              r_ij.distance);
     205//          cutoff * function_prefactor(
     206//              beta,
     207//              function_zeta(r_ij))
     208//          * function_smoother(
     209//              -params[B],
     210//              params[mu],
     211//              r_ij.distance)/params[B];
    179212      return results_t(1, result);
    180213      break;
    … …  
    185218      const double result = (cutoff == 0.) ?
    186219          0. :
    187           -r_ij.distance * cutoff * params[lambda] * (
    188           function_prefactor(
    189               alpha,
    190               function_eta(r_ij))
    191           * function_smoother(
    192               params[A],
    193               params[lambda],
    194               r_ij.distance)
    195           );
     220          -r_ij.distance * cutoff *
     221              function_prefactor(
     222                  alpha,
     223                  function_eta(r_ij))
     224              * function_smoother(
     225                  params[A],
     226                  params[lambda],
     227                  r_ij.distance);
    196228      return results_t(1, result);
    197229      break;
    … …  
    228260//          0. :
    229261//          function_smoother(
    230 //              -params[A],
     262//              params[A],
    231263//              params[lambda],
    232264//              r_ij.distance)
    … …  
    294326    {
    295327      const double zeta = function_zeta(r_ij);
     328      if (zeta == 0.)
     329        break;
    296330      const double temp =
    297           pow(params[beta]*zeta, params[n]);
    298       const double cutoff = function_cutoff(r_ij.distance);
    299       const double result = (cutoff == 0.) ?
    300           0. : - cutoff *
     331          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     332      const double cutoff = function_cutoff(r_ij.distance);
     333      const double result = (cutoff == 0.) ?
     334          0. : cutoff *
    301335          function_smoother(
    302336              -params[B],
    … …  
    306340              params[beta],
    307341              zeta)
    308           * (temp/zeta)/(1+(temp/zeta));
    309       double factor = 0.;
    310       std::vector<arguments_t> triples = triplefunction(r_ij, S);
    311       for (std::vector<arguments_t>::const_iterator iter = triples.begin();
    312           iter != triples.end(); ++iter) {
    313         ASSERT( iter->size() == 2,
    314             "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain exactly two distances.");
    315         const argument_t &r_ik = (*iter)[0];
    316         const argument_t &r_jk = (*iter)[1];
    317         factor +=
    318             function_cutoff(r_ik.distance)
    319             * (params[c]/Helpers::pow(params[c],2))
    320             * (1. - 1./(1.+(params[h]-function_theta(r_ij.distance, r_ik.distance, r_jk.distance))/Helpers::pow(params[c],2)));
    321       }
     342           * (-1.) * temp / (1.+temp*zeta);
     343      double factor = function_derivative_c(r_ij);
    322344      return results_t(1, result*factor);
    323345      break;
    … …  
    327349      const double zeta = function_zeta(r_ij);
    328350      const double temp =
    329           pow(params[beta]*zeta, params[n]);
    330       const double cutoff = function_cutoff(r_ij.distance);
    331       const double result = (cutoff == 0.) ?
    332           0. : - cutoff *
     351          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     352      const double cutoff = function_cutoff(r_ij.distance);
     353      const double result = (cutoff == 0.) ?
     354          0. : cutoff *
    333355          function_smoother(
    334356              -params[B],
    … …  
    338360              params[beta],
    339361              zeta)
    340           * (temp/zeta)/(1+(temp/zeta));
    341       double factor = 0.;
    342       std::vector<arguments_t> triples = triplefunction(r_ij, S);
    343       for (std::vector<arguments_t>::const_iterator iter = triples.begin();
    344           iter != triples.end(); ++iter) {
    345         ASSERT( iter->size() == 2,
    346             "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain exactly two distances.");
    347         const argument_t &r_ik = (*iter)[0];
    348         const argument_t &r_jk = (*iter)[1];
    349         factor +=
    350             function_cutoff(r_ik.distance)
    351             * (-Helpers::pow(params[c],2)/Helpers::pow(params[d],3)
    352                 - Helpers::pow(params[c],2) * params[d]
    353                  / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance),2),2)
    354             );
    355       }
     362           * (-1.) * temp / (1.+temp*zeta);
     363      double factor = function_derivative_d(r_ij);
    356364      return results_t(1, result*factor);
    357365      break;
    … …  
    361369      const double zeta = function_zeta(r_ij);
    362370      const double temp =
    363           pow(params[beta]*zeta, params[n]);
    364       const double cutoff = function_cutoff(r_ij.distance);
    365       const double result = (cutoff == 0.) ?
    366           0. : - cutoff *
     371          pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     372      const double cutoff = function_cutoff(r_ij.distance);
     373      const double result = (cutoff == 0.) ?
     374          0. : cutoff *
    367375          function_smoother(
    368376              -params[B],
    … …  
    372380              params[beta],
    373381              zeta)
    374           * (temp/zeta)/(1+(temp/zeta));
    375       double factor = 0.;
    376       std::vector<arguments_t> triples = triplefunction(r_ij, S);
    377       for (std::vector<arguments_t>::const_iterator iter = triples.begin();
    378           iter != triples.end(); ++iter) {
    379         ASSERT( iter->size() == 2,
    380             "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain exactly two distances.");
    381         const argument_t &r_ik = (*iter)[0];
    382         const argument_t &r_jk = (*iter)[1];
    383         const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
    384         factor +=
    385             function_cutoff(r_ik.distance)
    386             * (-1.*Helpers::pow(params[c],2)*tempangle)
    387             / Helpers::pow(Helpers::pow(params[d],2)+Helpers::pow(tempangle,2),2)
    388             ;
    389       }
     382           * (-1.) * temp / (1.+temp*zeta);
     383      double factor = function_derivative_h(r_ij);
    390384      return results_t(1, result*factor);
    391385      break;
    … …  
    395389  }
    396390  return results_t(1, 0.);
     391}
     392
     393ManyBodyPotential_Tersoff::result_t
     394ManyBodyPotential_Tersoff::function_derivative_c(
     395    const argument_t &r_ij
     396  ) const
     397{
     398  double result = 0.;
     399  std::vector<arguments_t> triples = triplefunction(r_ij, S);
     400  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
     401      iter != triples.end(); ++iter) {
     402    ASSERT( iter->size() == 2,
     403        "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
     404    const argument_t &r_ik = (*iter)[0];
     405    const argument_t &r_jk = (*iter)[1];
     406    const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
     407    const double cutoff = function_cutoff(r_ik.distance);
     408    result += (cutoff == 0.) ?
     409        0. : cutoff * omega * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
     410            params[c]/Helpers::pow(params[d],2)
     411            - params[c] / ( Helpers::pow(params[d],2) + Helpers::pow(tempangle,2) )
     412        );
     413  }
     414  return result;
     415}
     416
     417ManyBodyPotential_Tersoff::result_t
     418ManyBodyPotential_Tersoff::function_derivative_d(
     419    const argument_t &r_ij
     420  ) const
     421{
     422  double result = 0.;
     423  std::vector<arguments_t> triples = triplefunction(r_ij, S);
     424  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
     425      iter != triples.end(); ++iter) {
     426    ASSERT( iter->size() == 2,
     427        "ManyBodyPotential_Tersoff::function_derivative_d() - the triples result must contain exactly two distances.");
     428    const argument_t &r_ik = (*iter)[0];
     429    const argument_t &r_jk = (*iter)[1];
     430    const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
     431    const double cutoff = function_cutoff(r_ik.distance);
     432    result += (cutoff == 0.) ?
     433        0. : cutoff * omega  * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
     434          - Helpers::pow(params[c],2)/Helpers::pow(params[d],3)
     435          + Helpers::pow(params[c],2) * params[d]
     436            / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(tempangle,2),2)
     437        );
     438  }
     439  return result;
     440}
     441
     442ManyBodyPotential_Tersoff::result_t
     443ManyBodyPotential_Tersoff::function_derivative_h(
     444    const argument_t &r_ij
     445  ) const
     446{
     447  double result = 0.;
     448  std::vector<arguments_t> triples = triplefunction(r_ij, S);
     449  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
     450      iter != triples.end(); ++iter) {
     451    ASSERT( iter->size() == 2,
     452        "ManyBodyPotential_Tersoff::function_derivative_h() - the triples result must contain exactly two distances.");
     453    const argument_t &r_ik = (*iter)[0];
     454    const argument_t &r_jk = (*iter)[1];
     455    const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
     456    const double cutoff = function_cutoff(r_ik.distance);
     457    result += (cutoff == 0.) ?
     458        0. : cutoff * omega  * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
     459          ( Helpers::pow(params[c],2)*tempangle )
     460            / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(tempangle,2),2)
     461        );
     462  }
     463  return result;
    397464}
    398465
    … …  
    516583{
    517584//  Info info(__func__);
    518   const double angle = function_theta(r_ij, r_ik, r_jk);
    519 //  LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
    520585  const double result =
    521586      1.
    522587      + (Helpers::pow(params[c]/params[d], 2))
    523588      - Helpers::pow(params[c], 2)/(Helpers::pow(params[d], 2) +
    524           Helpers::pow(params[h] - angle,2));
    525 
    526 //    LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
    527   return result;
    528 }
    529 
     589          Helpers::pow(params[h] - function_theta(r_ij, r_ik, r_jk),2));
     590
     591//  LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
     592  return result;
     593}
     594

  • src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp

    rf904d5 rca8d82  
    186186
    187187private:
     188  result_t
     189  function_derivative_c(
     190      const argument_t &r_ij
     191    ) const;
     192
     193  result_t
     194  function_derivative_d(
     195      const argument_t &r_ij
     196    ) const;
     197
     198  result_t
     199  function_derivative_h(
     200      const argument_t &r_ij
     201    ) const;
     202
     203public:
    188204  enum parameter_enum_t {
    189205    A,
    … …  
    223239  void setParameters(const parameters_t &_params)
    224240  {
    225     ASSERT( _params.size() == getParameterDimension(),
    226         "ManyBodyPotential_Tersoff::setParameters() - we need exactly "
     241    ASSERT( _params.size() <= getParameterDimension(),
     242        "ManyBodyPotential_Tersoff::setParameters() - we need not more than "
    227243        +toString(getParameterDimension())+" parameters.");
    228     params = _params;
     244    for (size_t i=0; i< _params.size(); ++i)
     245      params[i] = _params[i];
    229246  }
    230247
Note: See TracChangeset for help on using the changeset viewer.