forked from lijiext/lammps
Improve C++-ness, eliminate some macros
- fm_exp moved to math_special (exp2 was already there) - use std::min/max template instead of macros - use memory->create for dynamic arrays (still 1-indexed with macro) - remove _ from function names, adjust method visibility
This commit is contained in:
Sebastian Hütter
parent
fea28d8028
commit
9f852f5f58
|
|
@ -1,135 +0,0 @@
|
|||
extern "C" {
|
||||
/*
|
||||
Copyright (c) 2012,2013 Axel Kohlmeyer <akohlmey@gmail.com>
|
||||
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 <organization> 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 <COPYRIGHT HOLDER> 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 <math.h>
|
||||
#include <stdint.h>
|
||||
|
||||
/* 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:
|
||||
*/
|
||||
}
|
||||
|
|
@ -2,31 +2,25 @@
|
|||
#define LMP_MEAM_H
|
||||
|
||||
#include <stdlib.h>
|
||||
#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]
|
||||
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
|
@ -1,5 +1,6 @@
|
|||
#include "meam.h"
|
||||
#include <math.h>
|
||||
#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) {
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
#include "meam.h"
|
||||
#include <math.h>
|
||||
#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];
|
||||
|
|
|
|||
|
|
@ -1,5 +1,7 @@
|
|||
#include "meam.h"
|
||||
#include <math.h>
|
||||
#include <algorithm>
|
||||
#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;
|
||||
|
|
|
|||
|
|
@ -1,5 +1,7 @@
|
|||
#include "meam.h"
|
||||
#include <math.h>
|
||||
#include <algorithm>
|
||||
#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 +
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -1,4 +1,5 @@
|
|||
#include "meam.h"
|
||||
#include <algorithm>
|
||||
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;
|
||||
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Reference in New Issue