source: src/Potentials/Specifics/ManyBodyPotential_Tersoff.cpp@ 55fe788

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
Last change on this file since 55fe788 was 55fe788, checked in by Frederik Heber <heber@…>, 12 years ago

All specific Potentials now assert that arguments have charges matching with internal ParticleTypes_t.

  • this ensures that we don't have to pick indices.first == 0 or something similar spooky anymore.
  • TESTFIX: Modified Potentials/Specifics/unittests that now need to specify the (correct) types in each argument_t.
  • Property mode set to 100644
File size: 20.2 KB
Line 
1/*
2 * Project: MoleCuilder
3 * Description: creates and alters molecular systems
4 * Copyright (C) 2012 University of Bonn. All rights reserved.
5 * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
6 *
7 *
8 * This file is part of MoleCuilder.
9 *
10 * MoleCuilder is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 2 of the License, or
13 * (at your option) any later version.
14 *
15 * MoleCuilder is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
22 */
23
24/*
25 * ManyBodyPotential_Tersoff.cpp
26 *
27 * Created on: Sep 26, 2012
28 * Author: heber
29 */
30
31
32// include config.h
33#ifdef HAVE_CONFIG_H
34#include <config.h>
35#endif
36
37#include "CodePatterns/MemDebug.hpp"
38
39#include "ManyBodyPotential_Tersoff.hpp"
40
41#include <boost/assign/list_of.hpp> // for 'map_list_of()'
42#include <boost/bind.hpp>
43#include <cmath>
44#include <string>
45
46#include "CodePatterns/Assert.hpp"
47//#include "CodePatterns/Info.hpp"
48#include "CodePatterns/Log.hpp"
49
50#include "Potentials/helpers.hpp"
51
52// static definitions
53const ManyBodyPotential_Tersoff::ParameterNames_t
54ManyBodyPotential_Tersoff::ParameterNames =
55 boost::assign::list_of<std::string>
56 ("A")
57 ("B")
58 ("lambda")
59 ("mu")
60 ("beta")
61 ("n")
62 ("c")
63 ("d")
64 ("h")
65 ("offset")
66// ("R")
67// ("S")
68// ("lambda3")
69// ("alpha")
70// ("chi")
71// ("omega")
72 ;
73const std::string ManyBodyPotential_Tersoff::potential_token("tersoff");
74
75ManyBodyPotential_Tersoff::ManyBodyPotential_Tersoff(
76 const ParticleTypes_t &_ParticleTypes,
77 boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction
78 ) :
79 SerializablePotential(_ParticleTypes),
80 params(parameters_t(MAXPARAMS, 0.)),
81 R(3.2),
82 S(3.5),
83 lambda3(0.),
84 alpha(0.),
85 chi(1.),
86 omega(1.),
87 triplefunction(_triplefunction)
88{}
89
90ManyBodyPotential_Tersoff::ManyBodyPotential_Tersoff(
91 const ParticleTypes_t &_ParticleTypes,
92 const double &_R,
93 const double &_S,
94 const double &_A,
95 const double &_B,
96 const double &_lambda,
97 const double &_mu,
98 const double &_lambda3,
99 const double &_alpha,
100 const double &_beta,
101 const double &_chi,
102 const double &_omega,
103 const double &_n,
104 const double &_c,
105 const double &_d,
106 const double &_h,
107 const double &_offset,
108 boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
109 SerializablePotential(_ParticleTypes),
110 params(parameters_t(MAXPARAMS, 0.)),
111 R(_R),
112 S(_S),
113 lambda3(_lambda3),
114 alpha(_alpha),
115 chi(_chi),
116 omega(_mu),
117 triplefunction(_triplefunction)
118{
119// Info info(__func__);
120// R = _R;
121// S = _S;
122 params[A] = _A;
123 params[B] = _B;
124 params[lambda] = _lambda;
125 params[mu] = _mu;
126// lambda3 = _lambda3;
127// alpha = _alpha;
128 params[beta] = _beta;
129// chi = _chi;
130// omega = _omega;
131 params[n] = _n;
132 params[c] = _c;
133 params[d] = _d;
134 params[h] = _h;
135 params[offset] = _offset;
136}
137
138void ManyBodyPotential_Tersoff::setParameters(const parameters_t &_params)
139{
140 const size_t paramsDim = _params.size();
141 ASSERT( paramsDim <= getParameterDimension(),
142 "ManyBodyPotential_Tersoff::setParameters() - we need not more than "
143 +toString(getParameterDimension())+" parameters.");
144 for (size_t i=0; i< paramsDim; ++i)
145 params[i] = _params[i];
146
147#ifndef NDEBUG
148 parameters_t check_params(getParameters());
149 check_params.resize(paramsDim); // truncate to same size
150 ASSERT( check_params == _params,
151 "ManyBodyPotential_Tersoff::setParameters() - failed, mismatch in to be set "
152 +toString(_params)+" and set "+toString(check_params)+" params.");
153#endif
154}
155
156ManyBodyPotential_Tersoff::results_t
157ManyBodyPotential_Tersoff::operator()(
158 const arguments_t &arguments
159 ) const
160{
161// Info info(__func__);
162 double result = 0.;
163 for(arguments_t::const_iterator argiter = arguments.begin();
164 argiter != arguments.end();
165 ++argiter) {
166 const argument_t &r_ij = *argiter;
167 ASSERT( arguments[0].types.first == getParticleTypes()[0],
168 "ManyBodyPotential_Tersoff::operator() - first charge "
169 +toString(arguments[0].types.first)+" is not "+toString(getParticleTypes()[0]));
170 ASSERT( arguments[0].types.second == getParticleTypes()[1],
171 "ManyBodyPotential_Tersoff::operator() - second charge "
172 +toString(arguments[0].types.second)+" is not "+toString(getParticleTypes()[1]));
173
174 const double cutoff = function_cutoff(r_ij.distance);
175 const double temp = (cutoff == 0.) ?
176 0. :
177 cutoff * (
178 function_prefactor(
179 alpha,
180 function_eta(r_ij))
181 * function_smoother(
182 params[A],
183 params[lambda],
184 r_ij.distance)
185 +
186 function_prefactor(
187 params[beta],
188 function_zeta(r_ij))
189 * function_smoother(
190 -params[B],
191 params[mu],
192 r_ij.distance)
193 );
194 result += temp;
195 }
196// LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
197 return std::vector<result_t>(1, params[offset]+result);
198}
199
200ManyBodyPotential_Tersoff::derivative_components_t
201ManyBodyPotential_Tersoff::derivative(
202 const arguments_t &arguments
203 ) const
204{
205// Info info(__func__);
206 return ManyBodyPotential_Tersoff::derivative_components_t();
207}
208
209ManyBodyPotential_Tersoff::results_t
210ManyBodyPotential_Tersoff::parameter_derivative(
211 const arguments_t &arguments,
212 const size_t index
213 ) const
214{
215// Info info(__func__);
216// ASSERT( arguments.size() == 1,
217// "ManyBodyPotential_Tersoff::parameter_derivative() - requires exactly one argument.");
218 if (index == offset)
219 return results_t(1, 1.);
220
221 double result = 0.;
222 for(arguments_t::const_iterator argiter = arguments.begin();
223 argiter != arguments.end();
224 ++argiter) {
225 const argument_t &r_ij = *argiter;
226 ASSERT( arguments[0].types.first == getParticleTypes()[0],
227 "ManyBodyPotential_Tersoff::operator() - first charge "
228 +toString(arguments[0].types.first)+" is not "+toString(getParticleTypes()[0]));
229 ASSERT( arguments[0].types.second == getParticleTypes()[1],
230 "ManyBodyPotential_Tersoff::operator() - second charge "
231 +toString(arguments[0].types.second)+" is not "+toString(getParticleTypes()[1]));
232
233 switch (index) {
234// case R:
235// {
236// result += 0.;
237// break;
238// }
239// case S:
240// {
241// result += 0.;
242// break;
243// }
244 case A:
245 {
246 const double cutoff = function_cutoff(r_ij.distance);
247 result += (cutoff == 0.) ?
248 0. :
249 cutoff *
250 function_prefactor(
251 alpha,
252 function_eta(r_ij))
253 * function_smoother(
254 1.,
255 params[lambda],
256 r_ij.distance);
257// cutoff * function_prefactor(
258// alpha,
259// function_eta(r_ij))
260// * function_smoother(
261// 1.,
262// params[mu],
263// r_ij.distance);
264 break;
265 }
266 case B:
267 {
268 const double cutoff = function_cutoff(r_ij.distance);
269 result += (cutoff == 0.) ?
270 0. :
271 cutoff * function_prefactor(
272 params[beta],
273 function_zeta(r_ij))
274 * function_smoother(
275 -1.,
276 params[mu],
277 r_ij.distance);
278// cutoff * function_prefactor(
279// beta,
280// function_zeta(r_ij))
281// * function_smoother(
282// -params[B],
283// params[mu],
284// r_ij.distance)/params[B];
285 break;
286 }
287 case lambda:
288 {
289 const double cutoff = function_cutoff(r_ij.distance);
290 result += (cutoff == 0.) ?
291 0. :
292 -r_ij.distance * cutoff *
293 function_prefactor(
294 alpha,
295 function_eta(r_ij))
296 * function_smoother(
297 params[A],
298 params[lambda],
299 r_ij.distance);
300 break;
301 }
302 case mu:
303 {
304 const double cutoff = function_cutoff(r_ij.distance);
305 result += (cutoff == 0.) ?
306 0. :
307 -r_ij.distance * cutoff *(
308 function_prefactor(
309 params[beta],
310 function_zeta(r_ij))
311 * function_smoother(
312 -params[B],
313 params[mu],
314 r_ij.distance)
315 );
316 break;
317 }
318// case lambda3:
319// {
320// result += 0.;
321// break;
322// }
323// case alpha:
324// {
325// const double temp =
326// pow(alpha*function_eta(r_ij), params[n]);
327// const double cutoff = function_cutoff(r_ij.distance);
328// result += (cutoff == 0.) || (alpha == 0. )?
329// 0. :
330// function_smoother(
331// params[A],
332// params[lambda],
333// r_ij.distance)
334// * (-.5) * alpha * (temp/alpha)
335// / (1. + temp)
336// ;
337// break;
338// }
339// case chi:
340// {
341// result += 0.;
342// break;
343// }
344// case omega:
345// {
346// result += 0.;
347// break;
348// }
349 case beta:
350 {
351 const double temp =
352 pow(params[beta]*function_zeta(r_ij), params[n]);
353 const double cutoff = function_cutoff(r_ij.distance);
354 result += (cutoff == 0.) || (params[beta] == 0. )?
355 0. : cutoff *
356 function_smoother(
357 -params[B],
358 params[mu],
359 r_ij.distance)
360 * (-.5)
361 * function_prefactor(
362 params[beta],
363 function_zeta(r_ij))
364 * (temp/params[beta])
365 / (1. + temp)
366 ;
367 break;
368 }
369 case n:
370 {
371 const double zeta = function_zeta(r_ij);
372 const double temp = pow( params[beta]*zeta , params[n]);
373 const double cutoff = function_cutoff(r_ij.distance);
374 const double tempres = ((cutoff == 0.) || (zeta == 0.)) ? // zeta must be caught if zero due to log
375 0. : .5 * cutoff *
376 function_smoother(
377 -params[B],
378 params[mu],
379 r_ij.distance)
380 * function_prefactor(
381 params[beta],
382 function_zeta(r_ij))
383 * ( log(1.+temp)/(params[n]*params[n]) - temp
384 * (log(function_zeta(r_ij)) + log(params[beta]))
385 /(params[n]*(1.+temp)))
386 ;
387// if (tempres != tempres)
388// LOG(2, "DEBUG: tempres is NaN.");
389// LOG(2, "DEBUG: Adding " << tempres << " for p.d. w.r.t n, temp=" << temp << ", cutoff=" << cutoff);
390 result += tempres;
391 break;
392 }
393 case c:
394 {
395 const double zeta = function_zeta(r_ij);
396 if (zeta == 0.)
397 break;
398 const double temp =
399 pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
400 const double cutoff = function_cutoff(r_ij.distance);
401 const double tempres = (cutoff == 0.) ?
402 0. : cutoff *
403 function_smoother(
404 -params[B],
405 params[mu],
406 r_ij.distance)
407 * function_prefactor(
408 params[beta],
409 zeta)
410 * (-1.) * temp / (1.+temp*zeta);
411 double factor = function_derivative_c(r_ij);
412 result += tempres*factor;
413 if (result != result)
414 ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
415 break;
416 }
417 case d:
418 {
419 const double zeta = function_zeta(r_ij);
420 const double temp =
421 pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
422 const double cutoff = function_cutoff(r_ij.distance);
423 const double tempres = (cutoff == 0.) ?
424 0. : cutoff *
425 function_smoother(
426 -params[B],
427 params[mu],
428 r_ij.distance)
429 * function_prefactor(
430 params[beta],
431 zeta)
432 * (-1.) * temp / (1.+temp*zeta);
433 double factor = function_derivative_d(r_ij);
434 result += tempres*factor;
435 if (result != result)
436 ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
437 break;
438 }
439 case h:
440 {
441 const double zeta = function_zeta(r_ij);
442 const double temp =
443 pow(zeta, params[n]-1.) * pow(params[beta],params[n]);
444 const double cutoff = function_cutoff(r_ij.distance);
445 const double tempres = (cutoff == 0.) ?
446 0. : cutoff *
447 function_smoother(
448 -params[B],
449 params[mu],
450 r_ij.distance)
451 * function_prefactor(
452 params[beta],
453 zeta)
454 * (-1.) * temp / (1.+temp*zeta);
455 double factor = function_derivative_h(r_ij);
456 result += tempres*factor;
457 if (result != result)
458 ELOG(1, "result is NaN, zeta=" << zeta << ", temp=" << temp << ", cutoff=" << cutoff << ", tempres=" << tempres << ", factor=" << factor);
459 break;
460 }
461 case offset:
462 result += 1.;
463 break;
464 default:
465 break;
466 }
467 if (result != result)
468 ELOG(1, "result is NaN.");
469 }
470 return results_t(1,-result);
471}
472
473ManyBodyPotential_Tersoff::result_t
474ManyBodyPotential_Tersoff::function_derivative_c(
475 const argument_t &r_ij
476 ) const
477{
478 double result = 0.;
479 std::vector<arguments_t> triples = triplefunction(r_ij, S);
480 for (std::vector<arguments_t>::const_iterator iter = triples.begin();
481 iter != triples.end(); ++iter) {
482 ASSERT( iter->size() == 2,
483 "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
484 const argument_t &r_ik = (*iter)[0];
485 const argument_t &r_jk = (*iter)[1];
486 const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
487 const double cutoff = function_cutoff(r_ik.distance);
488 result += (cutoff == 0.) ?
489 0. : cutoff * omega * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
490 params[c]/Helpers::pow(params[d],2)
491 - params[c] / ( Helpers::pow(params[d],2) + Helpers::pow(tempangle,2) )
492 );
493 }
494 return result;
495}
496
497ManyBodyPotential_Tersoff::result_t
498ManyBodyPotential_Tersoff::function_derivative_d(
499 const argument_t &r_ij
500 ) const
501{
502 double result = 0.;
503 std::vector<arguments_t> triples = triplefunction(r_ij, S);
504 for (std::vector<arguments_t>::const_iterator iter = triples.begin();
505 iter != triples.end(); ++iter) {
506 ASSERT( iter->size() == 2,
507 "ManyBodyPotential_Tersoff::function_derivative_d() - the triples result must contain exactly two distances.");
508 const argument_t &r_ik = (*iter)[0];
509 const argument_t &r_jk = (*iter)[1];
510 const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
511 const double cutoff = function_cutoff(r_ik.distance);
512 result += (cutoff == 0.) ?
513 0. : cutoff * omega * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
514 - Helpers::pow(params[c],2)/Helpers::pow(params[d],3)
515 + Helpers::pow(params[c],2) * params[d]
516 / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(tempangle,2),2)
517 );
518 }
519 return result;
520}
521
522ManyBodyPotential_Tersoff::result_t
523ManyBodyPotential_Tersoff::function_derivative_h(
524 const argument_t &r_ij
525 ) const
526{
527 double result = 0.;
528 std::vector<arguments_t> triples = triplefunction(r_ij, S);
529 for (std::vector<arguments_t>::const_iterator iter = triples.begin();
530 iter != triples.end(); ++iter) {
531 ASSERT( iter->size() == 2,
532 "ManyBodyPotential_Tersoff::function_derivative_h() - the triples result must contain exactly two distances.");
533 const argument_t &r_ik = (*iter)[0];
534 const argument_t &r_jk = (*iter)[1];
535 const double tempangle = params[h] - function_theta(r_ij.distance, r_ik.distance, r_jk.distance);
536 const double cutoff = function_cutoff(r_ik.distance);
537 result += (cutoff == 0.) ?
538 0. : cutoff * omega * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3)) * (
539 ( Helpers::pow(params[c],2)*tempangle )
540 / Helpers::pow(Helpers::pow(params[d],2) + Helpers::pow(tempangle,2),2)
541 );
542 }
543 return result;
544}
545
546ManyBodyPotential_Tersoff::result_t
547ManyBodyPotential_Tersoff::function_cutoff(
548 const double &distance
549 ) const
550{
551// Info info(__func__);
552 double result = 0.;
553 if (distance < R)
554 result = 1.;
555 else if (distance > S)
556 result = 0.;
557 else {
558 result = (0.5 + 0.5 * cos( M_PI * (distance - R)/(S-R)));
559 }
560// LOG(2, "DEBUG: function_cutoff(" << distance << ") = " << result);
561 return result;
562}
563
564ManyBodyPotential_Tersoff::result_t
565ManyBodyPotential_Tersoff::function_prefactor(
566 const double &alpha,
567 const double &eta
568 ) const
569{
570// Info info(__func__);
571 const double result = chi * pow(
572 (1. + pow(alpha * eta, params[n])),
573 -1./(2.*params[n]));
574// LOG(2, "DEBUG: function_prefactor(" << alpha << "," << eta << ") = " << result);
575 return result;
576}
577
578ManyBodyPotential_Tersoff::result_t
579ManyBodyPotential_Tersoff::function_smoother(
580 const double &prefactor,
581 const double &lambda,
582 const double &distance
583 ) const
584{
585// Info info(__func__);
586 const double result = prefactor * exp(-lambda * distance);
587// LOG(2, "DEBUG: function_smoother(" << prefactor << "," << lambda << "," << distance << ") = " << result);
588 return result;
589}
590
591ManyBodyPotential_Tersoff::result_t
592ManyBodyPotential_Tersoff::function_eta(
593 const argument_t &r_ij
594 ) const
595{
596// Info info(__func__);
597 result_t result = 0.;
598
599 // get all triples within the cutoff
600 std::vector<arguments_t> triples = triplefunction(r_ij, S);
601 for (std::vector<arguments_t>::const_iterator iter = triples.begin();
602 iter != triples.end(); ++iter) {
603 ASSERT( iter->size() == 2,
604 "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain of exactly two distances.");
605 const argument_t &r_ik = (*iter)[0];
606 result += function_cutoff(r_ik.distance)
607 * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
608 }
609
610// LOG(2, "DEBUG: function_eta(" << r_ij.distance << ") = " << result);
611 return result;
612}
613
614ManyBodyPotential_Tersoff::result_t
615ManyBodyPotential_Tersoff::function_zeta(
616 const argument_t &r_ij
617 ) const
618{
619// Info info(__func__);
620 result_t result = 0.;
621
622 // get all triples within the cutoff
623 std::vector<arguments_t> triples = triplefunction(r_ij, S);
624 for (std::vector<arguments_t>::const_iterator iter = triples.begin();
625 iter != triples.end(); ++iter) {
626 ASSERT( iter->size() == 2,
627 "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain exactly two distances.");
628 const argument_t &r_ik = (*iter)[0];
629 const argument_t &r_jk = (*iter)[1];
630 result +=
631 function_cutoff(r_ik.distance)
632 * omega
633 * function_angle(r_ij.distance, r_ik.distance, r_jk.distance)
634 * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
635 }
636
637// LOG(2, "DEBUG: function_zeta(" << r_ij.distance << ") = " << result);
638 return result;
639}
640
641ManyBodyPotential_Tersoff::result_t
642ManyBodyPotential_Tersoff::function_theta(
643 const double &r_ij,
644 const double &r_ik,
645 const double &r_jk
646 ) const
647{
648 const double angle = Helpers::pow(r_ij,2) + Helpers::pow(r_ik,2) - Helpers::pow(r_jk,2);
649 const double divisor = 2.* r_ij * r_ik;
650 if (divisor != 0.) {
651 LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
652 return angle/divisor;
653 } else
654 return 0.;
655}
656
657ManyBodyPotential_Tersoff::result_t
658ManyBodyPotential_Tersoff::function_angle(
659 const double &r_ij,
660 const double &r_ik,
661 const double &r_jk
662 ) const
663{
664// Info info(__func__);
665 const double result =
666 1.
667 + (Helpers::pow(params[c]/params[d], 2))
668 - Helpers::pow(params[c], 2)/(Helpers::pow(params[d], 2) +
669 Helpers::pow(params[h] - function_theta(r_ij, r_ik, r_jk),2));
670
671// LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
672 return result;
673}
674
Note: See TracBrowser for help on using the repository browser.