// // edip_multi.hpp // // LGPL Version 2.1 HEADER START // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public // // License as published by the Free Software Foundation; either // version 2.1 of the License, or (at your option) any later version. // // This library is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU // Lesser General Public License for more details. // // You should have received a copy of the GNU Lesser General Public // License along with this library; if not, write to the Free Software // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, // MA 02110-1301 USA // // LGPL Version 2.1 HEADER END // // // Copyright (c) 2021, Regents of the University of Minnesota. // All rights reserved. // // Contributors: // Yaser Afshar // // Brief: This file is adapted from the LAMMPS software package // `lammps/src/MANYBODY/pair_edip_multi.cpp` // and it is rewritten and updated for the KIM-API by // Yaser Afshar // #ifndef EDIP_MULTI_HPP #define EDIP_MULTI_HPP #include #include #include #include #include #include "edip_param.hpp" #include "helper.hpp" /*! * \brief The EDIPMulti class * * This style compute a 3-body EDIP potential for modeling materials where it * can have advantages over other models such as the Stillinger-Weber or * Tersoff potentials. This edip style supports multi-element. * * Reference: * Justo, Jo\~ao F. et al., "Interatomic potential for silicon defects and * disordered phases," Phys. Rev. B, 58(5), 2539--2550, 1998. * */ class EDIPMulti { public: /*! * \brief Process and parse the edip/c parameter file * * \param parameter_file_pointer FILE pointer to the opened file * \param max_line_size maximum line size * \param element_name name of the unique element * \return int 0|false if everything goes well and 1|true if it fails */ int ProcessParameterFile(std::FILE *const parameter_file_pointer, int const max_line_size, std::vector const &element_name); /*! * \brief Convert units of the parameters * * \param convert_length_factor length unit conversion factor * \param convert_energy_factor energy unit conversion factor */ void ConvertUnit(double const convert_length_factor, double const convert_energy_factor); /*! * \brief Complete the set up * * \param max_cutoff */ void CompleteSetup(double *max_cutoff); /*! * \brief Compute pair Vij, partial derivatives dVij(r,Z)/dr and dVij(r,Z)/dZ * (Eq.4) * * \param r Distance between pair i-j * \param z Effective coordination number * \param ijparam Index * \param fdr Partial derivatives dVij(r,Z)/dr * \param fZ dVij(r,Z)/dZ * \return Two-body term includes repulsive and attractive interactions */ inline double EdipPair(double const r, double const z, std::size_t const ijparam, double &fdr, double &fZ) const; /*! * \brief function fc(r) in calculating coordination Z and derivative fc'(r) * * \param r Distance between pair i-j * \param ijparam Index * \param fdr fc'(r) * \return fc(r) */ inline double EdipFc(double const r, std::size_t const ijparam, double &fdr) const; /*! * \brief function tau(Z) and its derivative tau'(Z) * * \param z Effective coordination number * \param ijparam Index * \param fdZ tau'(Z) * \return tau(Z) */ inline double EdipTau(double const z, std::size_t const ijparam, double &fdZ) const; /*! * \brief function h(l,Z) and its partial derivatives dh(l,Z)/dl and * dh(l,Z)/dZ * * \param l cosThetaijk * \param z Effective coordination number * \param ijparam Index * \param fdl partial derivatives dh(l,Z)/dl * \param fdZ partial derivatives dh(l,Z)/dZ * \return function h(l,Z) */ inline double EdipH(double const l, double const z, std::size_t const ijparam, double &fdl, double &fdZ) const; /*! * \brief cut-off function for Vij and its derivative fcut2'(r) * * \param r Distance * \param ijparam Index * \param fdr derivative fcut2'(r) * \return cut-off function */ inline double EdipFcut2(double const r, std::size_t const ijparam, double &fdr) const; /*! * \brief cut-off function for Vijk and its derivative fcut3'(r) * * \param r Distance * \param ijparam Index * \param fdr derivative fcut3'(r) * \return cut-off function */ inline double EdipFcut3(double const r, std::size_t const ijparam, double &fdr) const; public: /*! Parameter set for an I-J-K interaction */ std::vector params_; /*! Length unit conversion factor */ double convert_length_factor_{1.0}; /*! Energy unit conversion factor */ double convert_energy_factor_{1.0}; // max cutoff for all elements double cutmax_{0.0}; /*! Helper array, keep three-body idex */ Array3D elem3param_; /*! Temporary array */ std::vector pre_force_coord_data_; }; inline double EDIPMulti::EdipPair(double const r, double const z, std::size_t const ijparam, double &fdr, double &fZ) const { double const A = params_[ijparam].A; double const B = params_[ijparam].B; double const rho = params_[ijparam].rho; double const beta = params_[ijparam].beta; double const v1 = std::pow(B / r, rho); double const v2 = std::exp(-beta * z * z); double v4; double const v3 = EdipFcut2(r, ijparam, v4); fdr = A * (v1 - v2) * v4 + A * (-rho * v1 / r) * v3; fZ = A * (2 * beta * z * v2) * v3; double const eng = A * (v1 - v2) * v3; return eng; } inline double EDIPMulti::EdipFc(double const r, std::size_t const ijparam, double &fdr) const { double const cutoffC = params_[ijparam].cutoffC; double const r_cutoffC = r - cutoffC; if (r_cutoffC < 1E-6) { fdr = 0.0; return 1.0; } double const cutoffA = params_[ijparam].cutoffA; double const r_cutoffA = r - cutoffA; if (r_cutoffA > -1E-6) { fdr = 0.0; return 0.0; } double const cutoffA_cutoffC = cutoffA - cutoffC; double const x = cutoffA_cutoffC / r_cutoffC; double const v1 = x * x * x; double const v2 = 1.0 / (1.0 - v1); double const alpha = params_[ijparam].alpha; double const f = std::exp(alpha * v2); fdr = (3.0 * x * v1 / cutoffA_cutoffC) * (-alpha * v2 * v2) * f; return f; } // function tau(Z) and its derivative tau'(Z) inline double EDIPMulti::EdipTau(double const z, std::size_t const ijparam, double &fdZ) const { double const u1 = params_[ijparam].u1; double const u2 = params_[ijparam].u2; double const u3 = params_[ijparam].u3; double const u4 = params_[ijparam].u4; double const v1 = std::exp(-u4 * z); double const v2 = std::exp(-2.0 * u4 * z); fdZ = -u2 * u3 * u4 * v1 + 2.0 * u2 * u4 * v2; double const f = u1 + u2 * u3 * v1 - u2 * v2; return f; } // function h(l,Z) and its partial derivatives dh(l,Z)/dl and dh(l,Z)/dZ inline double EDIPMulti::EdipH(double const l, double const z, std::size_t const ijparam, double &fdl, double &fdZ) const { double const lambda = params_[ijparam].lambda; double const eta = params_[ijparam].eta; double const Q0 = params_[ijparam].Q0; double const mu = params_[ijparam].mu; // function Q(Z) double const Q = Q0 * std::exp(-mu * z); // derivative Q'(Z) double const QdZ = -mu * Q; double TaudZ; double const Tau = EdipTau(z, ijparam, TaudZ); double const l_Tau = l + Tau; double const v1 = l_Tau * l_Tau; double const u2 = Q * v1; double const v2 = std::exp(-u2); double const f = lambda * (1 - v2 + eta * u2); // df/du2 double const v3 = lambda * (v2 + eta); // du2/dl double const du2l = Q * 2 * (l + Tau); fdl = v3 * du2l; // du2/dZ double const du2Z = QdZ * v1 + Q * 2 * (l + Tau) * TaudZ; fdZ = v3 * du2Z; return f; } // cut-off function for Vij and its derivative fcut2'(r) inline double EDIPMulti::EdipFcut2(double const r, std::size_t const ijparam, double &fdr) const { double const cutoffA = params_[ijparam].cutoffA; double const r_cutoffA = r - cutoffA; if (r_cutoffA > -1E-6) { fdr = 0.0; return 0.0; } double const sigma = params_[ijparam].sigma; double const v1 = 1.0 / r_cutoffA; double const sv = sigma * v1; double const f = std::exp(sv); fdr = (-sv * v1) * f; return f; } // cut-off function for Vijk and its derivative fcut3'(r) inline double EDIPMulti::EdipFcut3(double const r, std::size_t const ijparam, double &fdr) const { double const cutoffA = params_[ijparam].cutoffA; double const r_cutoffA = r - cutoffA; if (r_cutoffA > -1E-6) { fdr = 0.0; return 0.0; } double const gamma = params_[ijparam].gamma; double const v1 = 1.0 / r_cutoffA; double const gv = gamma * v1; double const f = std::exp(gv); fdr = (-gv * v1) * f; return f; } #endif // EDIP_MULTI