diff --git a/src/USER-MEAMC/fm_exp.cpp b/src/USER-MEAMC/fm_exp.cpp deleted file mode 100644 index 253af868f7..0000000000 --- a/src/USER-MEAMC/fm_exp.cpp +++ /dev/null @@ -1,135 +0,0 @@ -extern "C" { -/* - Copyright (c) 2012,2013 Axel Kohlmeyer - All rights reserved. - - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions - are met: - - * Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright - notice, this list of conditions and the following disclaimer in the - documentation and/or other materials provided with the distribution. - * Neither the name of the nor the - names of its contributors may be used to endorse or promote products - derived from this software without specific prior written permission. - -THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -ARE DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY -DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND -ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF -THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -*/ - -/* faster versions of 2**x, e**x, and 10**x in single and double precision. - * - * Based on the Cephes math library 2.8 - */ - -#include -#include - -/* internal definitions for the fastermath library */ - -/* IEEE 754 double precision floating point data manipulation */ -typedef union -{ - double f; - uint64_t u; - struct {int32_t i0,i1;}; -} udi_t; -#define FM_DOUBLE_BIAS 1023 -#define FM_DOUBLE_EMASK 2146435072 -#define FM_DOUBLE_MBITS 20 -#define FM_DOUBLE_MMASK 1048575 -#define FM_DOUBLE_EZERO 1072693248 - -/* generate 2**num in floating point by bitshifting */ -#define FM_DOUBLE_INIT_EXP(var,num) \ - var.i0 = 0; \ - var.i1 = (((int) num) + FM_DOUBLE_BIAS) << 20 - -/* double precision constants */ -#define FM_DOUBLE_LOG2OFE 1.4426950408889634074 -#define FM_DOUBLE_LOGEOF2 6.9314718055994530942e-1 -#define FM_DOUBLE_LOG2OF10 3.32192809488736234789 -#define FM_DOUBLE_LOG10OF2 3.0102999566398119521e-1 -#define FM_DOUBLE_LOG10OFE 4.3429448190325182765e-1 -#define FM_DOUBLE_SQRT2 1.41421356237309504880 -#define FM_DOUBLE_SQRTH 0.70710678118654752440 - -/* optimizer friendly implementation of exp2(x). - * - * strategy: - * - * split argument into an integer part and a fraction: - * ipart = floor(x+0.5); - * fpart = x - ipart; - * - * compute exp2(ipart) from setting the ieee754 exponent - * compute exp2(fpart) using a pade' approximation for x in [-0.5;0.5[ - * - * the result becomes: exp2(x) = exp2(ipart) * exp2(fpart) - */ - -static const double fm_exp2_q[] = { -/* 1.00000000000000000000e0, */ - 2.33184211722314911771e2, - 4.36821166879210612817e3 -}; -static const double fm_exp2_p[] = { - 2.30933477057345225087e-2, - 2.02020656693165307700e1, - 1.51390680115615096133e3 -}; - -static double fm_exp2(double x) -{ - double ipart, fpart, px, qx; - udi_t epart; - - ipart = floor(x+0.5); - fpart = x - ipart; - FM_DOUBLE_INIT_EXP(epart,ipart); - - x = fpart*fpart; - - px = fm_exp2_p[0]; - px = px*x + fm_exp2_p[1]; - qx = x + fm_exp2_q[0]; - px = px*x + fm_exp2_p[2]; - qx = qx*x + fm_exp2_q[1]; - - px = px * fpart; - - x = 1.0 + 2.0*(px/(qx-px)); - return epart.f*x; -} - -double fm_exp(double x) -{ -#if defined(__BYTE_ORDER__) -#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ - return fm_exp2(FM_DOUBLE_LOG2OFE * x); -#endif -#endif - return exp(x); -} - -/* - * Local Variables: - * mode: c - * compile-command: "make -C .." - * c-basic-offset: 4 - * fill-column: 76 - * indent-tabs-mode: nil - * End: - */ -} diff --git a/src/USER-MEAMC/meam.h b/src/USER-MEAMC/meam.h index 8f652c8d1e..0590d2fe05 100644 --- a/src/USER-MEAMC/meam.h +++ b/src/USER-MEAMC/meam.h @@ -2,31 +2,25 @@ #define LMP_MEAM_H #include +#include "memory.h" #define maxelt 5 -extern "C" { - double fm_exp(double); -} - namespace LAMMPS_NS { typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t; -struct allocatable_double_2d { - allocatable_double_2d() : dim1(0),dim2(0),ptr(0) {}; - int dim1, dim2; - double* ptr; -}; - class MEAM { public: - MEAM(){}; + MEAM(Memory *mem) : + memory(mem) {}; ~MEAM() { - meam_cleanup_(); + meam_cleanup(); } private: + Memory *&memory; + // cutforce = force cutoff // cutforcesq = force cutoff squared @@ -89,10 +83,10 @@ class MEAM { int eltind[maxelt + 1][maxelt + 1]; int neltypes; - allocatable_double_2d phir; // [:,:] + double **phir; - allocatable_double_2d phirar, phirar1, phirar2, phirar3, phirar4, phirar5, - phirar6; // [:,:] + double **phirar, **phirar1, **phirar2, **phirar3, **phirar4, **phirar5, + **phirar6; double attrac_meam[maxelt + 1][maxelt + 1], repuls_meam[maxelt + 1][maxelt + 1]; @@ -108,13 +102,10 @@ class MEAM { int nr, nrar; double dr, rdrar; - public: + protected: void meam_checkindex(int, int, int, int*, int*); - void meam_setup_param_(int*, double*, int*, int*, int*); - void meam_dens_final_(int*, int*, int*, int*, int*, double*, double*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); void G_gam(double, int, double, double*, int*); void dG_gam(double, int, double, double*, double*); - void meam_dens_init_(int*, int*, int*, int*, int*, double*, int*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); void getscreen(int, int, double*, double*, double*, double*, int, int*, int, int*, int, int*, int*); void screen(int, int, int, double*, double, double*, int, int*, int, int*, int*); void calc_rho1(int, int, int, int*, int*, double*, int, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*); @@ -124,9 +115,6 @@ class MEAM { void dCfunc(double, double, double, double*); void dCfunc2(double, double, double, double*, double*); - void meam_force_(int*, int*, int*, int*, int*, int*, double*, double*, int*, int*, int*, double*, int*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); - void meam_cleanup_(); - void meam_setup_done_(double*); void alloyparams(); void compute_pair_meam(); double phi_meam(double, int, int); @@ -141,13 +129,17 @@ class MEAM { double erose(double, double, double, double, double, double, int); void interpolate_meam(int); double compute_phi(double, int, int); - void meam_setup_global_(int*, int*, double*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); - + public: + void meam_setup_global(int*, int*, double*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); + void meam_setup_param(int*, double*, int*, int*, int*); + void meam_setup_done(double*); + void meam_dens_init(int*, int*, int*, int*, int*, double*, int*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); + void meam_dens_final(int*, int*, int*, int*, int*, double*, double*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); + void meam_force(int*, int*, int*, int*, int*, int*, double*, double*, int*, int*, int*, double*, int*, int*, int*, int*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*); + void meam_cleanup(); }; // Functions we need for compat -#define MIN(A,B) ((A) < (B) ? (A) : (B)) -#define MAX(A,B) ((A) > (B) ? (A) : (B)) #define iszero(f) (fabs(f) < 1e-20) @@ -171,7 +163,7 @@ Fortran Array Semantics in C. - Multi-Dimensional MUST be declared in reverse, so that the order when accessing is the same as in Fortran - arrays that are passed externally via the meam_* functions must use the arr*v() functions below (or be used with 0-based indexing) - - allocatable arrays (only global phir*) are actually a struct, and must be accessed with arr2() + - allocatable arrays must be accessed with arr2() */ // we receive a pointer to the first element, and F dimensions is ptr(a,b,c) @@ -194,19 +186,9 @@ Fortran Array Semantics in C. ptr[(i - 1) + (j - 1) * (DIM1__##ptr) + \ (k - 1) * (DIM1__##ptr) * (DIM2__##ptr)] -// allocatable arrays via special type -#define allocate_2d(arr, cols, rows) \ - arr.dim1 = cols; \ - arr.dim2 = rows; \ - arr.ptr = (double*)calloc((size_t)(cols) * (size_t)(rows), sizeof(double)); -#define allocated(a) (a.ptr != NULL) -#define meam_deallocate(a) \ - ({ \ - free(a.ptr); \ - a.ptr = NULL; \ - }) +// allocatable arrays // access data with same index as used in fortran (1-based) -#define arr2(arr, i, j) arr.ptr[(arr.dim1) * (j - 1) + (i - 1)] +#define arr2(arr, i, j) arr[j-1][i-1] }; diff --git a/src/USER-MEAMC/meam_cleanup.cpp b/src/USER-MEAMC/meam_cleanup.cpp index 1b1dc45958..a0184451e8 100644 --- a/src/USER-MEAMC/meam_cleanup.cpp +++ b/src/USER-MEAMC/meam_cleanup.cpp @@ -3,15 +3,15 @@ using namespace LAMMPS_NS; void -MEAM::meam_cleanup_(void) +MEAM::meam_cleanup(void) { - meam_deallocate(this->phirar6); - meam_deallocate(this->phirar5); - meam_deallocate(this->phirar4); - meam_deallocate(this->phirar3); - meam_deallocate(this->phirar2); - meam_deallocate(this->phirar1); - meam_deallocate(this->phirar); - meam_deallocate(this->phir); + memory->destroy(this->phirar6); + memory->destroy(this->phirar5); + memory->destroy(this->phirar4); + memory->destroy(this->phirar3); + memory->destroy(this->phirar2); + memory->destroy(this->phirar1); + memory->destroy(this->phirar); + memory->destroy(this->phir); } \ No newline at end of file diff --git a/src/USER-MEAMC/meam_dens_final.cpp b/src/USER-MEAMC/meam_dens_final.cpp index 778787a062..0f27c795c3 100644 --- a/src/USER-MEAMC/meam_dens_final.cpp +++ b/src/USER-MEAMC/meam_dens_final.cpp @@ -1,5 +1,6 @@ #include "meam.h" #include +#include "math_special.h" using namespace LAMMPS_NS; // Extern "C" declaration has the form: @@ -20,7 +21,7 @@ using namespace LAMMPS_NS; // void -MEAM::meam_dens_final_(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, +MEAM::meam_dens_final(int* nlocal, int* nmax, int* eflag_either, int* eflag_global, int* eflag_atom, double* eng_vdwl, double* eatom, int* ntype, int* type, int* fmap, double* Arho1, double* Arho2, double* Arho2b, double* Arho3, double* Arho3b, double* t_ave, @@ -221,7 +222,7 @@ MEAM::G_gam(double Gamma, int ibar, double gsmooth_factor, double* G, int* error *G = sqrt(1.0 + Gamma); } } else if (ibar == 1) { - *G = fm_exp(Gamma / 2.0); + *G = MathSpecial::fm_exp(Gamma / 2.0); } else if (ibar == 3) { *G = 2.0 / (1.0 + exp(-Gamma)); } else if (ibar == -5) { @@ -242,9 +243,9 @@ MEAM::dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* d { // Compute G(Gamma) and dG(gamma) based on selection flag ibar: // 0 => G = sqrt(1+Gamma) - // 1 => G = fm_exp(Gamma/2) + // 1 => G = MathSpecial::fm_exp(Gamma/2) // 2 => not implemented - // 3 => G = 2/(1+fm_exp(-Gamma)) + // 3 => G = 2/(1+MathSpecial::fm_exp(-Gamma)) // 4 => G = sqrt(1+Gamma) // -5 => G = +-sqrt(abs(1+Gamma)) double gsmooth_switchpoint; @@ -264,10 +265,10 @@ MEAM::dG_gam(double Gamma, int ibar, double gsmooth_factor, double* G, double* d *dG = 1.0 / (2.0 * *G); } } else if (ibar == 1) { - *G = fm_exp(Gamma / 2.0); + *G = MathSpecial::fm_exp(Gamma / 2.0); *dG = *G / 2.0; } else if (ibar == 3) { - *G = 2.0 / (1.0 + fm_exp(-Gamma)); + *G = 2.0 / (1.0 + MathSpecial::fm_exp(-Gamma)); *dG = *G * (2.0 - *G) / 2; } else if (ibar == -5) { if ((1.0 + Gamma) >= 0) { diff --git a/src/USER-MEAMC/meam_dens_init.cpp b/src/USER-MEAMC/meam_dens_init.cpp index f03bc5a128..98f9f958f7 100644 --- a/src/USER-MEAMC/meam_dens_init.cpp +++ b/src/USER-MEAMC/meam_dens_init.cpp @@ -1,5 +1,6 @@ #include "meam.h" #include +#include "math_special.h" using namespace LAMMPS_NS; // Extern "C" declaration has the form: @@ -21,7 +22,7 @@ using namespace LAMMPS_NS; // void -MEAM::meam_dens_init_(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x, +MEAM::meam_dens_init(int* i, int* nmax, int* ntype, int* type, int* fmap, double* x, int* numneigh, int* firstneigh, int* numneigh_full, int* firstneigh_full, double* scrfcn, double* dscrfcn, double* fcpair, double* rho0, double* arho1, double* arho2, @@ -209,14 +210,14 @@ MEAM::calc_rho1(int i, int nmax, int ntype, int* type, int* fmap, double* x, aj = rij / this->re_meam[eltj][eltj] - 1.0; ro0i = this->rho0_meam[elti]; ro0j = this->rho0_meam[eltj]; - rhoa0j = ro0j * fm_exp(-this->beta0_meam[eltj] * aj) * sij; - rhoa1j = ro0j * fm_exp(-this->beta1_meam[eltj] * aj) * sij; - rhoa2j = ro0j * fm_exp(-this->beta2_meam[eltj] * aj) * sij; - rhoa3j = ro0j * fm_exp(-this->beta3_meam[eltj] * aj) * sij; - rhoa0i = ro0i * fm_exp(-this->beta0_meam[elti] * ai) * sij; - rhoa1i = ro0i * fm_exp(-this->beta1_meam[elti] * ai) * sij; - rhoa2i = ro0i * fm_exp(-this->beta2_meam[elti] * ai) * sij; - rhoa3i = ro0i * fm_exp(-this->beta3_meam[elti] * ai) * sij; + rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj) * sij; + rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj) * sij; + rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj) * sij; + rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj) * sij; + rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai) * sij; + rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai) * sij; + rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai) * sij; + rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai) * sij; if (this->ialloy == 1) { rhoa1j = rhoa1j * this->t1_meam[eltj]; rhoa2j = rhoa2j * this->t2_meam[eltj]; diff --git a/src/USER-MEAMC/meam_force.cpp b/src/USER-MEAMC/meam_force.cpp index f5d4e75799..3b2d0d8b80 100644 --- a/src/USER-MEAMC/meam_force.cpp +++ b/src/USER-MEAMC/meam_force.cpp @@ -1,5 +1,7 @@ #include "meam.h" #include +#include +#include "math_special.h" using namespace LAMMPS_NS; // Extern "C" declaration has the form: @@ -24,7 +26,7 @@ using namespace LAMMPS_NS; // void -MEAM::meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, +MEAM::meam_force(int* iptr, int* nmax, int* eflag_either, int* eflag_global, int* eflag_atom, int* vflag_atom, double* eng_vdwl, double* eatom, int* ntype, int* type, int* fmap, double* x, int* numneigh, int* firstneigh, int* numneigh_full, int* firstneigh_full, @@ -113,9 +115,9 @@ MEAM::meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, ind = this->eltind[elti][eltj]; pp = rij * this->rdrar + 1.0; kk = (int)pp; - kk = MIN(kk, this->nrar - 1); + kk = std::min(kk, this->nrar - 1); pp = pp - kk; - pp = MIN(pp, 1.0); + pp = std::min(pp, 1.0); phi = ((arr2(this->phirar3, kk, ind) * pp + arr2(this->phirar2, kk, ind)) * pp + @@ -144,26 +146,26 @@ MEAM::meam_force_(int* iptr, int* nmax, int* eflag_either, int* eflag_global, invrei = 1.0 / this->re_meam[elti][elti]; ai = rij * invrei - 1.0; ro0i = this->rho0_meam[elti]; - rhoa0i = ro0i * fm_exp(-this->beta0_meam[elti] * ai); + rhoa0i = ro0i * MathSpecial::fm_exp(-this->beta0_meam[elti] * ai); drhoa0i = -this->beta0_meam[elti] * invrei * rhoa0i; - rhoa1i = ro0i * fm_exp(-this->beta1_meam[elti] * ai); + rhoa1i = ro0i * MathSpecial::fm_exp(-this->beta1_meam[elti] * ai); drhoa1i = -this->beta1_meam[elti] * invrei * rhoa1i; - rhoa2i = ro0i * fm_exp(-this->beta2_meam[elti] * ai); + rhoa2i = ro0i * MathSpecial::fm_exp(-this->beta2_meam[elti] * ai); drhoa2i = -this->beta2_meam[elti] * invrei * rhoa2i; - rhoa3i = ro0i * fm_exp(-this->beta3_meam[elti] * ai); + rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai); drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i; if (elti != eltj) { invrej = 1.0 / this->re_meam[eltj][eltj]; aj = rij * invrej - 1.0; ro0j = this->rho0_meam[eltj]; - rhoa0j = ro0j * fm_exp(-this->beta0_meam[eltj] * aj); + rhoa0j = ro0j * MathSpecial::fm_exp(-this->beta0_meam[eltj] * aj); drhoa0j = -this->beta0_meam[eltj] * invrej * rhoa0j; - rhoa1j = ro0j * fm_exp(-this->beta1_meam[eltj] * aj); + rhoa1j = ro0j * MathSpecial::fm_exp(-this->beta1_meam[eltj] * aj); drhoa1j = -this->beta1_meam[eltj] * invrej * rhoa1j; - rhoa2j = ro0j * fm_exp(-this->beta2_meam[eltj] * aj); + rhoa2j = ro0j * MathSpecial::fm_exp(-this->beta2_meam[eltj] * aj); drhoa2j = -this->beta2_meam[eltj] * invrej * rhoa2j; - rhoa3j = ro0j * fm_exp(-this->beta3_meam[eltj] * aj); + rhoa3j = ro0j * MathSpecial::fm_exp(-this->beta3_meam[eltj] * aj); drhoa3j = -this->beta3_meam[eltj] * invrej * rhoa3j; } else { rhoa0j = rhoa0i; diff --git a/src/USER-MEAMC/meam_setup_done.cpp b/src/USER-MEAMC/meam_setup_done.cpp index bea6962de5..5e1ad1f4d2 100644 --- a/src/USER-MEAMC/meam_setup_done.cpp +++ b/src/USER-MEAMC/meam_setup_done.cpp @@ -1,5 +1,7 @@ #include "meam.h" #include +#include +#include "math_special.h" using namespace LAMMPS_NS; // Declaration in pair_meam.h: @@ -12,7 +14,7 @@ using namespace LAMMPS_NS; // void -MEAM::meam_setup_done_(double* cutmax) +MEAM::meam_setup_done(double* cutmax) { int nv2, nv3, m, n, p; @@ -160,7 +162,7 @@ MEAM::alloyparams(void) for (k = 1; k <= this->neltypes; k++) { eb = (this->Cmax_meam[i][j][k] * this->Cmax_meam[i][j][k]) / (4.0 * (this->Cmax_meam[i][j][k] - 1.0)); - this->ebound_meam[i][j] = MAX(this->ebound_meam[i][j], eb); + this->ebound_meam[i][j] = std::max(this->ebound_meam[i][j], eb); } } } @@ -183,43 +185,51 @@ MEAM::compute_pair_meam(void) double C, s111, s112, s221, S11, S22; // check for previously allocated arrays and free them - if (allocated(this->phir)) - meam_deallocate(this->phir); - if (allocated(this->phirar)) - meam_deallocate(this->phirar); - if (allocated(this->phirar1)) - meam_deallocate(this->phirar1); - if (allocated(this->phirar2)) - meam_deallocate(this->phirar2); - if (allocated(this->phirar3)) - meam_deallocate(this->phirar3); - if (allocated(this->phirar4)) - meam_deallocate(this->phirar4); - if (allocated(this->phirar5)) - meam_deallocate(this->phirar5); - if (allocated(this->phirar6)) - meam_deallocate(this->phirar6); + if (this->phir != NULL) + memory->destroy(this->phir); + if (this->phirar != NULL) + memory->destroy(this->phirar); + if (this->phirar1 != NULL) + memory->destroy(this->phirar1); + if (this->phirar2 != NULL) + memory->destroy(this->phirar2); + if (this->phirar3 != NULL) + memory->destroy(this->phirar3); + if (this->phirar4 != NULL) + memory->destroy(this->phirar4); + if (this->phirar5 != NULL) + memory->destroy(this->phirar5); + if (this->phirar6 != NULL) + memory->destroy(this->phirar6); // allocate memory for array that defines the potential - allocate_2d(this->phir, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); + memory->create(this->phir, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phir"); // allocate coeff memory - allocate_2d(this->phirar, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); - allocate_2d(this->phirar1, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); - allocate_2d(this->phirar2, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); - allocate_2d(this->phirar3, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); - allocate_2d(this->phirar4, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); - allocate_2d(this->phirar5, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); - allocate_2d(this->phirar6, this->nr, - (this->neltypes * (this->neltypes + 1)) / 2); + memory->create(this->phirar, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phirar"); + memory->create(this->phirar1, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phirar1"); + memory->create(this->phirar2, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phirar2"); + memory->create(this->phirar3, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phirar3"); + memory->create(this->phirar4, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phirar4"); + memory->create(this->phirar5, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phirar5"); + memory->create(this->phirar6, this->nr, + (this->neltypes * (this->neltypes + 1)) / 2, + "pair:phirar6"); // loop over pairs of element types nv2 = 0; @@ -600,7 +610,7 @@ MEAM::compute_reference_density(void) get_Zij2(&Z2, &arat, &scrn, this->lattce_meam[a][a], this->Cmin_meam[a][a][a], this->Cmax_meam[a][a][a]); rho0_2nn = - this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * (arat - 1)); + this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1)); rho0 = rho0 + Z2 * rho0_2nn * scrn; } @@ -667,8 +677,8 @@ MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, } else { a1 = r / this->re_meam[a][a] - 1.0; a2 = r / this->re_meam[b][b] - 1.0; - rhoa01 = this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * a1); - rhoa02 = this->rho0_meam[b] * fm_exp(-this->beta0_meam[b] * a2); + rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1); + rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2); if (latt == L12) { rho01 = 8 * rhoa01 + 4 * rhoa02; *t11av = (8 * t11 * rhoa01 + 4 * t12 * rhoa02) / rho01; @@ -803,14 +813,14 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* a1 = r / this->re_meam[a][a] - 1.0; a2 = r / this->re_meam[b][b] - 1.0; - rhoa01 = this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * a1); - rhoa11 = this->rho0_meam[a] * fm_exp(-this->beta1_meam[a] * a1); - rhoa21 = this->rho0_meam[a] * fm_exp(-this->beta2_meam[a] * a1); - rhoa31 = this->rho0_meam[a] * fm_exp(-this->beta3_meam[a] * a1); - rhoa02 = this->rho0_meam[b] * fm_exp(-this->beta0_meam[b] * a2); - rhoa12 = this->rho0_meam[b] * fm_exp(-this->beta1_meam[b] * a2); - rhoa22 = this->rho0_meam[b] * fm_exp(-this->beta2_meam[b] * a2); - rhoa32 = this->rho0_meam[b] * fm_exp(-this->beta3_meam[b] * a2); + rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1); + rhoa11 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1); + rhoa21 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1); + rhoa31 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1); + rhoa02 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2); + rhoa12 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta1_meam[b] * a2); + rhoa22 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta2_meam[b] * a2); + rhoa32 = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta3_meam[b] * a2); lat = this->lattce_meam[a][b]; @@ -889,8 +899,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* a1 = arat * r / this->re_meam[a][a] - 1.0; a2 = arat * r / this->re_meam[b][b] - 1.0; - rhoa01nn = this->rho0_meam[a] * fm_exp(-this->beta0_meam[a] * a1); - rhoa02nn = this->rho0_meam[b] * fm_exp(-this->beta0_meam[b] * a2); + rhoa01nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1); + rhoa02nn = this->rho0_meam[b] * MathSpecial::fm_exp(-this->beta0_meam[b] * a2); if (lat == L12) { // As usual, L12 thinks it's special; we need to be careful computing @@ -935,7 +945,7 @@ MEAM::zbl(double r, int z1, int z2) double result = 0.0; x = r / a; for (i = 0; i <= 3; i++) { - result = result + c[i] * fm_exp(-d[i] * x); + result = result + c[i] * MathSpecial::fm_exp(-d[i] * x); } if (r > 0.0) result = result * z1 * z2 / r * cc; @@ -962,12 +972,12 @@ MEAM::erose(double r, double re, double alpha, double Ec, double repuls, if (form == 1) result = -Ec * (1 + astar + (-attrac + repuls / r) * pow(astar, 3)) * - fm_exp(-astar); + MathSpecial::fm_exp(-astar); else if (form == 2) - result = -Ec * (1 + astar + a3 * pow(astar, 3)) * fm_exp(-astar); + result = -Ec * (1 + astar + a3 * pow(astar, 3)) * MathSpecial::fm_exp(-astar); else result = - -Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * fm_exp(-astar); + -Ec * (1 + astar + a3 * pow(astar, 3) / (r / re)) * MathSpecial::fm_exp(-astar); } return result; } @@ -1042,9 +1052,9 @@ MEAM::compute_phi(double rij, int elti, int eltj) ind = this->eltind[elti][eltj]; pp = rij * this->rdrar + 1.0; kk = (int)pp; - kk = MIN(kk, this->nrar - 1); + kk = std::min(kk, this->nrar - 1); pp = pp - kk; - pp = MIN(pp, 1.0); + pp = std::min(pp, 1.0); double result = ((arr2(this->phirar3, kk, ind) * pp + arr2(this->phirar2, kk, ind)) * pp + diff --git a/src/USER-MEAMC/meam_setup_global.cpp b/src/USER-MEAMC/meam_setup_global.cpp index 6d390e2e59..2df434d60f 100644 --- a/src/USER-MEAMC/meam_setup_global.cpp +++ b/src/USER-MEAMC/meam_setup_global.cpp @@ -18,7 +18,7 @@ using namespace LAMMPS_NS; // void -MEAM::meam_setup_global_(int* nelt, int* lat, double* z, int* ielement, double* atwt, +MEAM::meam_setup_global(int* nelt, int* lat, double* z, int* ielement, double* atwt, double* alpha, double* b0, double* b1, double* b2, double* b3, double* alat, double* esub, double* asub, double* t0, double* t1, double* t2, double* t3, diff --git a/src/USER-MEAMC/meam_setup_param.cpp b/src/USER-MEAMC/meam_setup_param.cpp index a04e5f8183..0af0cbb457 100644 --- a/src/USER-MEAMC/meam_setup_param.cpp +++ b/src/USER-MEAMC/meam_setup_param.cpp @@ -1,4 +1,5 @@ #include "meam.h" +#include using namespace LAMMPS_NS; // @@ -58,7 +59,7 @@ MEAM::meam_checkindex(int num, int lim, int nidx, int* idx /*idx(3)*/, int* ierr // 20 = bkgd_dyn void -MEAM::meam_setup_param_(int* which_p, double* value_p, int* nindex_p, +MEAM::meam_setup_param(int* which_p, double* value_p, int* nindex_p, int* index /*index(3)*/, int* errorflag) { //: index[0..2] @@ -149,8 +150,8 @@ MEAM::meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - i1 = MIN(index[0], index[1]); - i2 = MAX(index[0], index[1]); + i1 = std::min(index[0], index[1]); + i2 = std::max(index[0], index[1]); this->nn2_meam[i1][i2] = (int)value; break; @@ -218,8 +219,8 @@ MEAM::meam_setup_param_(int* which_p, double* value_p, int* nindex_p, meam_checkindex(2, maxelt, nindex, index, errorflag); if (*errorflag != 0) return; - i1 = MIN(index[0], index[1]); - i2 = MAX(index[0], index[1]); + i1 = std::min(index[0], index[1]); + i2 = std::max(index[0], index[1]); this->zbl_meam[i1][i2] = (int)value; break; diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index 94b0df7f7b..9f975bc4b4 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -64,7 +64,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp) nelements = 0; elements = NULL; mass = NULL; - meam_inst = new MEAM; + meam_inst = new MEAM(memory); // set comm size needed by this Pair @@ -238,7 +238,7 @@ void PairMEAMC::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_inst->meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0], + meam_inst->meam_dens_init(&ifort,&nmax,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], @@ -255,7 +255,7 @@ void PairMEAMC::compute(int eflag, int vflag) comm->reverse_comm_pair(this); - meam_inst->meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, + meam_inst->meam_dens_final(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, &eng_vdwl,eatom,&ntype,type,fmap, &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1, @@ -280,7 +280,7 @@ void PairMEAMC::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_inst->meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom, + meam_inst->meam_force(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom, &vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], @@ -369,7 +369,7 @@ void PairMEAMC::coeff(int narg, char **arg) // tell MEAM package that setup is done read_files(arg[2],arg[2+nelements+1]); - meam_inst->meam_setup_done_(&cutmax); + meam_inst->meam_setup_done(&cutmax); // read args that map atom types to MEAM elements // map[i] = which element the Ith atom type is, -1 if not mapped @@ -601,7 +601,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) // pass element parameters to MEAM package - meam_inst->meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, + meam_inst->meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, alat,esub,asub,t0,t1,t2,t3,rozero,ibar); // set element masses @@ -717,7 +717,7 @@ void PairMEAMC::read_files(char *globalfile, char *userfile) // pass single setting to MEAM package int errorflag = 0; - meam_inst->meam_setup_param_(&which,&value,&nindex,index,&errorflag); + meam_inst->meam_setup_param(&which,&value,&nindex,index,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); diff --git a/src/math_special.cpp b/src/math_special.cpp index 39487dd386..6777123ac8 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -508,6 +508,9 @@ static const double fm_exp2_p[] = { 1.51390680115615096133e3 }; +/* double precision constants */ +#define FM_DOUBLE_LOG2OFE 1.4426950408889634074 + double MathSpecial::exp2_x86(double x) { double ipart, fpart, px, qx; @@ -531,3 +534,14 @@ double MathSpecial::exp2_x86(double x) x = 1.0 + 2.0*(px/(qx-px)); return epart.f*x; } + +double MathSpecial::fm_exp(double x) +{ +#if defined(__BYTE_ORDER__) +#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__ + return exp2_x86(FM_DOUBLE_LOG2OFE * x); +#endif +#endif + return ::exp(x); +} + diff --git a/src/math_special.h b/src/math_special.h index e4b4998d54..8cd328f5fc 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -27,6 +27,9 @@ namespace MathSpecial { // fast 2**x function without argument checks for little endian CPUs extern double exp2_x86(double x); +// fast e**x function for little endian CPUs, falls back to libc on other platforms + extern double fm_exp(double x); + // scaled error function complement exp(x*x)*erfc(x) for coul/long styles static inline double my_erfcx(const double x)