Changeset e1fe7e for src/Potentials
- Timestamp:
- Jun 27, 2014, 9:32:55 PM (11 years ago)
- 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)
- Location:
- src/Potentials
- Files:
-
- 27 edited
Legend:
- Unmodified
- Added
- Removed
-
src/Potentials/CompoundPotential.cpp
r16227a re1fe7e 228 228 { 229 229 arguments_by_model_t partial_args; 230 // go through each model and have it filter out its arguments 230 // go through each model and have it filter out its arguments, this already 231 // returns a list of tuples associated with the specific model 231 232 for(models_t::const_iterator modeliter = models.begin(); 232 233 modeliter != models.end(); ++modeliter) { 233 234 FunctionModel::filter_t filterfunction = (*modeliter)->getSpecificFilter(); 234 arguments_t tempargs = filterfunction(arguments);235 list_of_arguments_t tempargs = filterfunction(arguments); 235 236 // then split up all the bunches, too. 236 arguments_t::const_iterator advanceiter = tempargs.begin(); 237 for (arguments_t::const_iterator argiter = tempargs.begin(); 238 argiter != tempargs.end(); argiter = advanceiter) { 239 advanceiter += (*modeliter)->getSpecificArgumentCount(); 237 for (list_of_arguments_t::const_iterator argiter = tempargs.begin(); 238 argiter != tempargs.end(); ++argiter) { 239 const arguments_t &args = *argiter; 240 240 partial_args.push_back( 241 241 std::make_pair( 242 242 *modeliter, 243 arg uments_t(argiter, advanceiter)243 args 244 244 ) 245 245 ); … … 251 251 252 252 CompoundPotential::arguments_by_model_t CompoundPotential::splitUpArgumentsByModels( 253 const arguments_t &arguments) const253 const list_of_arguments_t &listarguments) const 254 254 { 255 255 arguments_by_model_t partial_args; 256 arguments_t::const_iterator argiter = arguments.begin();257 256 particletypes_per_model_t::const_iterator typesiter = particletypes_per_model.begin(); 258 257 models_t::const_iterator modeliter = models.begin(); 259 258 260 // add constant model (which is always first model) with empty args if present259 /// add constant model (which is always first model) with empty args if present 261 260 if (typesiter->empty()) { 262 261 partial_args.push_back( … … 266 265 ++typesiter; 267 266 } 267 268 268 // then check other models 269 while (argiter != arguments.end()) { 269 /// we only have to check whether the current model still matches or whether 270 /// have to use the next model. 271 for (list_of_arguments_t::const_iterator argiter = listarguments.begin(); 272 argiter != listarguments.end(); ++argiter) { 273 const arguments_t &arguments = *argiter; 270 274 if (typesiter+1 != particletypes_per_model.end()) { 271 275 // check whether next argument bunch is for same model or different one … … 275 279 276 280 // we always expect N(N-1)/2 distances for N particle types 277 arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2); 278 arguments_t::const_iterator nextenditer = argiter+(nexttypes.size()*(nexttypes.size()-1)/2); 279 arguments_t args(argiter, enditer); 280 arguments_t nextargs(argiter, nextenditer); 281 if (types.size() < nexttypes.size()) { 282 if (areValidArguments(nexttypes, nextargs)) { 281 // check first from sizes alone 282 const size_t tuplesize = types.size()*(types.size()-1)/2; 283 const size_t nexttuplesize = nexttypes.size()*(nexttypes.size()-1)/2; 284 if ((tuplesize != nexttuplesize)) { 285 if ((arguments.size() == tuplesize) && areValidArguments(types, arguments)) { 286 // only former still matches, don't increment 287 partial_args.push_back( 288 std::make_pair(*modeliter, arguments) 289 ); 290 } else if ((arguments.size() == nexttuplesize) && areValidArguments(nexttypes, arguments)) { 283 291 // latter matches, increment 284 292 ++typesiter; 285 293 partial_args.push_back( 286 std::make_pair(*(++modeliter), arguments _t(argiter, nextenditer))294 std::make_pair(*(++modeliter), arguments) 287 295 ); 288 argiter = nextenditer; 289 } else if (areValidArguments(types, args)) { 290 // only former matches, don't increment 291 partial_args.push_back( 292 std::make_pair(*modeliter, arguments_t(argiter, enditer)) 293 ); 294 argiter = enditer; 295 } else 296 } else { 296 297 ASSERT(0, 297 "CompoundPotential::splitUpArgumentsByModels() - neither type matches its argument bunch."); 298 } else { 299 if (areValidArguments(types, args)) { 300 // only former matches, don't increment 301 partial_args.push_back( 302 std::make_pair(*modeliter, arguments_t(argiter, enditer)) 303 ); 304 argiter = enditer; 305 } else if (areValidArguments(nexttypes, nextargs)) { 306 // latter matches, increment 307 ++typesiter; 308 partial_args.push_back( 309 std::make_pair(*(++modeliter), arguments_t(argiter, nextenditer)) 310 ); 311 argiter = nextenditer; 298 "CompoundPotential::splitUpArgumentsByModels() - neither this model nor next model match (size) with current tuple."); 299 } 300 } else { // same size, now we have to check the types individually 301 size_t encodeValidity = 0; 302 encodeValidity += 1*areValidArguments(types, arguments); 303 encodeValidity += 2*areValidArguments(nexttypes, arguments); 304 305 switch (encodeValidity) { 306 case 1: 307 // only former still matches, don't increment 308 partial_args.push_back( 309 std::make_pair(*modeliter, arguments) 310 ); 311 break; 312 case 2: 313 ++typesiter; 314 partial_args.push_back( 315 std::make_pair(*(++modeliter), arguments) 316 ); 317 break; 318 case 0: 319 case 3: 320 default: 321 ASSERT(0, 322 "CompoundPotential::splitUpArgumentsByModels() - neither this model nor next model match (type) with current tuple."); 323 break; 312 324 } 313 325 } 314 326 } else { 315 327 const SerializablePotential::ParticleTypes_t &types = *typesiter; 316 // we always expect N(N-1)/2 distances for N particle types 317 arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2); 318 partial_args.push_back( 319 std::make_pair(*modeliter, arguments_t(argiter, enditer)) 320 ); 321 argiter = enditer; 328 if (areValidArguments(types, arguments)) { 329 // only former matches, don't increment 330 partial_args.push_back( 331 std::make_pair(*modeliter, arguments) 332 ); 333 } else { 334 ASSERT(0, 335 "CompoundPotential::splitUpArgumentsByModels() - last model does not match with current tuple."); 336 } 322 337 } 323 338 } … … 326 341 } 327 342 328 CompoundPotential::results_t CompoundPotential::operator()(const arguments_t &arguments) const 343 CompoundPotential::results_t CompoundPotential::operator()( 344 const list_of_arguments_t &listarguments) const 329 345 { 330 346 /// first, we have to split up the given arguments 331 347 arguments_by_model_t partial_args = 332 splitUpArgumentsByModels( arguments);348 splitUpArgumentsByModels(listarguments); 333 349 // print split up argument list for debugging 334 350 if (DoLog(4)) { … … 348 364 iter != partial_args.end(); ++iter) { 349 365 partial_results.push_back( 350 (*iter->first)(iter->second) 366 (*iter->first)( 367 FunctionModel::list_of_arguments_t(1, iter->second)) 351 368 ); 352 369 } … … 366 383 } 367 384 368 CompoundPotential::results_t CompoundPotential::parameter_derivative(const arguments_t &arguments, const size_t index) const 385 CompoundPotential::results_t CompoundPotential::parameter_derivative( 386 const list_of_arguments_t &listarguments, 387 const size_t index) const 369 388 { 370 389 // first, we have to split up the given arguments 371 390 arguments_by_model_t partial_args = 372 splitUpArgumentsByModels( arguments);391 splitUpArgumentsByModels(listarguments); 373 392 // then, with each bunch of arguments, we call the specific model 374 393 // get parameter dimensions per model … … 399 418 const size_t indexbase = (iter == dimensions.begin()) ? 0 : *(iter-1); 400 419 CompoundPotential::results_t results = 401 model->parameter_derivative(argiter->second, index-indexbase); 420 model->parameter_derivative( 421 FunctionModel::list_of_arguments_t(1, argiter->second), index-indexbase); 402 422 403 423 // either set results or add … … 462 482 // create initial returnfunction 463 483 FunctionModel::filter_t returnfunction = 464 boost::bind(&Helpers::returnEmpty Arguments);484 boost::bind(&Helpers::returnEmptyListArguments); 465 485 466 486 // every following fragments combines its arguments with the initial function … … 468 488 modeliter != models.end(); ++modeliter) { 469 489 returnfunction = 470 boost::bind(&Extractors::concatenate Arguments,490 boost::bind(&Extractors::concatenateListOfArguments, 471 491 boost::bind(returnfunction, _1), 472 492 boost::bind((*modeliter)->getSpecificFilter(), _1) -
src/Potentials/CompoundPotential.hpp
r16227a re1fe7e 93 93 /** Evaluates the harmonic potential function for the given arguments. 94 94 * 95 * @param arguments single distance95 * @param listarguments list of tuples of distances 96 96 * @return value of the potential function 97 97 */ 98 results_t operator()(const arguments_t &arguments) const;98 results_t operator()(const list_of_arguments_t &listarguments) const; 99 99 100 100 /** Evaluates the derivative of the function with the given \a arguments 101 101 * with respect to a specific parameter indicated by \a index. 102 102 * 103 * \param arguments setof arguments as input variables to the function103 * \param arguments list of sets of arguments as input variables to the function 104 104 * \param index derivative of which parameter 105 105 * \return result vector containing the derivative with respect to the given 106 106 * input 107 107 */ 108 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;108 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 109 109 110 110 /** States whether lower and upper boundaries should be used to constraint … … 153 153 /** Helper function to split up concatenated arguments for operator() calls. 154 154 * 155 * \param arguments arguments to split up155 * \param listarguments list of tuples of distances to assign to a model each 156 156 * \return vector of partial arguments with associated model 157 157 */ 158 arguments_by_model_t splitUpArgumentsByModels(const arguments_t &arguments) const;158 arguments_by_model_t splitUpArgumentsByModels(const list_of_arguments_t &listarguments) const; 159 159 160 160 /** Helper function to split up total list of arguments for operator() calls. -
src/Potentials/EmpiricalPotential.hpp
r16227a re1fe7e 60 60 * parameters. 61 61 * 62 * \param arguments setof arguments as input variables to the function62 * \param listarguments list of sets of arguments as input variables to the function 63 63 * \return result of the function 64 64 */ 65 virtual results_t operator()(const arguments_t &arguments) const=0;65 virtual results_t operator()(const list_of_arguments_t &listarguments) const=0; 66 66 67 67 /** Evaluates the derivative of the function with the given \a arguments 68 68 * for each component. 69 69 * 70 * \param arguments setof arguments as input variables to the function70 * \param listarguments list of sets of arguments as input variables to the function 71 71 * \return result vector containing the derivative with respect to each 72 72 * input comonent of the function 73 73 */ 74 virtual derivative_components_t derivative(const arguments_t &arguments) const=0;74 virtual derivative_components_t derivative(const list_of_arguments_t &listarguments) const=0; 75 75 76 76 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/ConstantPotential.cpp
r16227a re1fe7e 112 112 ConstantPotential::results_t 113 113 ConstantPotential::operator()( 114 const arguments_t &arguments114 const list_of_arguments_t &listarguments 115 115 ) const 116 116 { 117 ASSERT( arguments.size() == 0, 118 "ConstantPotential::operator() - requires no argument."); 119 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 120 arguments, getParticleTypes()), 121 "ConstantPotential::operator() - types don't match with ones in arguments."); 122 const result_t result = params[energy_offset]; 123 return std::vector<result_t>(1, result); 117 // directly set result as energy constant as independent of arg list 118 result_t result = params[energy_offset]; 119 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 120 iter != listarguments.end(); ++iter) { 121 const arguments_t &arguments = *iter; 122 ASSERT( arguments.size() == 0, 123 "ConstantPotential::operator() - requires no argument."); 124 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 125 arguments, getParticleTypes()), 126 "ConstantPotential::operator() - types don't match with ones in arguments."); 127 result += 0.; 128 } 129 return results_t(1, result); 124 130 } 125 131 126 132 ConstantPotential::derivative_components_t 127 133 ConstantPotential::derivative( 128 const arguments_t &arguments134 const list_of_arguments_t &listarguments 129 135 ) const 130 136 { 131 ASSERT( arguments.size() == 0, 132 "ConstantPotential::operator() - requires no argument."); 133 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 134 arguments, getParticleTypes()), 135 "ConstantPotential::operator() - types don't match with ones in arguments."); 136 derivative_components_t result(1, 0.); 137 return result; 137 result_t result = 0.; 138 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 139 iter != listarguments.end(); ++iter) { 140 const arguments_t &arguments = *iter; 141 ASSERT( arguments.size() == 0, 142 "ConstantPotential::operator() - requires no argument."); 143 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 144 arguments, getParticleTypes()), 145 "ConstantPotential::operator() - types don't match with ones in arguments."); 146 result += 0.; 147 } 148 return derivative_components_t(1, result); 138 149 } 139 150 140 151 ConstantPotential::results_t 141 152 ConstantPotential::parameter_derivative( 142 const arguments_t &arguments,153 const list_of_arguments_t &listarguments, 143 154 const size_t index 144 155 ) const 145 156 { 146 ASSERT( arguments.size() == 0, 147 "ConstantPotential::parameter_derivative() - requires no argument."); 148 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 149 arguments, getParticleTypes()), 150 "ConstantPotential::operator() - types don't match with ones in arguments."); 151 switch (index) { 152 case energy_offset: 153 { 154 // Maple result: 1 155 const result_t result = +1.; 156 return std::vector<result_t>(1, result); 157 break; 158 } 159 default: 160 break; 157 // Maple result: 1 158 result_t result = 1.; // energy_offset 159 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 160 iter != listarguments.end(); ++iter) { 161 const arguments_t &arguments = *iter; 162 ASSERT( arguments.size() == 0, 163 "ConstantPotential::parameter_derivative() - requires no argument."); 164 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 165 arguments, getParticleTypes()), 166 "ConstantPotential::operator() - types don't match with ones in arguments."); 167 // switch (index) { 168 // case energy_offset: 169 // { 170 // // Maple result: 1 171 // result = +1.; 172 // break; 173 // } 174 // default: 175 // break; 176 // } 161 177 } 162 return std::vector<result_t>(1, 0.);178 return results_t(1, result); 163 179 } 164 180 -
src/Potentials/Specifics/ConstantPotential.hpp
r16227a re1fe7e 36 36 friend class PotentialFactory; 37 37 // some repeated typedefs to avoid ambiguities 38 typedef FunctionModel::list_of_arguments_t list_of_arguments_t; 38 39 typedef FunctionModel::arguments_t arguments_t; 39 40 typedef FunctionModel::result_t result_t; … … 90 91 /** Evaluates the harmonic potential function for the given arguments. 91 92 * 92 * @param arguments single distance93 * @param listarguments empty list 93 94 * @return value of the potential function 94 95 */ 95 results_t operator()(const arguments_t &arguments) const;96 results_t operator()(const list_of_arguments_t &listarguments) const; 96 97 97 98 /** Evaluates the derivative of the potential function. 98 99 * 99 * @param arguments single distance100 * @param listarguments empty list 100 101 * @return vector with derivative with respect to the input degrees of freedom 101 102 */ 102 derivative_components_t derivative(const arguments_t &arguments) const;103 derivative_components_t derivative(const list_of_arguments_t &listarguments) const; 103 104 104 105 /** Evaluates the derivative of the function with the given \a arguments 105 106 * with respect to a specific parameter indicated by \a index. 106 107 * 107 * \param arguments set of arguments as input variables to the function108 * \param listarguments empty list 108 109 * \param index derivative of which parameter 109 110 * \return result vector containing the derivative with respect to the given 110 111 * input 111 112 */ 112 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;113 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 113 114 114 115 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/FourBodyPotential_Torsion.cpp
r16227a re1fe7e 140 140 FourBodyPotential_Torsion::results_t 141 141 FourBodyPotential_Torsion::operator()( 142 const arguments_t &arguments142 const list_of_arguments_t &listarguments 143 143 ) const 144 144 { 145 ASSERT( arguments.size() == getSpecificArgumentCount(), 146 "FourBodyPotential_Torsion::operator() - requires exactly three arguments."); 147 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 148 arguments, getParticleTypes()), 149 "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments."); 150 const argument_t &r_ij = arguments[0]; // 01 151 const argument_t &r_ik = arguments[1]; // 02 152 const argument_t &r_il = arguments[2]; // 03 153 const argument_t &r_jk = arguments[3]; // 12 154 const argument_t &r_jl = arguments[4]; // 13 155 const argument_t &r_kl = arguments[5]; // 23 156 const result_t result = 157 params[spring_constant] 158 * Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 ); 159 return std::vector<result_t>(1, result); 145 result_t result = 0.; 146 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 147 iter != listarguments.end(); ++iter) { 148 const arguments_t &arguments = *iter; 149 ASSERT( arguments.size() == getSpecificArgumentCount(), 150 "FourBodyPotential_Torsion::operator() - requires exactly three arguments."); 151 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 152 arguments, getParticleTypes()), 153 "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments."); 154 const argument_t &r_ij = arguments[0]; // 01 155 const argument_t &r_ik = arguments[1]; // 02 156 const argument_t &r_il = arguments[2]; // 03 157 const argument_t &r_jk = arguments[3]; // 12 158 const argument_t &r_jl = arguments[4]; // 13 159 const argument_t &r_kl = arguments[5]; // 23 160 result += 161 params[spring_constant] 162 * Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) 163 - params[equilibrium_distance], 2 ); 164 } 165 return results_t(1, result); 160 166 } 161 167 162 168 FourBodyPotential_Torsion::derivative_components_t 163 169 FourBodyPotential_Torsion::derivative( 164 const arguments_t &arguments170 const list_of_arguments_t &listarguments 165 171 ) const 166 172 { 167 ASSERT( arguments.size() == getSpecificArgumentCount(), 168 "FourBodyPotential_Torsion::operator() - requires exactly three arguments."); 169 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 170 arguments, getParticleTypes()), 171 "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments."); 172 derivative_components_t result; 173 const argument_t &r_ij = arguments[0]; // 01 174 const argument_t &r_ik = arguments[1]; // 02 175 const argument_t &r_il = arguments[2]; // 03 176 const argument_t &r_jk = arguments[3]; // 12 177 const argument_t &r_jl = arguments[4]; // 13 178 const argument_t &r_kl = arguments[5]; // 23 179 result.push_back( 2. * params[spring_constant] * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]) ); 180 ASSERT( result.size() == 1, 181 "FourBodyPotential_Torsion::operator() - we did not create exactly one component."); 182 return result; 173 result_t result = 0.; 174 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 175 iter != listarguments.end(); ++iter) { 176 const arguments_t &arguments = *iter; 177 ASSERT( arguments.size() == getSpecificArgumentCount(), 178 "FourBodyPotential_Torsion::operator() - requires exactly three arguments."); 179 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 180 arguments, getParticleTypes()), 181 "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments."); 182 const argument_t &r_ij = arguments[0]; // 01 183 const argument_t &r_ik = arguments[1]; // 02 184 const argument_t &r_il = arguments[2]; // 03 185 const argument_t &r_jk = arguments[3]; // 12 186 const argument_t &r_jl = arguments[4]; // 13 187 const argument_t &r_kl = arguments[5]; // 23 188 result += 189 2. * params[spring_constant] * 190 ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) 191 - params[equilibrium_distance]); 192 } 193 return derivative_components_t(1, result); 183 194 } 184 195 185 196 FourBodyPotential_Torsion::results_t 186 197 FourBodyPotential_Torsion::parameter_derivative( 187 const arguments_t &arguments,198 const list_of_arguments_t &listarguments, 188 199 const size_t index 189 200 ) const 190 201 { 191 ASSERT( arguments.size() == getSpecificArgumentCount(), 192 "FourBodyPotential_Torsion::parameter_derivative() - requires exactly three arguments."); 193 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 194 arguments, getParticleTypes()), 195 "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments."); 196 const argument_t &r_ij = arguments[0]; // 01 197 const argument_t &r_ik = arguments[1]; // 02 198 const argument_t &r_il = arguments[2]; // 03 199 const argument_t &r_jk = arguments[3]; // 12 200 const argument_t &r_jl = arguments[4]; // 13 201 const argument_t &r_kl = arguments[5]; // 23 202 switch (index) { 203 case spring_constant: 204 { 205 const result_t result = 206 Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 ); 207 return std::vector<result_t>(1, result); 208 break; 202 result_t result = 0.; 203 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 204 iter != listarguments.end(); ++iter) { 205 const arguments_t &arguments = *iter; 206 ASSERT( arguments.size() == getSpecificArgumentCount(), 207 "FourBodyPotential_Torsion::parameter_derivative() - requires exactly three arguments."); 208 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 209 arguments, getParticleTypes()), 210 "FourBodyPotential_Torsion::operator() - types don't match with ones in arguments."); 211 const argument_t &r_ij = arguments[0]; // 01 212 const argument_t &r_ik = arguments[1]; // 02 213 const argument_t &r_il = arguments[2]; // 03 214 const argument_t &r_jk = arguments[3]; // 12 215 const argument_t &r_jl = arguments[4]; // 13 216 const argument_t &r_kl = arguments[5]; // 23 217 switch (index) { 218 case spring_constant: 219 { 220 result += 221 Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 ); 222 break; 223 } 224 case equilibrium_distance: 225 { 226 result += 227 -2. * params[spring_constant] 228 * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]); 229 break; 230 } 231 default: 232 ASSERT(0, "FourBodyPotential_Torsion::parameter_derivative() - derivative to unknown parameter desired."); 233 break; 209 234 } 210 case equilibrium_distance:211 {212 const result_t result =213 -2. * params[spring_constant]214 * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]);215 return std::vector<result_t>(1, result);216 break;217 }218 default:219 ASSERT(0, "FourBodyPotential_Torsion::parameter_derivative() - derivative to unknown parameter desired.");220 break;221 235 } 222 return std::vector<result_t>(1);236 return results_t(1, result); 223 237 } 224 238 -
src/Potentials/Specifics/FourBodyPotential_Torsion.hpp
r16227a re1fe7e 37 37 friend class PotentialFactory; 38 38 // some repeated typedefs to avoid ambiguities 39 typedef FunctionModel::list_of_arguments_t list_of_arguments_t; 39 40 typedef FunctionModel::arguments_t arguments_t; 40 41 typedef FunctionModel::result_t result_t; … … 92 93 /** Evaluates the harmonic potential function for the given arguments. 93 94 * 94 * @param arguments single distance95 * @param listarguments list of six distances 95 96 * @return value of the potential function 96 97 */ 97 results_t operator()(const arguments_t &arguments) const;98 results_t operator()(const list_of_arguments_t &listarguments) const; 98 99 99 100 /** Evaluates the derivative of the potential function. 100 101 * 101 * @param arguments single distance102 * @param listarguments list of six distances 102 103 * @return vector with derivative with respect to the input degrees of freedom 103 104 */ 104 derivative_components_t derivative(const arguments_t &arguments) const;105 derivative_components_t derivative(const list_of_arguments_t &listarguments) const; 105 106 106 107 /** Evaluates the derivative of the function with the given \a arguments 107 108 * with respect to a specific parameter indicated by \a index. 108 109 * 109 * \param arguments set of arguments as input variables to the function110 * \param listarguments list of six distances 110 111 * \param index derivative of which parameter 111 112 * \return result vector containing the derivative with respect to the given 112 113 * input 113 114 */ 114 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;115 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 115 116 116 117 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp
r16227a re1fe7e 184 184 ManyBodyPotential_Tersoff::results_t 185 185 ManyBodyPotential_Tersoff::operator()( 186 const arguments_t &arguments186 const list_of_arguments_t &listarguments 187 187 ) const 188 188 { 189 189 // 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 229 ManyBodyPotential_Tersoff::derivative_components_t 230 ManyBodyPotential_Tersoff::derivative( 231 const list_of_arguments_t &listarguments 232 ) const 233 { 234 // Info info(__func__); 235 return derivative_components_t(); 236 } 237 238 ManyBodyPotential_Tersoff::results_t 239 ManyBodyPotential_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)) 206 301 * 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 *( 211 335 function_prefactor( 212 336 params[beta], … … 217 341 r_ij.distance) 218 342 ); 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; 285 491 } 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."); 306 494 } 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 log395 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]) - temp404 * (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.");487 495 } 488 496 return results_t(1,-result); -
src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp
r16227a re1fe7e 38 38 friend class PotentialFactory; 39 39 // some repeated typedefs to avoid ambiguities 40 typedef FunctionModel::list_of_arguments_t list_of_arguments_t; 40 41 typedef FunctionModel::arguments_t arguments_t; 41 42 typedef FunctionModel::result_t result_t; … … 106 107 /** Evaluates the Tersoff potential for the given arguments. 107 108 * 108 * @param arguments single distance109 * @param listarguments list of distances 109 110 * @return value of the potential function 110 111 */ 111 results_t operator()(const arguments_t &arguments) const;112 results_t operator()(const list_of_arguments_t &listarguments) const; 112 113 113 114 /** Evaluates the derivative of the Tersoff potential with respect to the 114 115 * input variables. 115 116 * 116 * @param arguments single distance117 * @param listarguments list of distances 117 118 * @return vector with components of the derivative 118 119 */ 119 derivative_components_t derivative(const arguments_t &arguments) const;120 derivative_components_t derivative(const list_of_arguments_t &listarguments) const; 120 121 121 122 /** Evaluates the derivative of the function with the given \a arguments 122 123 * with respect to a specific parameter indicated by \a index. 123 124 * 124 * \param arguments set of arguments as input variables to the function125 * \param listarguments list of distances 125 126 * \param index derivative of which parameter 126 127 * \return result vector containing the derivative with respect to the given 127 128 * input 128 129 */ 129 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;130 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 130 131 131 132 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/PairPotential_Harmonic.cpp
r16227a re1fe7e 117 117 PairPotential_Harmonic::results_t 118 118 PairPotential_Harmonic::operator()( 119 const arguments_t &arguments119 const list_of_arguments_t &listarguments 120 120 ) const 121 121 { 122 ASSERT( arguments.size() == 1, 123 "PairPotential_Harmonic::operator() - requires exactly one argument."); 124 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 125 arguments, getParticleTypes()), 126 "PairPotential_Harmonic::operator() - types don't match with ones in arguments."); 127 const argument_t &r_ij = arguments[0]; 128 const result_t result = 129 params[spring_constant] 130 * Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 ); 131 return std::vector<result_t>(1, result); 122 result_t result = 0.; 123 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 124 iter != listarguments.end(); ++iter) { 125 const arguments_t &arguments = *iter; 126 ASSERT( arguments.size() == 1, 127 "PairPotential_Harmonic::operator() - requires exactly one argument."); 128 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 129 arguments, getParticleTypes()), 130 "PairPotential_Harmonic::operator() - types don't match with ones in arguments."); 131 const argument_t &r_ij = arguments[0]; 132 result += 133 params[spring_constant] 134 * Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 ); 135 } 136 return results_t(1, result); 132 137 } 133 138 134 139 PairPotential_Harmonic::derivative_components_t 135 140 PairPotential_Harmonic::derivative( 136 const arguments_t &arguments141 const list_of_arguments_t &listarguments 137 142 ) const 138 143 { 139 ASSERT( arguments.size() == 1, 140 "PairPotential_Harmonic::operator() - requires exactly one argument."); 141 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 142 arguments, getParticleTypes()), 143 "PairPotential_Harmonic::operator() - types don't match with ones in arguments."); 144 derivative_components_t result; 145 const argument_t &r_ij = arguments[0]; 146 result.push_back( 2. * params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]) ); 147 ASSERT( result.size() == 1, 148 "PairPotential_Harmonic::operator() - we did not create exactly one component."); 149 return result; 144 result_t result = 0.; 145 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 146 iter != listarguments.end(); ++iter) { 147 const arguments_t &arguments = *iter; 148 ASSERT( arguments.size() == 1, 149 "PairPotential_Harmonic::operator() - requires exactly one argument."); 150 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 151 arguments, getParticleTypes()), 152 "PairPotential_Harmonic::operator() - types don't match with ones in arguments."); 153 const argument_t &r_ij = arguments[0]; 154 result += 155 2. * params[spring_constant] * 156 ( r_ij.distance - params[equilibrium_distance]); 157 } 158 return derivative_components_t(1, result); 150 159 } 151 160 152 161 PairPotential_Harmonic::results_t 153 162 PairPotential_Harmonic::parameter_derivative( 154 const arguments_t &arguments,163 const list_of_arguments_t &listarguments, 155 164 const size_t index 156 165 ) const 157 166 { 158 ASSERT( arguments.size() == 1, 159 "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument."); 160 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 161 arguments, getParticleTypes()), 162 "PairPotential_Harmonic::operator() - types don't match with ones in arguments."); 163 const argument_t &r_ij = arguments[0]; 164 switch (index) { 165 case spring_constant: 166 { 167 const result_t result = 168 Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 ); 169 return std::vector<result_t>(1, result); 170 break; 167 result_t result = 0.; 168 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 169 iter != listarguments.end(); ++iter) { 170 const arguments_t &arguments = *iter; 171 ASSERT( arguments.size() == 1, 172 "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument."); 173 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 174 arguments, getParticleTypes()), 175 "PairPotential_Harmonic::operator() - types don't match with ones in arguments."); 176 const argument_t &r_ij = arguments[0]; 177 switch (index) { 178 case spring_constant: 179 { 180 result += 181 Helpers::pow( r_ij.distance - params[equilibrium_distance], 2 ); 182 break; 183 } 184 case equilibrium_distance: 185 { 186 result += 187 -2. * params[spring_constant] 188 * ( r_ij.distance - params[equilibrium_distance]); 189 break; 190 } 191 default: 192 ASSERT(0, "PairPotential_Harmonic::parameter_derivative() - derivative to unknown parameter desired."); 193 break; 171 194 } 172 case equilibrium_distance:173 {174 const result_t result =175 -2. * params[spring_constant]176 * ( r_ij.distance - params[equilibrium_distance]);177 return std::vector<result_t>(1, result);178 break;179 }180 default:181 ASSERT(0, "PairPotential_Harmonic::parameter_derivative() - derivative to unknown parameter desired.");182 break;183 195 } 184 185 return PairPotential_Harmonic::results_t(1, 0.); 196 return results_t(1, result); 186 197 } 187 198 -
src/Potentials/Specifics/PairPotential_Harmonic.hpp
r16227a re1fe7e 35 35 friend class PotentialFactory; 36 36 // some repeated typedefs to avoid ambiguities 37 typedef FunctionModel::list_of_arguments_t list_of_arguments_t; 37 38 typedef FunctionModel::arguments_t arguments_t; 38 39 typedef FunctionModel::result_t result_t; … … 90 91 /** Evaluates the harmonic potential function for the given arguments. 91 92 * 92 * @param arguments single distance93 * @param listarguments list of single distances 93 94 * @return value of the potential function 94 95 */ 95 results_t operator()(const arguments_t &arguments) const;96 results_t operator()(const list_of_arguments_t &listarguments) const; 96 97 97 98 /** Evaluates the derivative of the potential function. 98 99 * 99 * @param arguments single distance100 * @param listarguments list of single distances 100 101 * @return vector with derivative with respect to the input degrees of freedom 101 102 */ 102 derivative_components_t derivative(const arguments_t &arguments) const;103 derivative_components_t derivative(const list_of_arguments_t &listarguments) const; 103 104 104 105 /** Evaluates the derivative of the function with the given \a arguments 105 106 * with respect to a specific parameter indicated by \a index. 106 107 * 107 * \param arguments set of arguments as input variables to the function108 * \param listarguments list of single distances 108 109 * \param index derivative of which parameter 109 110 * \return result vector containing the derivative with respect to the given 110 111 * input 111 112 */ 112 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;113 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 113 114 114 115 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/PairPotential_LennardJones.cpp
r16227a re1fe7e 123 123 PairPotential_LennardJones::results_t 124 124 PairPotential_LennardJones::operator()( 125 const arguments_t &arguments125 const list_of_arguments_t &listarguments 126 126 ) const 127 127 { 128 ASSERT( arguments.size() == 1, 129 "PairPotential_LennardJones::operator() - requires one argument."); 130 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 131 arguments, getParticleTypes()), 132 "PairPotential_LennardJones::operator() - types don't match with ones in arguments."); 133 const double &r = arguments[0].distance; 134 const double temp = Helpers::pow(params[sigma]/r, 6); 135 const result_t result = 4.*params[epsilon] * (temp*temp - temp); 136 return std::vector<result_t>(1, result); 128 result_t result = 0.; 129 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 130 iter != listarguments.end(); ++iter) { 131 const arguments_t &arguments = *iter; 132 ASSERT( arguments.size() == 1, 133 "PairPotential_LennardJones::operator() - requires one argument."); 134 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 135 arguments, getParticleTypes()), 136 "PairPotential_LennardJones::operator() - types don't match with ones in arguments."); 137 const double &r = arguments[0].distance; 138 const double temp = Helpers::pow(params[sigma]/r, 6); 139 result += 4.*params[epsilon] * (temp*temp - temp); 140 } 141 return results_t(1, result); 137 142 } 138 143 139 144 PairPotential_LennardJones::derivative_components_t 140 145 PairPotential_LennardJones::derivative( 141 const arguments_t &arguments146 const list_of_arguments_t &listarguments 142 147 ) const 143 148 { 144 ASSERT( arguments.size() == 1,145 "PairPotential_LennardJones::operator() - requires no argument.");146 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(147 arguments, getParticleTypes()),148 "PairPotential_LennardJones::operator() - types don't match with ones in arguments.");149 const double &r = arguments[0].distance;150 149 const double sigma6 = Helpers::pow(params[sigma], 6); 151 const result_t result = 152 4.*params[epsilon] * ( 153 sigma6*sigma6*(-12.) / Helpers::pow(r,13) 154 - sigma6*(-6.) /Helpers::pow(r,7) 155 ); 156 derivative_components_t results(1, result); 157 return results; 150 result_t result = 0.; 151 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 152 iter != listarguments.end(); ++iter) { 153 const arguments_t &arguments = *iter; 154 ASSERT( arguments.size() == 1, 155 "PairPotential_LennardJones::operator() - requires no argument."); 156 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 157 arguments, getParticleTypes()), 158 "PairPotential_LennardJones::operator() - types don't match with ones in arguments."); 159 const double &r = arguments[0].distance; 160 result += 161 4.*params[epsilon] * ( 162 sigma6*sigma6*(-12.) / Helpers::pow(r,13) 163 - sigma6*(-6.) /Helpers::pow(r,7) 164 ); 165 } 166 return derivative_components_t(1, result); 158 167 } 159 168 160 169 PairPotential_LennardJones::results_t 161 170 PairPotential_LennardJones::parameter_derivative( 162 const arguments_t &arguments,171 const list_of_arguments_t &listarguments, 163 172 const size_t index 164 173 ) const 165 174 { 166 ASSERT( arguments.size() == 1, 167 "PairPotential_LennardJones::parameter_derivative() - requires no argument."); 168 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 169 arguments, getParticleTypes()), 170 "PairPotential_LennardJones::operator() - types don't match with ones in arguments."); 171 const double &r = arguments[0].distance; 172 switch (index) { 173 case epsilon: 174 { 175 const double temp = Helpers::pow(params[sigma]/r, 6); 176 const result_t result = 4. * (temp*temp - temp); 177 return std::vector<result_t>(1, result); 178 break; 175 result_t result = 0.; 176 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 177 iter != listarguments.end(); ++iter) { 178 const arguments_t &arguments = *iter; 179 ASSERT( arguments.size() == 1, 180 "PairPotential_LennardJones::parameter_derivative() - requires no argument."); 181 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 182 arguments, getParticleTypes()), 183 "PairPotential_LennardJones::operator() - types don't match with ones in arguments."); 184 const double &r = arguments[0].distance; 185 switch (index) { 186 case epsilon: 187 { 188 const double temp = Helpers::pow(params[sigma]/r, 6); 189 result += 4. * (temp*temp - temp); 190 break; 191 } 192 case sigma: 193 { 194 const double r6 = Helpers::pow(r, 6); 195 result += 196 4.*params[epsilon] * ( 197 12. * Helpers::pow(params[sigma],11)/(r6*r6) 198 - 6. * Helpers::pow(params[sigma],5)/r6 199 ); 200 break; 201 } 202 default: 203 break; 179 204 } 180 case sigma:181 {182 const double r6 = Helpers::pow(r, 6);183 const result_t result =184 4.*params[epsilon] * (185 12. * Helpers::pow(params[sigma],11)/(r6*r6)186 - 6. * Helpers::pow(params[sigma],5)/r6187 );188 return std::vector<result_t>(1, result);189 break;190 }191 default:192 break;193 205 } 194 return std::vector<result_t>(1, 0.);206 return results_t(1, result); 195 207 } 196 208 -
src/Potentials/Specifics/PairPotential_LennardJones.hpp
r16227a re1fe7e 36 36 friend class PotentialFactory; 37 37 // some repeated typedefs to avoid ambiguities 38 typedef FunctionModel::list_of_arguments_t list_of_arguments_t; 38 39 typedef FunctionModel::arguments_t arguments_t; 39 40 typedef FunctionModel::result_t result_t; … … 92 93 /** Evaluates the harmonic potential function for the given arguments. 93 94 * 94 * @param arguments single distance95 * @param listarguments list of single distances 95 96 * @return value of the potential function 96 97 */ 97 results_t operator()(const arguments_t &arguments) const;98 results_t operator()(const list_of_arguments_t &listarguments) const; 98 99 99 100 /** Evaluates the derivative of the potential function. 100 101 * 101 * @param arguments single distance102 * @param listarguments list of single distances 102 103 * @return vector with derivative with respect to the input degrees of freedom 103 104 */ 104 derivative_components_t derivative(const arguments_t &arguments) const;105 derivative_components_t derivative(const list_of_arguments_t &listarguments) const; 105 106 106 107 /** Evaluates the derivative of the function with the given \a arguments 107 108 * with respect to a specific parameter indicated by \a index. 108 109 * 109 * \param arguments set of arguments as input variables to the function110 * \param listarguments list of single distances 110 111 * \param index derivative of which parameter 111 112 * \return result vector containing the derivative with respect to the given 112 113 * input 113 114 */ 114 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;115 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 115 116 116 117 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/PairPotential_Morse.cpp
r16227a re1fe7e 124 124 PairPotential_Morse::results_t 125 125 PairPotential_Morse::operator()( 126 const arguments_t &arguments126 const list_of_arguments_t &listarguments 127 127 ) const 128 128 { 129 ASSERT( arguments.size() == 1, 130 "PairPotential_Morse::operator() - requires exactly one argument."); 131 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 132 arguments, getParticleTypes()), 133 "PairPotential_Morse::operator() - types don't match with ones in arguments."); 134 const argument_t &r_ij = arguments[0]; 135 // Maple: f(r,D,k,R,c) := D * (1 - exp(-k*(r-R)))^(2)+c 136 const result_t result = 137 params[dissociation_energy] * Helpers::pow( 1. 138 - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2); 129 result_t result = 0.; 130 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 131 iter != listarguments.end(); ++iter) { 132 const arguments_t &arguments = *iter; 133 ASSERT( arguments.size() == 1, 134 "PairPotential_Morse::operator() - requires exactly one argument."); 135 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 136 arguments, getParticleTypes()), 137 "PairPotential_Morse::operator() - types don't match with ones in arguments."); 138 const argument_t &r_ij = arguments[0]; 139 // Maple: f(r,D,k,R,c) := D * (1 - exp(-k*(r-R)))^(2)+c 140 result += 141 params[dissociation_energy] * Helpers::pow( 1. 142 - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2); 143 } 139 144 return std::vector<result_t>(1, result); 140 145 } … … 142 147 PairPotential_Morse::derivative_components_t 143 148 PairPotential_Morse::derivative( 144 const arguments_t &arguments149 const list_of_arguments_t &listarguments 145 150 ) const 146 151 { 147 ASSERT( arguments.size() == 1, 148 "PairPotential_Morse::operator() - requires exactly one argument."); 149 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 150 arguments, getParticleTypes()), 151 "PairPotential_Morse::operator() - types don't match with ones in arguments."); 152 derivative_components_t result; 153 const argument_t &r_ij = arguments[0]; 154 // Maple result: 2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R)) 155 result.push_back( 156 2. * params[dissociation_energy] 157 * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))) 158 * (- params[spring_constant]) * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])) 159 ); 160 ASSERT( result.size() == 1, 161 "PairPotential_Morse::operator() - we did not create exactly one component."); 162 return result; 152 result_t result = 0.; 153 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 154 iter != listarguments.end(); ++iter) { 155 const arguments_t &arguments = *iter; 156 ASSERT( arguments.size() == 1, 157 "PairPotential_Morse::operator() - requires exactly one argument."); 158 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 159 arguments, getParticleTypes()), 160 "PairPotential_Morse::operator() - types don't match with ones in arguments."); 161 const argument_t &r_ij = arguments[0]; 162 // Maple result: 2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R)) 163 result += 164 2. * params[dissociation_energy] 165 * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))) 166 * (- params[spring_constant]) * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])); 167 } 168 return derivative_components_t(1, result); 163 169 } 164 170 165 171 PairPotential_Morse::results_t 166 172 PairPotential_Morse::parameter_derivative( 167 const arguments_t &arguments,173 const list_of_arguments_t &listarguments, 168 174 const size_t index 169 175 ) const 170 176 { 171 ASSERT( arguments.size() == 1, 172 "PairPotential_Morse::parameter_derivative() - requires exactly one argument."); 173 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 174 arguments, getParticleTypes()), 175 "PairPotential_Morse::operator() - types don't match with ones in arguments."); 176 const argument_t &r_ij = arguments[0]; 177 switch (index) { 178 case spring_constant: 179 { 180 // Maple result: -2*D*(1-exp(-k*(r-R)))*(-r+R)*exp(-k*(r-R)) 181 const result_t result = 182 - 2. * params[dissociation_energy] 183 * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))) 184 * (- r_ij.distance + params[equilibrium_distance]) 185 * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])) 186 ; 187 return std::vector<result_t>(1, result); 188 break; 177 result_t result = 0.; 178 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 179 iter != listarguments.end(); ++iter) { 180 const arguments_t &arguments = *iter; 181 ASSERT( arguments.size() == 1, 182 "PairPotential_Morse::parameter_derivative() - requires exactly one argument."); 183 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 184 arguments, getParticleTypes()), 185 "PairPotential_Morse::operator() - types don't match with ones in arguments."); 186 const argument_t &r_ij = arguments[0]; 187 switch (index) { 188 case spring_constant: 189 { 190 // Maple result: -2*D*(1-exp(-k*(r-R)))*(-r+R)*exp(-k*(r-R)) 191 result += 192 - 2. * params[dissociation_energy] 193 * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))) 194 * (- r_ij.distance + params[equilibrium_distance]) 195 * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])) 196 ; 197 break; 198 } 199 case equilibrium_distance: 200 { 201 // Maple result: -2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R)) 202 result += 203 - 2. * params[dissociation_energy] 204 * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))) 205 * params[spring_constant] * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])) 206 ; 207 break; 208 } 209 case dissociation_energy: 210 { 211 // Maple result: (1-exp(-k*(r-R)))^2 212 result += 213 Helpers::pow(1. 214 - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2); 215 break; 216 } 217 default: 218 ASSERT(0, "PairPotential_Morse::parameter_derivative() - derivative to unknown parameter desired."); 219 break; 189 220 } 190 case equilibrium_distance:191 {192 // Maple result: -2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))193 const result_t result =194 - 2. * params[dissociation_energy]195 * ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))196 * params[spring_constant] * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))197 ;198 return std::vector<result_t>(1, result);199 break;200 }201 case dissociation_energy:202 {203 // Maple result: (1-exp(-k*(r-R)))^2204 const result_t result =205 Helpers::pow(1.206 - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);207 return std::vector<result_t>(1, result);208 break;209 }210 default:211 ASSERT(0, "PairPotential_Morse::parameter_derivative() - derivative to unknown parameter desired.");212 break;213 221 } 214 return std::vector<result_t>(1, 0.);222 return results_t(1, result); 215 223 } 216 224 -
src/Potentials/Specifics/PairPotential_Morse.hpp
r16227a re1fe7e 35 35 friend class PotentialFactory; 36 36 // some repeated typedefs to avoid ambiguities 37 typedef FunctionModel::list_of_arguments_t list_of_arguments_t; 37 38 typedef FunctionModel::arguments_t arguments_t; 38 39 typedef FunctionModel::result_t result_t; … … 91 92 /** Evaluates the harmonic potential function for the given arguments. 92 93 * 93 * @param arguments single distance94 * @param listarguments list of single distances 94 95 * @return value of the potential function 95 96 */ 96 results_t operator()(const arguments_t &arguments) const;97 results_t operator()(const list_of_arguments_t &listarguments) const; 97 98 98 99 /** Evaluates the derivative of the potential function. 99 100 * 100 * @param arguments single distance101 * @param listarguments list of single distances 101 102 * @return vector with derivative with respect to the input degrees of freedom 102 103 */ 103 derivative_components_t derivative(const arguments_t &arguments) const;104 derivative_components_t derivative(const list_of_arguments_t &listarguments) const; 104 105 105 106 /** Evaluates the derivative of the function with the given \a arguments 106 107 * with respect to a specific parameter indicated by \a index. 107 108 * 108 * \param arguments set of arguments as input variables to the function109 * \param listarguments list of single distances 109 110 * \param index derivative of which parameter 110 111 * \return result vector containing the derivative with respect to the given 111 112 * input 112 113 */ 113 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;114 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 114 115 115 116 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/ThreeBodyPotential_Angle.cpp
r16227a re1fe7e 136 136 ThreeBodyPotential_Angle::results_t 137 137 ThreeBodyPotential_Angle::operator()( 138 const arguments_t &arguments138 const list_of_arguments_t &listarguments 139 139 ) const 140 140 { 141 ASSERT( arguments.size() == 3, 142 "ThreeBodyPotential_Angle::operator() - requires exactly three arguments."); 143 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 144 arguments, getParticleTypes()), 145 "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments."); 146 const argument_t &r_ij = arguments[0]; // 01 147 const argument_t &r_jk = arguments[2]; // 12 148 const argument_t &r_ik = arguments[1]; // 02 149 const result_t result = 150 params[spring_constant] 151 * Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 ); 152 return std::vector<result_t>(1, result); 141 result_t result = 0.; 142 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 143 iter != listarguments.end(); ++iter) { 144 const arguments_t &arguments = *iter; 145 ASSERT( arguments.size() == 3, 146 "ThreeBodyPotential_Angle::operator() - requires exactly three arguments."); 147 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 148 arguments, getParticleTypes()), 149 "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments."); 150 const argument_t &r_ij = arguments[0]; // 01 151 const argument_t &r_jk = arguments[2]; // 12 152 const argument_t &r_ik = arguments[1]; // 02 153 result += 154 params[spring_constant] 155 * Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) 156 - params[equilibrium_distance], 2 ); 157 } 158 return results_t(1, result); 153 159 } 154 160 155 161 ThreeBodyPotential_Angle::derivative_components_t 156 162 ThreeBodyPotential_Angle::derivative( 157 const arguments_t &arguments163 const list_of_arguments_t &listarguments 158 164 ) const 159 165 { 160 ASSERT( arguments.size() == 3, 161 "ThreeBodyPotential_Angle::operator() - requires exactly three arguments."); 162 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 163 arguments, getParticleTypes()), 164 "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments."); 165 derivative_components_t result; 166 const argument_t &r_ij = arguments[0]; //01 167 const argument_t &r_jk = arguments[2]; //12 168 const argument_t &r_ik = arguments[1]; //02 169 result.push_back( 2. * params[spring_constant] * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]) ); 170 ASSERT( result.size() == 1, 171 "ThreeBodyPotential_Angle::operator() - we did not create exactly one component."); 172 return result; 166 result_t result = 0.; 167 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 168 iter != listarguments.end(); ++iter) { 169 const arguments_t &arguments = *iter; 170 ASSERT( arguments.size() == 3, 171 "ThreeBodyPotential_Angle::operator() - requires exactly three arguments."); 172 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 173 arguments, getParticleTypes()), 174 "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments."); 175 const argument_t &r_ij = arguments[0]; //01 176 const argument_t &r_jk = arguments[2]; //12 177 const argument_t &r_ik = arguments[1]; //02 178 result += 179 2. * params[spring_constant] * 180 ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) 181 - params[equilibrium_distance]); 182 } 183 return derivative_components_t(1, result); 173 184 } 174 185 175 186 ThreeBodyPotential_Angle::results_t 176 187 ThreeBodyPotential_Angle::parameter_derivative( 177 const arguments_t &arguments,188 const list_of_arguments_t &listarguments, 178 189 const size_t index 179 190 ) const 180 191 { 181 ASSERT( arguments.size() == 3, 182 "ThreeBodyPotential_Angle::parameter_derivative() - requires exactly three arguments."); 183 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 184 arguments, getParticleTypes()), 185 "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments."); 186 const argument_t &r_ij = arguments[0]; //01 187 const argument_t &r_jk = arguments[2]; //12 188 const argument_t &r_ik = arguments[1]; //02 189 switch (index) { 190 case spring_constant: 191 { 192 const result_t result = 193 Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 ); 194 return std::vector<result_t>(1, result); 195 break; 192 result_t result = 0.; 193 for(list_of_arguments_t::const_iterator iter = listarguments.begin(); 194 iter != listarguments.end(); ++iter) { 195 const arguments_t &arguments = *iter; 196 ASSERT( arguments.size() == 3, 197 "ThreeBodyPotential_Angle::parameter_derivative() - requires exactly three arguments."); 198 ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes( 199 arguments, getParticleTypes()), 200 "ThreeBodyPotential_Angle::operator() - types don't match with ones in arguments."); 201 const argument_t &r_ij = arguments[0]; //01 202 const argument_t &r_jk = arguments[2]; //12 203 const argument_t &r_ik = arguments[1]; //02 204 switch (index) { 205 case spring_constant: 206 { 207 result += 208 Helpers::pow( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance], 2 ); 209 break; 210 } 211 case equilibrium_distance: 212 { 213 result += 214 -2. * params[spring_constant] 215 * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]); 216 break; 217 } 218 default: 219 ASSERT(0, "ThreeBodyPotential_Angle::parameter_derivative() - derivative to unknown parameter desired."); 220 break; 196 221 } 197 case equilibrium_distance:198 {199 const result_t result =200 -2. * params[spring_constant]201 * ( function_theta(r_ij.distance, r_jk.distance, r_ik.distance) - params[equilibrium_distance]);202 return std::vector<result_t>(1, result);203 break;204 }205 default:206 ASSERT(0, "ThreeBodyPotential_Angle::parameter_derivative() - derivative to unknown parameter desired.");207 break;208 222 } 209 return std::vector<result_t>(1);223 return results_t(1, result); 210 224 } 211 225 -
src/Potentials/Specifics/ThreeBodyPotential_Angle.hpp
r16227a re1fe7e 35 35 friend class PotentialFactory; 36 36 // some repeated typedefs to avoid ambiguities 37 typedef FunctionModel::list_of_arguments_t list_of_arguments_t; 37 38 typedef FunctionModel::arguments_t arguments_t; 38 39 typedef FunctionModel::result_t result_t; … … 90 91 /** Evaluates the harmonic potential function for the given arguments. 91 92 * 92 * @param arguments single distance93 * @param listarguments list of three distances 93 94 * @return value of the potential function 94 95 */ 95 results_t operator()(const arguments_t &arguments) const;96 results_t operator()(const list_of_arguments_t &listarguments) const; 96 97 97 98 /** Evaluates the derivative of the potential function. 98 99 * 99 * @param arguments single distance100 * @param listarguments list of three distances 100 101 * @return vector with derivative with respect to the input degrees of freedom 101 102 */ 102 derivative_components_t derivative(const arguments_t &arguments) const;103 derivative_components_t derivative(const list_of_arguments_t &listarguments) const; 103 104 104 105 /** Evaluates the derivative of the function with the given \a arguments 105 106 * with respect to a specific parameter indicated by \a index. 106 107 * 107 * \param arguments set of arguments as input variables to the function108 * \param listarguments list of three distances 108 109 * \param index derivative of which parameter 109 110 * \return result vector containing the derivative with respect to the given 110 111 * input 111 112 */ 112 results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;113 results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const; 113 114 114 115 /** Returns the functor that converts argument_s into the -
src/Potentials/Specifics/unittests/ConstantPotentialUnitTest.cpp
r16227a re1fe7e 82 82 Helpers::isEqual( 83 83 offset, 84 constant( FunctionModel:: arguments_t() )[0]84 constant( FunctionModel::list_of_arguments_t() )[0] 85 85 ) 86 86 ); … … 97 97 0., 98 98 constant.derivative( 99 FunctionModel:: arguments_t()99 FunctionModel::list_of_arguments_t() 100 100 )[0] 101 101 ) … … 114 114 1., 115 115 constant.parameter_derivative( 116 FunctionModel:: arguments_t(),116 FunctionModel::list_of_arguments_t(), 117 117 0 118 118 )[0] -
src/Potentials/Specifics/unittests/FourBodyPotential_ImproperUnitTest.cpp
r16227a re1fe7e 121 121 for (size_t index = 0; index < input.size(); ++index) { 122 122 const FourBodyPotential_Improper::results_t result = 123 angle( input[index]);123 angle( FunctionModel::list_of_arguments_t(1, input[index]) ); 124 124 CPPUNIT_ASSERT( 125 125 Helpers::isEqual( … … 144 144 0., 145 145 angle.derivative( 146 input[5]146 FunctionModel::list_of_arguments_t(1, input[5]) 147 147 )[0], 148 148 10. … … 165 165 0., 166 166 angle.parameter_derivative( 167 input[5],167 FunctionModel::list_of_arguments_t(1, input[5]), 168 168 0 169 169 )[0], … … 175 175 0., 176 176 angle.parameter_derivative( 177 input[5],177 FunctionModel::list_of_arguments_t(1, input[5]), 178 178 1 179 179 )[0], -
src/Potentials/Specifics/unittests/FourBodyPotential_TorsionUnitTest.cpp
r16227a re1fe7e 121 121 for (size_t index = 0; index < input.size(); ++index) { 122 122 const FourBodyPotential_Torsion::results_t result = 123 angle( input[index]);123 angle( FunctionModel::list_of_arguments_t(1, input[index]) ); 124 124 CPPUNIT_ASSERT( 125 125 Helpers::isEqual( … … 144 144 0., 145 145 angle.derivative( 146 input[5]146 FunctionModel::list_of_arguments_t(1, input[5]) 147 147 )[0], 148 148 10. … … 165 165 0., 166 166 angle.parameter_derivative( 167 input[5],167 FunctionModel::list_of_arguments_t(1, input[5]), 168 168 0 169 169 )[0], … … 175 175 0., 176 176 angle.parameter_derivative( 177 input[5],177 FunctionModel::list_of_arguments_t(1, input[5]), 178 178 1 179 179 )[0], -
src/Potentials/Specifics/unittests/ManyBodyPotential_TersoffUnitTest.cpp
r16227a re1fe7e 299 299 arg.globalid = index; // this is needed for the triplefunction to the configuration 300 300 FunctionModel::arguments_t args(1,arg); 301 const ManyBodyPotential_Tersoff::results_t res = tersoff(args); 301 const ManyBodyPotential_Tersoff::results_t res = 302 tersoff( FunctionModel::list_of_arguments_t(1, args) ); 302 303 temp += res[0]; 303 304 } … … 336 337 // 0., 337 338 // tersoff.derivative( 338 // FunctionModel:: arguments_t(1,argument_t(1.))339 // FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))) 339 340 // )[0] 340 341 // ) … … 361 362 // 0., 362 363 // tersoff.parameter_derivative( 363 // FunctionModel:: arguments_t(1,argument_t(1.)),364 // FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))), 364 365 // 0 365 366 // )[0] … … 370 371 // 0., 371 372 // tersoff.parameter_derivative( 372 // FunctionModel:: arguments_t(1,argument_t(1.)),373 // FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))), 373 374 // 1 374 375 // )[0] … … 379 380 // 1., 380 381 // tersoff.parameter_derivative( 381 // FunctionModel:: arguments_t(1,argument_t(1.)),382 // FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,argument_t(1.))), 382 383 // 2 383 384 // )[0] -
src/Potentials/Specifics/unittests/PairPotential_HarmonicUnitTest.cpp
r16227a re1fe7e 102 102 Helpers::isEqual( 103 103 output[index], 104 harmonic( FunctionModel:: arguments_t(1,arg) )[0]104 harmonic( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0] 105 105 ) 106 106 ); … … 122 122 0., 123 123 harmonic.derivative( 124 FunctionModel:: arguments_t(1,arg)124 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) 125 125 )[0] 126 126 ) … … 143 143 0., 144 144 harmonic.parameter_derivative( 145 FunctionModel:: arguments_t(1,arg),145 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)), 146 146 0 147 147 )[0] … … 152 152 0., 153 153 harmonic.parameter_derivative( 154 FunctionModel:: arguments_t(1,arg),154 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)), 155 155 1 156 156 )[0] -
src/Potentials/Specifics/unittests/PairPotential_LennardJonesUnitTest.cpp
r16227a re1fe7e 118 118 Helpers::isEqual( 119 119 output[index], 120 lj( FunctionModel:: arguments_t(1,arg) )[0],120 lj( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0], 121 121 1.e-4/std::numeric_limits<double>::epsilon() // only compare four digits 122 122 ) … … 140 140 0., 141 141 lj.derivative( 142 FunctionModel:: arguments_t(1,arg)142 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) 143 143 )[0], 144 144 1.e+6 … … 162 162 -1., 163 163 lj.parameter_derivative( 164 FunctionModel:: arguments_t(1,arg),164 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)), 165 165 0 166 166 )[0], … … 172 172 0., 173 173 lj.parameter_derivative( 174 FunctionModel:: arguments_t(1,arg),174 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)), 175 175 1 176 176 )[0], -
src/Potentials/Specifics/unittests/PairPotential_MorseUnitTest.cpp
r16227a re1fe7e 120 120 Helpers::isEqual( 121 121 output[index], 122 Morse( FunctionModel:: arguments_t(1,arg) )[0],122 Morse( FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) )[0], 123 123 1.e-4/std::numeric_limits<double>::epsilon() // only compare four digits 124 124 ) … … 141 141 0., 142 142 Morse.derivative( 143 FunctionModel:: arguments_t(1,arg)143 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)) 144 144 )[0], 145 145 1.e+6 … … 163 163 0., 164 164 Morse.parameter_derivative( 165 FunctionModel:: arguments_t(1,arg),165 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)), 166 166 0 167 167 )[0], … … 173 173 0., 174 174 Morse.parameter_derivative( 175 FunctionModel:: arguments_t(1,arg),175 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)), 176 176 1 177 177 )[0], … … 183 183 0., 184 184 Morse.parameter_derivative( 185 FunctionModel:: arguments_t(1,arg),185 FunctionModel::list_of_arguments_t(1, FunctionModel::arguments_t(1,arg)), 186 186 2 187 187 )[0], -
src/Potentials/Specifics/unittests/ThreeBodyPotential_AngleUnitTest.cpp
r16227a re1fe7e 112 112 for (size_t index = 0; index < input.size(); ++index) { 113 113 const ThreeBodyPotential_Angle::results_t result = 114 angle( input[index]);114 angle( FunctionModel::list_of_arguments_t(1, input[index]) ); 115 115 CPPUNIT_ASSERT( 116 116 Helpers::isEqual( … … 135 135 0., 136 136 angle.derivative( 137 input[5]137 FunctionModel::list_of_arguments_t(1, input[5]) 138 138 )[0], 139 139 10. … … 156 156 0., 157 157 angle.parameter_derivative( 158 input[5],158 FunctionModel::list_of_arguments_t(1, input[5]), 159 159 0 160 160 )[0], … … 166 166 0., 167 167 angle.parameter_derivative( 168 input[5],168 FunctionModel::list_of_arguments_t(1, input[5]), 169 169 1 170 170 )[0], -
src/Potentials/helpers.hpp
r16227a re1fe7e 108 108 } 109 109 110 inline FunctionModel::list_of_arguments_t 111 returnEmptyListArguments() 112 { 113 return FunctionModel::list_of_arguments_t(); 114 } 115 110 116 }; /* namespace Helpers */ 111 117 -
src/Potentials/unittests/CompoundPotentialUnitTest.cpp
r16227a re1fe7e 159 159 argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]); 160 160 argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), input[index]); 161 FunctionModel:: arguments_targs;162 args += firstarg,secondarg;163 const double result = compound( args )[0];161 FunctionModel::list_of_arguments_t listargs; 162 listargs += FunctionModel::arguments_t(1,firstarg),FunctionModel::arguments_t(1,secondarg); 163 const double result = compound( listargs )[0]; 164 164 CPPUNIT_ASSERT( 165 165 Helpers::isEqual( … … 205 205 argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c); 206 206 argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), c); 207 FunctionModel:: arguments_targs;208 args += firstarg,secondarg;209 { 210 const double result = 211 compound.parameter_derivative( 212 args,207 FunctionModel::list_of_arguments_t listargs; 208 listargs += FunctionModel::arguments_t(1,firstarg),FunctionModel::arguments_t(1,secondarg); 209 { 210 const double result = 211 compound.parameter_derivative( 212 listargs, 213 213 0)[0]; 214 214 CPPUNIT_ASSERT( … … 223 223 const double result = 224 224 compound.parameter_derivative( 225 args,225 listargs, 226 226 1)[0]; 227 227 CPPUNIT_ASSERT( … … 236 236 const double result = 237 237 compound.parameter_derivative( 238 args,238 listargs, 239 239 2)[0]; 240 240 CPPUNIT_ASSERT( … … 249 249 const double result = 250 250 compound.parameter_derivative( 251 args,251 listargs, 252 252 3)[0]; 253 253 CPPUNIT_ASSERT(
Note:
See TracChangeset
for help on using the changeset viewer.