// // edip_c.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 // N. A. Marks & Sam McSweeney (Curtin University) EDIP/C code // and it is rewritten and updated for the KIM-API by // Yaser Afshar // #ifndef EDIP_C_HPP #define EDIP_C_HPP #include #include #include #include "edip_param_c.hpp" #include "fermi.hpp" #include "helper.hpp" #include "zbl.hpp" /*! * \brief Carbon EDIP * * Reference: * Marks, N. A., "Generalizing the environment-dependent interaction potential * for carbon," Phys. Rev.B, 63(3):035401, Dec 2000. */ class EDIPC { 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 * \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); /*! * \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 Grow local arrays if necessary * * \param number_of_particles Number of particles * \param total_number_of_neighbors Total number of neighbors */ void Resize(std::size_t const number_of_particles, std::size_t const total_number_of_neighbors); /*! * \brief Grow local arrays if necessary * * \param number_of_particles Number of particles */ void Resize(std::size_t const number_of_particles); /*! * \brief Complete the set up * * \param max_cutoff */ void CompleteSetup(double *max_cutoff); /*! * \brief Calculate Cutoff functions * * \param number_of_particles Number of particles */ void CutoffFunction(std::size_t const number_of_particles); /*! * \brief Calculate Switching functions * * \param number_of_particles Number of particles */ void SwitchFunction(std::size_t const number_of_particles); /*! * \brief Calculate generalized coordination Z_i for all atoms * * \param i Particle index */ void Coordination(std::size_t const i); /*! * \brief Calculate dihedral contribution to Z_i * * \param i Particle index */ void DihedralContribution(std::size_t const i); /*! * \brief Calculate repulsion contribution to Z_i * * \param i Particle index */ void RepulsionContribution(std::size_t const i); /*! * \brief Calculate linear contribution to Z_i * * \param i Particle index */ void LinearContribution(std::size_t const i); /*! * \brief Calculate the pair energy for a given (i,jj)-pair * * \param i Particle index * \param r_i_jj Distance between particle i and its jjth neighbor * \param u2_arg * \param u2_bond * \param u2_r5 * \param u2_part1 * \param u2_part2 * \param f22 * \return double The pair energy for a given (i,jj)-pair */ double CalcU2(std::size_t const i, double const r_i_jj, double &u2_arg, double &u2_bond, double &u2_r5, double &u2_part1, double &u2_part2, double &f22) const; /*! * \brief Calculate the pair force due to neighbours of i * * \param i Particle index * \param jj Neighbor index * \param kk Neighbor index * \param u2_arg * \param u2_bond * \param u2_r5 * \param u2_part1 * \param u2_part2 * \param f2 The pair force due to neighbours of i * \param dxr */ void CalcF2(std::size_t const i, std::size_t const jj, std::size_t const kk, double const u2_arg, double const u2_bond, double const u2_r5, double const u2_part1, double const u2_part2, VectorOfSizeDIM f2, VectorOfSizeDIM dxr) const; /*! * \brief Calculate pair energy term * * \tparam IsComputeEnergy ComputeEnergy flag * \tparam IsComputeForces ComputeForces flag * \tparam IsComputeParticleEnergy ComputeParticleEnergy flag * \tparam IsComputeVirial ComputeVirial flag * \tparam IsComputeParticleVirial ComputeParticleVirial flag * * \param i Particle index * \param particle_contributing Contribution array * \param energy Global energy value * \param forces Particles forces * \param particle_energy Particles energy * \param virial Global virial term * \param particle_virial Particles virial */ template void Pair(std::size_t const i, int const *const particle_contributing, double *const energy, VectorOfSizeDIM *const forces, double *const particle_energy, VectorOfSizeSix virial, VectorOfSizeSix *const particle_virial); /*! * \brief Calculate pair & zbl energy term * * \tparam IsComputeEnergy ComputeEnergy flag * \tparam IsComputeForces ComputeForces flag * \tparam IsComputeParticleEnergy ComputeParticleEnergy flag * \tparam IsComputeVirial ComputeVirial flag * \tparam IsComputeParticleVirial ComputeParticleVirial flag * * \param i Particle index * \param particle_contributing Contribution array * \param energy Global energy value * \param forces Particles forces * \param particle_energy Particles energy * \param virial Global virial term * \param particle_virial Particles virial */ template void PairZbl(std::size_t const i, int const *const particle_contributing, double *const energy, VectorOfSizeDIM *const forces, double *const particle_energy, VectorOfSizeSix virial, VectorOfSizeSix *const particle_virial); /*! * \brief Calculate tau and dtau * * \f$ \tau \f$ describes the variation in ideal bonding angle, favoring * angles of 180°, 120°, 109.5°, and 90° for coordi- nations 2, 3, 4, and 6, * respectively. * * \f$ \tau(Z) = 1 - \frac{Z}{12} \tanh \left[t_1 (Z - t_2)\right] \f$ * * \param i Particle index * \param dtau Derivative of tau * \return double tau */ double TauFunction(std::size_t const i, double &dtau) const; /*! * \brief Calculate the triple energy for a given (i,jj)-pair * * \param i Particle index * \param jj Neighbor index * \param kk Neighbor index * \param r_i_jj Distance between particle i and its jj neighbor * \param r_i_kk Distance between particle i and its kk neighbor * \param tau * \param dtau * \param dcosj * \param dcosk * \param zexpgam * \param dtheta * \param u3_arg11 * \param u3_arg22 * \param u3_argz * \param f33 * \return double */ double CalcU3(std::size_t const i, std::size_t const jj, std::size_t const kk, double const r_i_jj, double const r_i_kk, double const tau, double const dtau, VectorOfSizeDIM dcosj, VectorOfSizeDIM dcosk, double &zexpgam, double &dtheta, double &u3_arg11, double &u3_arg22, double &u3_argz, double &f33) const; /*! * \brief Calculate the triple force due to neighbours of i * * \param i Particle index * \param jj Neighbor index * \param kk Neighbor index * \param ll Neighbor index * \param dtau Derivative of tau * \param dcosj * \param dcosk * \param zexpgam * \param dtheta * \param u3_arg11 * \param u3_arg22 * \param u3_argz * \param f3 */ void CalcF3(std::size_t const i, std::size_t const jj, std::size_t const kk, std::size_t const ll, double const dtau, VectorOfSizeDIM dcosj, VectorOfSizeDIM dcosk, double const zexpgam, double const dtheta, double const u3_arg11, double const u3_arg22, double const u3_argz, VectorOfSizeDIM f3) const; /*! * \brief Calculate triple energy term * * \tparam IsComputeEnergy ComputeEnergy flag * \tparam IsComputeForces ComputeForces flag * \tparam IsComputeParticleEnergy ComputeParticleEnergy flag * \tparam IsComputeVirial ComputeVirial flag * \tparam IsComputeParticleVirial ComputeParticleVirial flag * * \param i Particle index * \param particle_contributing Contribution array * \param energy Global energy value * \param forces Particles forces * \param particle_energy Particles energy * \param virial Global virial term * \param particle_virial Particles virial */ template void Triple(std::size_t const i, int const *const particle_contributing, double *const energy, VectorOfSizeDIM *const forces, double *const particle_energy, VectorOfSizeSix virial, VectorOfSizeSix *const particle_virial); public: /*! Parameters for Carbon EDIP */ EDIPCParam params_; /*! * \brief zbl flag * * Flag indicating whether to use C interactions in combination with * the standard Ziegler-Biersack-Littmark (ZBL) potential. */ bool lzbl_{false}; /*! length unit conversion factor */ double convert_length_factor_{1.0}; /*! energy unit conversion factor */ double convert_energy_factor_{1.0}; /*! Number of neighbors for each atom */ std::vector number_of_neighbors_; /*! Particle index for each neighbor of a particle */ std::vector neighbors_to_particles_; /*! Offset (cumulative distance) of the first neighbor of a particle */ std::vector neighbors_offset_; /*! Distance between two atoms */ std::vector r_; /* Normal vector between two atoms */ std::vector N_; /*! Spherical coordination */ std::vector z_; /*! Derivative of z for every pair of atoms */ std::vector dz_; /*! p(r) for every pair of atoms */ std::vector p_; /*! dp/dr for every pair of atoms */ std::vector dp_; // Switching functions std::vector pi_; std::vector dpi_; std::vector pi2_; std::vector dpi2_; std::vector pi3_; std::vector dpi3_; /*! Generalized coordination */ std::vector Z_; std::vector Dzdx_; Array2D Dzdxx_; /*! Early detection of zero-force pairs */ std::vector finiteforce_; Array2D finiteforce2_; /*! ZBL object */ std::unique_ptr zbl_{nullptr}; /*! FERMI object */ std::unique_ptr fermi_{nullptr}; }; #endif // EDIP_C