Commit 55487047 authored by Sebastian Hütter's avatar Sebastian Hütter
Browse files

Move stateless functions to separate module, improve style 

- use static/const
- return instead of ptr-parameter, &ref if more than one return
- replace macros from header with inline functions
- remove useless/old comments
parent e0939ac7

src/USER-MEAMC/meam.h

+35 −26
Original line number Diff line number Diff line
@@ -2,6 +2,7 @@ 
#define LMP_MEAM_H



#include "memory.h"

#include <math.h>

#include <stdlib.h>



#define maxelt 5
@@ -107,10 +108,23 @@ public: 
  int maxneigh;

  double *scrfcn, *dscrfcn, *fcpair;



protected:

  // meam_funcs.cpp

  static double fcut(const double xi);

  static double dfcut(const double xi, double &dfc);



  static double dCfunc(const double rij2, const double rik2, const double rjk2);

  static void dCfunc2(const double rij2, const double rik2, const double rjk2, double &dCikj1, double &dCikj2);

  double G_gam(const double gamma, const int ibar, int &errorflag) const;

  double dG_gam(const double gamma, const int ibar, double &dG) const;

  static double zbl(const double r, const int z1, const int z2);

  static double erose(const double r, const double re, const double alpha, const double Ec, const double repuls, const double attrac, const int form);



  static void get_shpfcn(const lattice_t latt, double (&s)[3]);

  static int get_Zij(const lattice_t latt);

  static int get_Zij2(const lattice_t latt, const double cmin, const double cmax, double &a, double &S);

protected:

  void meam_checkindex(int, int, int, int*, int*);

  void G_gam(double, int, double*, int*);

  void dG_gam(double, int, double*, double*);

  void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,

                 int* firstneigh, int numneigh_full, int* firstneigh_full, int ntype, int* type, int* fmap);

  void screen(int i, int j, double** x, double rijsq, double* sij, int numneigh_full, int* firstneigh_full,
@@ -119,24 +133,15 @@ protected: 
                 double* scrfcn, double* fcpair);

  void dsij(int i, int j, int k, int jn, int numneigh, double rij2, double* dsij1, double* dsij2, int ntype,

            int* type, int* fmap, double** x, double* scrfcn, double* fcpair);

  void fcut(double, double*);

  void dfcut(double, double*, double*);

  void dCfunc(double, double, double, double*);

  void dCfunc2(double, double, double, double*, double*);



  void alloyparams();

  void compute_pair_meam();

  double phi_meam(double, int, int);

  void compute_reference_density();

  void get_shpfcn(double*, lattice_t);

  void get_tavref(double*, double*, double*, double*, double*, double*, double, double, double, double,

                  double, double, double, int, int, lattice_t);

  void get_Zij(int*, lattice_t);

  void get_Zij2(int*, double*, double*, lattice_t, double, double);

  void get_sijk(double, int, int, int, double*);

  void get_densref(double, int, int, double*, double*, double*, double*, double*, double*, double*, double*);

  double zbl(double, int, int);

  double erose(double, double, double, double, double, double, int);

  void interpolate_meam(int);

  double compute_phi(double, int, int);
@@ -160,20 +165,24 @@ public: 


// Functions we need for compat



#define iszero(f) (fabs(f) < 1e-20)

static inline bool iszero(const double f) {

  return fabs(f) < 1e-20;

}



#define setall2d(arr, v)                                                                                     \

  {                                                                                                          \

    for (int __i = 0; __i < maxelt; __i++)                                                                   \

      for (int __j = 0; __j < maxelt; __j++)                                                                 \

        arr[__i][__j] = v;                                                                                   \

template <typename TYPE, size_t maxi, size_t maxj>

static inline void setall2d(TYPE (&arr)[maxi][maxj], const TYPE v) {

  for (size_t i = 0; i < maxi; i++)

    for (size_t j = 0; j < maxj; j++)

      arr[i][j] = v;

}

#define setall3d(arr, v)                                                                                     \

  {                                                                                                          \

    for (int __i = 0; __i < maxelt; __i++)                                                                   \

      for (int __j = 0; __j < maxelt; __j++)                                                                 \

        for (int __k = 0; __k < maxelt; __k++)                                                               \

          arr[__i][__j][__k] = v;                                                                            \



template <typename TYPE, size_t maxi, size_t maxj, size_t maxk>

static inline void setall3d(TYPE (&arr)[maxi][maxj][maxk], const TYPE v) {

  for (size_t i = 0; i < maxi; i++)

    for (size_t j = 0; j < maxj; j++)

      for (size_t k = 0; k < maxk; k++)

        arr[i][j][k] = v;

}



};

#endif

src/USER-MEAMC/meam_dens_final.cpp

+5 −107
Original line number Diff line number Diff line 
#include "meam.h"

#include "math_special.h"

#include <math.h>



using namespace LAMMPS_NS;

// Extern "C" declaration has the form:

//

//  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 *);

//

// Call from pair_meam.cpp has the form:

//

//  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,

//	     dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);

//



void

MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_atom, double* eng_vdwl,
@@ -72,10 +55,10 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ 


      Z = this->Z_meam[elti];



      G_gam(gamma[i], this->ibar_meam[elti], &G, errorflag);

      G = G_gam(gamma[i], this->ibar_meam[elti], *errorflag);

      if (*errorflag != 0)

        return;

      get_shpfcn(shp, this->lattce_meam[elti][elti]);

      get_shpfcn(this->lattce_meam[elti][elti], shp);

      if (this->ibar_meam[elti] <= 0) {

        Gbar = 1.0;

        dGbar = 0.0;
@@ -86,7 +69,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ 
          gam = (this->t1_meam[elti] * shp[0] + this->t2_meam[elti] * shp[1] + this->t3_meam[elti] * shp[2]) /

                (Z * Z);

        }

        G_gam(gam, this->ibar_meam[elti], &Gbar, errorflag);

        Gbar = G_gam(gam, this->ibar_meam[elti], *errorflag);

      }

      rho[i] = rho0[i] * G;
@@ -96,7 +79,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ 
          dGbar = 0.0;

        } else {

          gam = (t_ave[i][0] * shp[0] + t_ave[i][1] * shp[1] + t_ave[i][2] * shp[2]) / (Z * Z);

          dG_gam(gam, this->ibar_meam[elti], &Gbar, &dGbar);

          Gbar = dG_gam(gam, this->ibar_meam[elti], dGbar);

        }

        rho_bkgd = this->rho0_meam[elti] * Z * Gbar;

      } else {
@@ -109,7 +92,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ 
      rhob = rho[i] / rho_bkgd;

      denom = 1.0 / rho_bkgd;



      dG_gam(gamma[i], this->ibar_meam[elti], &G, &dG);

      G = dG_gam(gamma[i], this->ibar_meam[elti], dG);



      dgamma1[i] = (G - 2 * dG * gamma[i]) * denom;
@@ -163,88 +146,3 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ 
  }

}

// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



void

MEAM::G_gam(double gamma, int ibar, double* G, int* errorflag)

{

  //     Compute G(gamma) based on selection flag ibar:

  //   0 => G = sqrt(1+gamma)

  //   1 => G = exp(gamma/2)

  //   2 => not implemented

  //   3 => G = 2/(1+exp(-gamma))

  //   4 => G = sqrt(1+gamma)

  //  -5 => G = +-sqrt(abs(1+gamma))



  double gsmooth_switchpoint;

  if (ibar == 0 || ibar == 4) {

    gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);

    if (gamma < gsmooth_switchpoint) {

      //         e.g. gsmooth_factor is 99, {:

      //         gsmooth_switchpoint = -0.99

      //         G = 0.01*(-0.99/gamma)**99

      *G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);

      *G = sqrt(*G);

    } else {

      *G = sqrt(1.0 + gamma);

    }

  } else if (ibar == 1) {

    *G = MathSpecial::fm_exp(gamma / 2.0);

  } else if (ibar == 3) {

    *G = 2.0 / (1.0 + exp(-gamma));

  } else if (ibar == -5) {

    if ((1.0 + gamma) >= 0) {

      *G = sqrt(1.0 + gamma);

    } else {

      *G = -sqrt(-1.0 - gamma);

    }

  } else {

    *errorflag = 1;

  }

}



// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



void

MEAM::dG_gam(double gamma, int ibar, double* G, double* dG)

{

  // Compute G(gamma) and dG(gamma) based on selection flag ibar:

  //   0 => G = sqrt(1+gamma)

  //   1 => G = MathSpecial::fm_exp(gamma/2)

  //   2 => not implemented

  //   3 => G = 2/(1+MathSpecial::fm_exp(-gamma))

  //   4 => G = sqrt(1+gamma)

  //  -5 => G = +-sqrt(abs(1+gamma))

  double gsmooth_switchpoint;



  if (ibar == 0 || ibar == 4) {

    gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor + 1);

    if (gamma < gsmooth_switchpoint) {

      //         e.g. gsmooth_factor is 99, {:

      //         gsmooth_switchpoint = -0.99

      //         G = 0.01*(-0.99/gamma)**99

      *G = 1 / (gsmooth_factor + 1) * pow((gsmooth_switchpoint / gamma), gsmooth_factor);

      *G = sqrt(*G);

      *dG = -gsmooth_factor * *G / (2.0 * gamma);

    } else {

      *G = sqrt(1.0 + gamma);

      *dG = 1.0 / (2.0 * *G);

    }

  } else if (ibar == 1) {

    *G = MathSpecial::fm_exp(gamma / 2.0);

    *dG = *G / 2.0;

  } else if (ibar == 3) {

    *G = 2.0 / (1.0 + MathSpecial::fm_exp(-gamma));

    *dG = *G * (2.0 - *G) / 2;

  } else if (ibar == -5) {

    if ((1.0 + gamma) >= 0) {

      *G = sqrt(1.0 + gamma);

      *dG = 1.0 / (2.0 * *G);

    } else {

      *G = -sqrt(-1.0 - gamma);

      *dG = -1.0 / (2.0 * *G);

    }

  }

}



// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

src/USER-MEAMC/meam_dens_init.cpp

+6 −112
Original line number Diff line number Diff line 
#include "meam.h"

#include "math_special.h"

#include <math.h>



using namespace LAMMPS_NS;
@@ -77,24 +76,6 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh) 
  }

}



//     Extern "C" declaration has the form:

//

//  void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *,

//  double *,

//		 int *, int *, int *, int *,

//		 double *, double *, double *, double *, double *, double *,

//		 double *, double *, double *, double *, double *, int *);

//

//

//     Call from pair_meam.cpp has the form:

//

//    meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0],

//	       &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],

//	       &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],

//	       rho0,&arho1[0][0],&arho2[0][0],arho2b,

//	       &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag);

//



void

MEAM::meam_dens_init(int i, int ntype, int* type, int* fmap, double** x, int numneigh, int* firstneigh,

                     int numneigh_full, int* firstneigh_full, int fnoffset, int* errorflag)
@@ -158,7 +139,7 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** 
        } else {

          rnorm = (this->rc_meam - rij) * drinv;

          screen(i, j, x, rij2, &sij, numneigh_full, firstneigh_full, ntype, type, fmap);

          dfcut(rnorm, &fc, &dfc);

          fc = dfcut(rnorm, dfc);

          fcij = fc;

          dfcij = dfc * drinv;

        }
@@ -213,9 +194,9 @@ MEAM::getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** 
          } else {

            delc = Cmax - Cmin;

            cikj = (cikj - Cmin) / delc;

            dfcut(cikj, &sikj, &dfikj);

            sikj = dfcut(cikj, dfikj);

            coef1 = dfikj / (delc * sikj);

            dCfunc(rij2, rik2, rjk2, &dCikj);

            dCikj = dCfunc(rij2, rik2, rjk2);

            dscrfcn[jn] = dscrfcn[jn] + coef1 * dCikj;

          }

        }
@@ -394,7 +375,7 @@ MEAM::screen(int i, int j, double** x, double rijsq, double* sij, int numneigh_f 
    } else {

      delc = Cmax - Cmin;

      cikj = (cikj - Cmin) / delc;

      fcut(cikj, &sikj);

      sikj = fcut(cikj);

    }

    *sij = *sij * sikj;

  }
@@ -449,8 +430,8 @@ MEAM::dsij(int i, int j, int k, int jn, int numneigh, double rij2, double* dsij1 
          cikj = (2.0 * (xik + xjk) + a - 2.0) / a;

          if (cikj >= Cmin && cikj <= Cmax) {

            cikj = (cikj - Cmin) / delc;

            dfcut(cikj, &sikj, &dfc);

            dCfunc2(rij2, rik2, rjk2, &dCikj1, &dCikj2);

            sikj = dfcut(cikj, dfc);

            dCfunc2(rij2, rik2, rjk2, dCikj1, dCikj2);

            a = sij / delc * dfc / sikj;

            *dsij1 = a * dCikj1;

            *dsij2 = a * dCikj2;
@@ -460,90 +441,3 @@ MEAM::dsij(int i, int j, int k, int jn, int numneigh, double rij2, double* dsij1 
    }

  }

}
 


// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



void

MEAM::fcut(double xi, double* fc)

{

  //     cutoff function

  double a;

  if (xi >= 1.0)

    *fc = 1.0;

  else if (xi <= 0.0)

    *fc = 0.0;

  else {

    a = 1.0 - xi;

    a = a * a;

    a = a * a;

    a = 1.0 - a;

    *fc = a * a;

    //     fc = xi

  }

}



// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



void

MEAM::dfcut(double xi, double* fc, double* dfc)

{

  //     cutoff function and its derivative

  double a, a3, a4;

  if (xi >= 1.0) {

    *fc = 1.0;

    *dfc = 0.0;

  } else if (xi <= 0.0) {

    *fc = 0.0;

    *dfc = 0.0;

  } else {

    a = 1.0 - xi;

    a3 = a * a * a;

    a4 = a * a3;

    *fc = pow((1.0 - a4), 2);

    *dfc = 8 * (1.0 - a4) * a3;

    //     fc = xi

    //     dfc = 1.d0

  }

}



// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



void

MEAM::dCfunc(double rij2, double rik2, double rjk2, double* dCikj)

{

  //     Inputs: rij,rij2,rik2,rjk2

  //     Outputs: dCikj = derivative of Cikj w.r.t. rij

  double rij4, a, b, denom;



  rij4 = rij2 * rij2;

  a = rik2 - rjk2;

  b = rik2 + rjk2;

  denom = rij4 - a * a;

  denom = denom * denom;

  *dCikj = -4 * (-2 * rij2 * a * a + rij4 * b + a * a * b) / denom;

}



// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



void

MEAM::dCfunc2(double rij2, double rik2, double rjk2, double* dCikj1, double* dCikj2)

{

  //     Inputs: rij,rij2,rik2,rjk2

  //     Outputs: dCikj1 = derivative of Cikj w.r.t. rik

  //     dCikj2 = derivative of Cikj w.r.t. rjk

  double rij4, rik4, rjk4, a, b, denom;



  rij4 = rij2 * rij2;

  rik4 = rik2 * rik2;

  rjk4 = rjk2 * rjk2;

  a = rik2 - rjk2;

  b = rik2 + rjk2;

  denom = rij4 - a * a;

  denom = denom * denom;

  *dCikj1 = 4 * rij2 * (rij4 + rik4 + 2 * rik2 * rjk2 - 3 * rjk4 - 2 * rij2 * a) / denom;

  *dCikj2 = 4 * rij2 * (rij4 - 3 * rik4 + 2 * rik2 * rjk2 + rjk4 + 2 * rij2 * a) / denom;



  (void)(b);

}



// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

src/USER-MEAMC/meam_force.cpp

+3 −23
Original line number Diff line number Diff line 
#include "meam.h"

#include "math_special.h"

#include <algorithm>

#include <math.h>



using namespace LAMMPS_NS;

//     Extern "C" declaration has the form:

//

//  void meam_force_(int *, int *, int *, 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

//*);

//

//     Call from pair_meam.cpp has the form:

//

//    meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom,

//              &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],

//	       &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],

//	       &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],

//	       dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,

//	       &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],

//	       &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag);

//





void

MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int vflag_atom, double* eng_vdwl,
@@ -313,8 +293,8 @@ MEAM::meam_force(int i, int eflag_either, int eflag_global, int eflag_atom, int 
          }



          //     Compute derivatives of total density wrt rij, sij and rij(3)

          get_shpfcn(shpi, this->lattce_meam[elti][elti]);

          get_shpfcn(shpj, this->lattce_meam[eltj][eltj]);

          get_shpfcn(this->lattce_meam[elti][elti], shpi);

          get_shpfcn(this->lattce_meam[eltj][eltj], shpj);

          drhodr1 = dgamma1[i] * drho0dr1 +

                    dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 +

                                  dt3dr1 * rho3[i] + t3i * drho3dr1) -

src/USER-MEAMC/meam_funcs.cpp

0 → 100644
+384 −0

File added.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE