Ignore:
Timestamp:
Jun 27, 2014, 9:32:55 PM (11 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:
550f2a
Parents:
16227a
git-author:
Frederik Heber <heber@…> (02/27/14 20:15:41)
git-committer:
Frederik Heber <heber@…> (06/27/14 21:32:55)
Message:

FunctionModel now uses list_of_arguments to split sequence of subsets of distances.

  • this fixes ambiguities with the set of distances: Imagine the distances within a water molecule as OH (A) and HH (B). We then may have a sequence of argument_t as AABAAB. And with the current implementation of CompoundPotential::splitUpArgumentsByModels() we would always choose the latter (and more complex) model. Hence, we make two calls to TriplePotential_Angle, instead of calls twice to PairPotential_Harmonic for A, one to PairPotential_Harmonic for B, and once to TriplePotential_Angle for AAB.
  • now, we new list looks like A,A,B,AAB where each tuple of distances can be uniquely associated with a specific potential.
  • changed signatures of EmpiricalPotential::operator(), ::derivative(), ::parameter_derivative(). This involved changing all of the current specific potentials and CompoundPotential.
  • TrainingData must discern between the InputVector_t (just all distances) and the FilteredInputVector_t (tuples of subsets of distances).
  • FunctionApproximation now has list_of_arguments_t as parameter to evaluate() and evaluate_derivative().
  • DOCU: docu change in TrainingData.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp

    r16227a re1fe7e  
    184184ManyBodyPotential_Tersoff::results_t
    185185ManyBodyPotential_Tersoff::operator()(
    186     const arguments_t &arguments
     186    const list_of_arguments_t &listarguments
    187187    ) const
    188188{
    189189//  Info info(__func__);
    190   double result = 0.;
    191   for(arguments_t::const_iterator argiter = arguments.begin();
    192       argiter != arguments.end();
    193       ++argiter) {
    194     const argument_t &r_ij = *argiter;
    195     ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    196         arguments_t(1, r_ij), getParticleTypes()),
    197         "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
    198  
    199     const double cutoff = function_cutoff(r_ij.distance);
    200     const double temp = (cutoff == 0.) ?
    201         0. :
    202         cutoff * (
    203             function_prefactor(
    204                 alpha,
    205                 function_eta(r_ij))
     190  result_t result = 0.;
     191  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     192      iter != listarguments.end(); ++iter) {
     193    const arguments_t &arguments = *iter;
     194    for(arguments_t::const_iterator argiter = arguments.begin();
     195        argiter != arguments.end();
     196        ++argiter) {
     197      const argument_t &r_ij = *argiter;
     198      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     199          arguments_t(1, r_ij), getParticleTypes()),
     200          "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
     201
     202      const double cutoff = function_cutoff(r_ij.distance);
     203      const double temp = (cutoff == 0.) ?
     204          0. :
     205          cutoff * (
     206              function_prefactor(
     207                  alpha,
     208                  function_eta(r_ij))
     209              * function_smoother(
     210                  params[A],
     211                  params[lambda],
     212                  r_ij.distance)
     213              +
     214              function_prefactor(
     215                  params[beta],
     216                  function_zeta(r_ij))
     217              * function_smoother(
     218                  -params[B],
     219                  params[mu],
     220                  r_ij.distance)
     221          );
     222      result += temp;
     223    }
     224  }
     225//  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
     226  return results_t(1, result);
     227}
     228
     229ManyBodyPotential_Tersoff::derivative_components_t
     230ManyBodyPotential_Tersoff::derivative(
     231    const list_of_arguments_t &listarguments
     232    ) const
     233{
     234//  Info info(__func__);
     235  return derivative_components_t();
     236}
     237
     238ManyBodyPotential_Tersoff::results_t
     239ManyBodyPotential_Tersoff::parameter_derivative(
     240    const list_of_arguments_t &listarguments,
     241    const size_t index
     242    ) const
     243{
     244//  Info info(__func__);
     245//  ASSERT( arguments.size() == 1,
     246//      "ManyBodyPotential_Tersoff::parameter_derivative() - requires exactly one argument.");
     247
     248  result_t result = 0.;
     249  for(list_of_arguments_t::const_iterator iter = listarguments.begin();
     250      iter != listarguments.end(); ++iter) {
     251    const arguments_t &arguments = *iter;
     252    for(arguments_t::const_iterator argiter = arguments.begin();
     253        argiter != arguments.end();
     254        ++argiter) {
     255      const argument_t &r_ij = *argiter;
     256      ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
     257          arguments_t(1, r_ij), getParticleTypes()),
     258          "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
     259
     260    switch (index) {
     261//      case R:
     262//      {
     263//        result += 0.;
     264//        break;
     265//      }
     266//      case S:
     267//      {
     268//        result += 0.;
     269//        break;
     270//      }
     271      case A:
     272      {
     273        const double cutoff = function_cutoff(r_ij.distance);
     274        result += (cutoff == 0.) ?
     275            0. :
     276            cutoff *
     277                function_prefactor(
     278                    alpha,
     279                    function_eta(r_ij))
     280                * function_smoother(
     281                    1.,
     282                    params[lambda],
     283                    r_ij.distance);
     284//            cutoff * function_prefactor(
     285//                alpha,
     286//                function_eta(r_ij))
     287//            * function_smoother(
     288//                1.,
     289//                params[mu],
     290//                r_ij.distance);
     291        break;
     292      }
     293      case B:
     294      {
     295        const double cutoff = function_cutoff(r_ij.distance);
     296        result += (cutoff == 0.) ?
     297            0. :
     298            cutoff * function_prefactor(
     299                params[beta],
     300                function_zeta(r_ij))
    206301            * function_smoother(
    207                 params[A],
    208                 params[lambda],
    209                 r_ij.distance)
    210             +
     302                -1.,
     303                params[mu],
     304                r_ij.distance);
     305//            cutoff * function_prefactor(
     306//                beta,
     307//                function_zeta(r_ij))
     308//            * function_smoother(
     309//                -params[B],
     310//                params[mu],
     311//                r_ij.distance)/params[B];
     312        break;
     313      }
     314      case lambda:
     315      {
     316        const double cutoff = function_cutoff(r_ij.distance);
     317        result += (cutoff == 0.) ?
     318            0. :
     319            -r_ij.distance * cutoff *
     320                function_prefactor(
     321                    alpha,
     322                    function_eta(r_ij))
     323                * function_smoother(
     324                    params[A],
     325                    params[lambda],
     326                    r_ij.distance);
     327        break;
     328      }
     329      case mu:
     330      {
     331        const double cutoff = function_cutoff(r_ij.distance);
     332        result += (cutoff == 0.) ?
     333            0. :
     334            -r_ij.distance * cutoff *(
    211335            function_prefactor(
    212336                params[beta],
    … …  
    217341                r_ij.distance)
    218342        );
    219     result += temp;
    220   }
    221 //  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
    222   return std::vector<result_t>(1, result);
    223 }
    224 
    225 ManyBodyPotential_Tersoff::derivative_components_t
    226 ManyBodyPotential_Tersoff::derivative(
    227     const arguments_t &arguments
    228     ) const
    229 {
    230 //  Info info(__func__);
    231   return ManyBodyPotential_Tersoff::derivative_components_t();
    232 }
    233 
    234 ManyBodyPotential_Tersoff::results_t
    235 ManyBodyPotential_Tersoff::parameter_derivative(
    236     const arguments_t &arguments,
    237     const size_t index
    238     ) const
    239 {
    240 //  Info info(__func__);
    241 //  ASSERT( arguments.size() == 1,
    242 //      "ManyBodyPotential_Tersoff::parameter_derivative() - requires exactly one argument.");
    243 
    244   double result = 0.;
    245   for(arguments_t::const_iterator argiter = arguments.begin();
    246       argiter != arguments.end();
    247       ++argiter) {
    248     const argument_t &r_ij = *argiter;
    249     ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
    250         arguments_t(1, r_ij), getParticleTypes()),
    251         "ManyBodyPotential_Tersoff::operator() - types don't match with ones in arguments.");
    252 
    253   switch (index) {
    254 //    case R:
    255 //    {
    256 //      result += 0.;
    257 //      break;
    258 //    }
    259 //    case S:
    260 //    {
    261 //      result += 0.;
    262 //      break;
    263 //    }
    264     case A:
    265     {
    266       const double cutoff = function_cutoff(r_ij.distance);
    267       result += (cutoff == 0.) ?
    268           0. :
    269           cutoff *
    270               function_prefactor(
    271                   alpha,
    272                   function_eta(r_ij))
    273               * function_smoother(
    274                   1.,
    275                   params[lambda],
    276                   r_ij.distance);
    277 //          cutoff * function_prefactor(
    278 //              alpha,
    279 //              function_eta(r_ij))
    280 //          * function_smoother(
    281 //              1.,
    282 //              params[mu],
    283 //              r_ij.distance);
    284       break;
     343        break;
     344      }
     345//      case lambda3:
     346//      {
     347//        result += 0.;
     348//        break;
     349//      }
     350//      case alpha:
     351//      {
     352//        const double temp =
     353//            pow(alpha*function_eta(r_ij), params[n]);
     354//        const double cutoff = function_cutoff(r_ij.distance);
     355//        result += (cutoff == 0.) || (alpha == 0. )?
     356//            0. :
     357//            function_smoother(
     358//                params[A],
     359//                params[lambda],
     360//                r_ij.distance)
     361//            * (-.5) * alpha * (temp/alpha)
     362//            / (1. + temp)
     363//            ;
     364//        break;
     365//      }
     366//      case chi:
     367//      {
     368//        result += 0.;
     369//        break;
     370//      }
     371//      case omega:
     372//      {
     373//        result += 0.;
     374//        break;
     375//      }
     376      case beta:
     377      {
     378        const double temp =
     379            pow(params[beta]*function_zeta(r_ij), params[n]);
     380        const double cutoff = function_cutoff(r_ij.distance);
     381        result += (cutoff == 0.) || (params[beta] == 0. )?
     382            0. : cutoff *
     383            function_smoother(
     384                -params[B],
     385                params[mu],
     386                r_ij.distance)
     387            * (-.5)
     388            * function_prefactor(
     389                params[beta],
     390                function_zeta(r_ij))
     391            * (temp/params[beta])
     392            / (1. + temp)
     393            ;
     394        break;
     395      }
     396      case n:
     397      {
     398        const double zeta = function_zeta(r_ij);
     399        const double temp = pow( params[beta]*zeta , params[n]);
     400        const double cutoff = function_cutoff(r_ij.distance);
     401        const double tempres = ((cutoff == 0.) || (zeta == 0.)) ? // zeta must be caught if zero due to log
     402            0. : .5 * cutoff *
     403            function_smoother(
     404                -params[B],
     405                params[mu],
     406                r_ij.distance)
     407            * function_prefactor(
     408                params[beta],
     409                function_zeta(r_ij))
     410            * ( log(1.+temp)/(params[n]*params[n]) - temp
     411                * (log(function_zeta(r_ij)) + log(params[beta]))
     412                /(params[n]*(1.+temp)))
     413            ;
     414//        if (tempres != tempres)
     415//    LOG(2, "DEBUG: tempres is NaN.");
     416//        LOG(2, "DEBUG: Adding " << tempres << " for p.d. w.r.t n, temp=" << temp << ", cutoff=" << cutoff);
     417        result += tempres;
     418        break;
     419      }
     420      case c:
     421      {
     422        const double zeta = function_zeta(r_ij);
     423        if (zeta == 0.)
     424          break;
     425        const double temp =
     426            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     427        const double cutoff = function_cutoff(r_ij.distance);
     428        const double tempres = (cutoff == 0.) ?
     429            0. : cutoff *
     430            function_smoother(
     431                -params[B],
     432                params[mu],
     433                r_ij.distance)
     434            * function_prefactor(
     435                params[beta],
     436                zeta)
     437             * (-1.) * temp / (1.+temp*zeta);
     438        double factor = function_derivative_c(r_ij);
     439        result += tempres*factor;
     440        if (result != result)
     441          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
     442        break;
     443      }
     444      case d:
     445      {
     446        const double zeta = function_zeta(r_ij);
     447        const double temp =
     448            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     449        const double cutoff = function_cutoff(r_ij.distance);
     450        const double tempres = (cutoff == 0.) ?
     451            0. : cutoff *
     452            function_smoother(
     453                -params[B],
     454                params[mu],
     455                r_ij.distance)
     456            * function_prefactor(
     457                params[beta],
     458                zeta)
     459             * (-1.) * temp / (1.+temp*zeta);
     460        double factor = function_derivative_d(r_ij);
     461        result += tempres*factor;
     462        if (result != result)
     463          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
     464        break;
     465      }
     466      case h:
     467      {
     468        const double zeta = function_zeta(r_ij);
     469        const double temp =
     470            pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
     471        const double cutoff = function_cutoff(r_ij.distance);
     472        const double tempres = (cutoff == 0.) ?
     473            0. : cutoff *
     474            function_smoother(
     475                -params[B],
     476                params[mu],
     477                r_ij.distance)
     478            * function_prefactor(
     479                params[beta],
     480                zeta)
     481             * (-1.) * temp / (1.+temp*zeta);
     482        double factor = function_derivative_h(r_ij);
     483        result += tempres*factor;
     484        if (result != result)
     485          ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
     486        break;
     487      }
     488      default:
     489        ASSERT(0, "ManyBodyPotential_Tersoff::parameter_derivative() - derivative to unknown parameter desired.");
     490        break;
    285491    }
    286     case B:
    287     {
    288       const double cutoff = function_cutoff(r_ij.distance);
    289       result += (cutoff == 0.) ?
    290           0. :
    291           cutoff * function_prefactor(
    292               params[beta],
    293               function_zeta(r_ij))
    294           * function_smoother(
    295               -1.,
    296               params[mu],
    297               r_ij.distance);
    298 //          cutoff * function_prefactor(
    299 //              beta,
    300 //              function_zeta(r_ij))
    301 //          * function_smoother(
    302 //              -params[B],
    303 //              params[mu],
    304 //              r_ij.distance)/params[B];
    305       break;
     492    if (result != result)
     493      ELOG(1, "result is NaN.");
    306494    }
    307     case lambda:
    308     {
    309       const double cutoff = function_cutoff(r_ij.distance);
    310       result += (cutoff == 0.) ?
    311           0. :
    312           -r_ij.distance * cutoff *
    313               function_prefactor(
    314                   alpha,
    315                   function_eta(r_ij))
    316               * function_smoother(
    317                   params[A],
    318                   params[lambda],
    319                   r_ij.distance);
    320       break;
    321     }
    322     case mu:
    323     {
    324       const double cutoff = function_cutoff(r_ij.distance);
    325       result += (cutoff == 0.) ?
    326           0. :
    327           -r_ij.distance * cutoff *(
    328           function_prefactor(
    329               params[beta],
    330               function_zeta(r_ij))
    331           * function_smoother(
    332               -params[B],
    333               params[mu],
    334               r_ij.distance)
    335       );
    336       break;
    337     }
    338 //    case lambda3:
    339 //    {
    340 //      result += 0.;
    341 //      break;
    342 //    }
    343 //    case alpha:
    344 //    {
    345 //      const double temp =
    346 //          pow(alpha*function_eta(r_ij), params[n]);
    347 //      const double cutoff = function_cutoff(r_ij.distance);
    348 //      result += (cutoff == 0.) || (alpha == 0. )?
    349 //          0. :
    350 //          function_smoother(
    351 //              params[A],
    352 //              params[lambda],
    353 //              r_ij.distance)
    354 //          * (-.5) * alpha * (temp/alpha)
    355 //          / (1. + temp)
    356 //          ;
    357 //      break;
    358 //    }
    359 //    case chi:
    360 //    {
    361 //      result += 0.;
    362 //      break;
    363 //    }
    364 //    case omega:
    365 //    {
    366 //      result += 0.;
    367 //      break;
    368 //    }
    369     case beta:
    370     {
    371       const double temp =
    372           pow(params[beta]*function_zeta(r_ij), params[n]);
    373       const double cutoff = function_cutoff(r_ij.distance);
    374       result += (cutoff == 0.) || (params[beta] == 0. )?
    375           0. : cutoff *
    376           function_smoother(
    377               -params[B],
    378               params[mu],
    379               r_ij.distance)
    380           * (-.5)
    381           * function_prefactor(
    382               params[beta],
    383               function_zeta(r_ij))
    384           * (temp/params[beta])
    385           / (1. + temp)
    386           ;
    387       break;
    388     }
    389     case n:
    390     {
    391       const double zeta = function_zeta(r_ij);
    392       const double temp = pow( params[beta]*zeta , params[n]);
    393       const double cutoff = function_cutoff(r_ij.distance);
    394       const double tempres = ((cutoff == 0.) || (zeta == 0.)) ? // zeta must be caught if zero due to log
    395           0. : .5 * cutoff *
    396           function_smoother(
    397               -params[B],
    398               params[mu],
    399               r_ij.distance)
    400           * function_prefactor(
    401               params[beta],
    402               function_zeta(r_ij))
    403           * ( log(1.+temp)/(params[n]*params[n]) - temp
    404               * (log(function_zeta(r_ij)) + log(params[beta]))
    405               /(params[n]*(1.+temp)))
    406           ;
    407 //      if (tempres != tempres)
    408 //      LOG(2, "DEBUG: tempres is NaN.");
    409 //      LOG(2, "DEBUG: Adding " << tempres << " for p.d. w.r.t n, temp=" << temp << ", cutoff=" << cutoff);
    410       result += tempres;
    411       break;
    412     }
    413     case c:
    414     {
    415       const double zeta = function_zeta(r_ij);
    416       if (zeta == 0.)
    417         break;
    418       const double temp =
    419           pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
    420       const double cutoff = function_cutoff(r_ij.distance);
    421       const double tempres = (cutoff == 0.) ?
    422           0. : cutoff *
    423           function_smoother(
    424               -params[B],
    425               params[mu],
    426               r_ij.distance)
    427           * function_prefactor(
    428               params[beta],
    429               zeta)
    430            * (-1.) * temp / (1.+temp*zeta);
    431       double factor = function_derivative_c(r_ij);
    432       result += tempres*factor;
    433       if (result != result)
    434         ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
    435       break;
    436     }
    437     case d:
    438     {
    439       const double zeta = function_zeta(r_ij);
    440       const double temp =
    441           pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
    442       const double cutoff = function_cutoff(r_ij.distance);
    443       const double tempres = (cutoff == 0.) ?
    444           0. : cutoff *
    445           function_smoother(
    446               -params[B],
    447               params[mu],
    448               r_ij.distance)
    449           * function_prefactor(
    450               params[beta],
    451               zeta)
    452            * (-1.) * temp / (1.+temp*zeta);
    453       double factor = function_derivative_d(r_ij);
    454       result += tempres*factor;
    455       if (result != result)
    456         ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
    457       break;
    458     }
    459     case h:
    460     {
    461       const double zeta = function_zeta(r_ij);
    462       const double temp =
    463           pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
    464       const double cutoff = function_cutoff(r_ij.distance);
    465       const double tempres = (cutoff == 0.) ?
    466           0. : cutoff *
    467           function_smoother(
    468               -params[B],
    469               params[mu],
    470               r_ij.distance)
    471           * function_prefactor(
    472               params[beta],
    473               zeta)
    474            * (-1.) * temp / (1.+temp*zeta);
    475       double factor = function_derivative_h(r_ij);
    476       result += tempres*factor;
    477       if (result != result)
    478         ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
    479       break;
    480     }
    481     default:
    482       ASSERT(0, "ManyBodyPotential_Tersoff::parameter_derivative() - derivative to unknown parameter desired.");
    483       break;
    484   }
    485   if (result != result)
    486     ELOG(1, "result is NaN.");
    487495  }
    488496  return results_t(1,-result);
Note: See TracChangeset for help on using the changeset viewer.