Commit 81fb0d61 authored by Sungkwang Mun's avatar Sungkwang Mun Committed by Sebastian Hütter
Browse files

* This commit includes the addition of new reference structures such as 

ch4: methane-like structure only for binary system.
dia3: diamond structure with primary 1NN and secondary 3NN inteation
tri: H2O-like structure that has an angle
zig: zigzag structure with a uniform angle
lin: linear structure (180 degree angle)

** tri, zig, and lin reference structures require angle information (in degree)
such as the following.
   theta = 109.5
parent ce05ed4c

src/USER-MEAMC/meam.h

+27 −5
Original line number Diff line number Diff line
@@ -3,13 +3,14 @@ 


#include <cmath>

#include <cstring>

#include "math_const.h"



#define maxelt 5



namespace LAMMPS_NS {

class Memory;



typedef enum { FCC, BCC, HCP, DIM, DIA, B1, C11, L12, B2 } lattice_t;

typedef enum { FCC, BCC, HCP, DIM, DIA, DIA3, B1, C11, L12, B2, CH4, LIN, ZIG, TRI } lattice_t;



class MEAM

{
@@ -63,6 +64,10 @@ private: 
  // nr,dr = pair function discretization parameters

  // nrar,rdrar = spline coeff array parameters



  // theta = angle between three atoms in line, zigzag, and trimer reference structures

  // stheta_meam = sin(theta/2) in radian used in line, zigzag, and trimer reference structures

  // ctheta_meam = cos(theta/2) in radian used in line, zigzag, and trimer reference structures



  double Ec_meam[maxelt][maxelt], re_meam[maxelt][maxelt];

  double A_meam[maxelt], alpha_meam[maxelt][maxelt], rho0_meam[maxelt];

  double delta_meam[maxelt][maxelt];
@@ -106,6 +111,10 @@ public: 
  int maxneigh;

  double *scrfcn, *dscrfcn, *fcpair;



  //angle for trimer, zigzag, line reference structures

  double stheta_meam[maxelt][maxelt];

  double ctheta_meam[maxelt][maxelt];



protected:

  // meam_funcs.cpp
@@ -189,8 +198,12 @@ protected: 
  double embedding(const double A, const double Ec, const double rhobar, double& dF) const;

  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_Zij2(const lattice_t latt, const double cmin, const double cmax, double &a, double &S);

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



  static int get_Zij2(const lattice_t latt, const double cmin, const double cmax,

                      const double sthe, double &a, double &S);

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



protected:

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

  void getscreen(int i, double* scrfcn, double* dscrfcn, double* fcpair, double** x, int numneigh,
@@ -201,6 +214,7 @@ protected: 
  void alloyparams();

  void compute_pair_meam();

  double phi_meam(double, int, int);

  const double phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat);

  void compute_reference_density();

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

                  double, double, double, int, int, lattice_t);
@@ -221,6 +235,10 @@ public: 
    else if (strcmp(str,"hcp") == 0) lat = HCP;

    else if (strcmp(str,"dim") == 0) lat = DIM;

    else if (strcmp(str,"dia") == 0) lat = DIA;

    else if (strcmp(str,"dia3") == 0) lat = DIA3;

    else if (strcmp(str,"lin") == 0) lat = LIN;

    else if (strcmp(str,"zig") == 0) lat = ZIG;

    else if (strcmp(str,"tri") == 0) lat = TRI;

    else {

      if (single)

        return false;
@@ -229,6 +247,10 @@ public: 
      else if (strcmp(str,"c11") == 0) lat = C11;

      else if (strcmp(str,"l12") == 0) lat = L12;

      else if (strcmp(str,"b2")  == 0) lat = B2;

      else if (strcmp(str,"ch4")  == 0) lat = CH4;

      else if (strcmp(str,"lin")  == 0) lat =LIN;

      else if (strcmp(str,"zig")  == 0) lat = ZIG;

      else if (strcmp(str,"tri")  == 0) lat = TRI;

      else return false;

    }

    return true;

src/USER-MEAMC/meam_dens_final.cpp

+3 −1
Original line number Diff line number Diff line
@@ -59,7 +59,9 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_ 
      G = G_gam(gamma[i], this->ibar_meam[elti], errorflag);

      if (errorflag != 0)

        return;

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



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



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

        Gbar = 1.0;

        dGbar = 0.0;

src/USER-MEAMC/meam_force.cpp

+3 −2
Original line number Diff line number Diff line
@@ -290,8 +290,9 @@ 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(this->lattce_meam[elti][elti], shpi);

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

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

        get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], 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

+82 −8
Original line number Diff line number Diff line
@@ -198,7 +198,7 @@ MEAM::erose(const double r, const double re, const double alpha, const double Ec 
// Shape factors for various configurations

//

void

MEAM::get_shpfcn(const lattice_t latt, double (&s)[3])

MEAM::get_shpfcn(const lattice_t latt, const double sthe, const double cthe, double (&s)[3])

{

  switch (latt) {

    case FCC:
@@ -214,7 +214,9 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3]) 
      s[1] = 0.0;

      s[2] = 1.0 / 3.0;

      break;

    case CH4: // CH4 actually needs shape factor for diamond for C, dimer for H

    case DIA:

    case DIA3:

      s[0] = 0.0;

      s[1] = 0.0;

      s[2] = 32.0 / 9.0;
@@ -222,8 +224,20 @@ MEAM::get_shpfcn(const lattice_t latt, double (&s)[3]) 
    case DIM:

      s[0] = 1.0;

      s[1] = 2.0 / 3.0;

      //        s(3) = 1.d0

      s[2] = 0.40;

      //        s(3) = 1.d0 // this should be 0.4 unless (1-legendre) is multiplied in the density calc.

      s[2] = 0.40; // this is (1-legendre) where legendre = 0.6 in dynamo is accounted.

      break;

    case LIN: //linear, theta being 180

      s[0] = 0.0;

      s[1] = 8.0 / 3.0; // 4*(co**4 + si**4 - 1.0/3.0) in zig become 4*(1-1/3)

      s[2] = 0.0;

      break;

    case ZIG: //zig-zag

    case TRI: //trimer e.g. H2O

      s[0] = 4.0*pow(cthe,2);

      s[1] = 4.0*(pow(cthe,4) + pow(sthe,4) - 1.0/3.0);

      s[2] = 4.0*(pow(cthe,2) * (3*pow(sthe,4) + pow(cthe,4)));

      s[2] = s[2] - 0.6*s[0]; //legend in dyn, 0.6 is default value.

      break;

    default:

      s[0] = 0.0;
@@ -245,6 +259,7 @@ MEAM::get_Zij(const lattice_t latt) 
    case HCP:

      return 12;

    case DIA:

    case DIA3:

      return 4;

    case DIM:

      return 1;
@@ -256,6 +271,12 @@ MEAM::get_Zij(const lattice_t latt) 
      return 12;

    case B2:

      return 8;

    case CH4: // DYNAMO currenly implemented this way while it needs two Z values, 4 and 1

      return 4;

    case LIN:

    case ZIG:

    case TRI:

      return 2;

      //        call error('Lattice not defined in get_Zij.')

  }

  return 0;
@@ -263,11 +284,13 @@ MEAM::get_Zij(const lattice_t latt) 


//-----------------------------------------------------------------------------

// Number of second neighbors for the reference structure

//   a = distance ratio R1/R2

//   S = second neighbor screening function

//   a = distance ratio R1/R2 (a2nn in dynamo)

//   numscr = number of atoms that screen the 2NN bond

//   S = second neighbor screening function (xfac, a part of b2nn in dynamo)

//

int

MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, double& a, double& S)

MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax,

                const double stheta, double& a, double& S)

{



  double C, x, sijk;
@@ -299,7 +322,7 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl 
    numscr = 2;

    break;



  case DIA:

  case DIA: // 2NN

    Zij2 = 12;

    a = sqrt(8.0 / 3.0);

    numscr = 1;
@@ -308,12 +331,35 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl 
    }

    break;



  case DIA3: // 3NN

    Zij2 = 12;

    a = sqrt(11.0 / 3.0);

    numscr = 4;

    if (cmin < 0.500001) {

        //          call error('can not do 2NN MEAM for dia')

    }

    break;



  case CH4: //does not have 2nn structure so it returns 0

  case LIN: //line

  case DIM:

    //        this really shouldn't be allowed; make sure screening is zero

    a = 1.0;

    S = 0.0;

    return 0;



  case TRI: //TRI

    Zij2 = 1;

    a = 2.0*stheta;

    numscr=2;

    break;



  case ZIG: //zig-zag

    Zij2 = 2;

    a = 2.0*stheta;

    numscr=1;

    break;



  case L12:

    Zij2 = 6;

    a = sqrt(2.0);
@@ -335,10 +381,38 @@ MEAM::get_Zij2(const lattice_t latt, const double cmin, const double cmax, doubl 
  }



  // Compute screening for each first neighbor

  if (latt==DIA3){ // special case for 3NN diamond structure

	C = 1.0;

  } else {

	C = 4.0 / (a * a) - 1.0;

  }

  x = (C - cmin) / (cmax - cmin);

  sijk = fcut(x);

  // There are numscr first neighbors screening the second neighbors

  S = MathSpecial::powint(sijk, numscr);

  return Zij2;

}



int

MEAM::get_Zij2_b2nn(const lattice_t latt, const double cmin, const double cmax, double &S)

{



  double x, sijk, C;

  int numscr = 0, Zij2 = 0;

  switch (latt) {

  case ZIG: //zig-zag for b11s and b22s term

  case TRI: //trimer for b11s

    Zij2 = 2;

    numscr=1;

    break;

  default:

    // unknown lattic flag in get Zij

    //        call error('Lattice not defined in get_Zij.')

    break;

  }

  C = 1.0;

  x = (C - cmin) / (cmax - cmin);

  sijk = fcut(x);

  S = MathSpecial::powint(sijk, numscr);

  return Zij2;

}
 
 No newline at end of file

src/USER-MEAMC/meam_setup_done.cpp

+152 −40
Original line number Diff line number Diff line
@@ -103,6 +103,9 @@ MEAM::alloyparams(void) 
        this->alpha_meam[i][j] = this->alpha_meam[j][i];

        this->lattce_meam[i][j] = this->lattce_meam[j][i];

        this->nn2_meam[i][j] = this->nn2_meam[j][i];

	    // theta for lin,tri,zig references

        this->stheta_meam[i][j] = this->stheta_meam[j][i];

        this->ctheta_meam[i][j] = this->ctheta_meam[j][i];

        // If i<j and term is unset, use default values (e.g. mean of i-i and

        // j-j)

      } else if (j > i) {
@@ -161,7 +164,7 @@ void 
MEAM::compute_pair_meam(void)

{



  double r;

  double r, b2nn, phi_val;

  int j, a, b, nv2;

  double astar, frac, phizbl;

  int n, Z1, Z2;
@@ -215,7 +218,7 @@ MEAM::compute_pair_meam(void) 
        if (this->nn2_meam[a][b] == 1) {

          Z1 = get_Zij(this->lattce_meam[a][b]);

          Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],

                   this->Cmax_meam[a][a][b], arat, scrn);

                     this->Cmax_meam[a][a][b], this->stheta_meam[a][b], arat, scrn);



          //     The B1, B2,  and L12 cases with NN2 have a trick to them; we need to

          //     compute the contributions from second nearest neighbors, like a-a
@@ -229,33 +232,26 @@ MEAM::compute_pair_meam(void) 
            phiaa = phi_meam(rarat, a, a);

            Z1 = get_Zij(this->lattce_meam[a][a]);

            Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],

                     this->Cmax_meam[a][a][a], arat, scrn);

            if (scrn > 0.0) {

              for (n = 1; n <= 10; n++) {

                phiaa = phiaa + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), a, a);

              }

            }

                     this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);

            phiaa+= phi_meam_series(scrn, Z1, Z2, a, a, rarat, arat);



            //               phi_bb

            phibb = phi_meam(rarat, b, b);

            Z1 = get_Zij(this->lattce_meam[b][b]);

            Z2 = get_Zij2(this->lattce_meam[b][b], this->Cmin_meam[b][b][b],

                     this->Cmax_meam[b][b][b], arat, scrn);

            if (scrn > 0.0) {

              for (n = 1; n <= 10; n++) {

                phibb = phibb + pow((-Z2 * scrn / Z1), n) * phi_meam(rarat * pow(arat, n), b, b);

              }

            }

                     this->Cmax_meam[b][b][b], this->stheta_meam[b][b], arat, scrn);

            phibb+= phi_meam_series(scrn, Z1, Z2, b, b, rarat, arat);



            if (this->lattce_meam[a][b] == B1 || this->lattce_meam[a][b] == B2 ||

                this->lattce_meam[a][b] == DIA) {

              //     Add contributions to the B1 or B2 potential

              Z1 = get_Zij(this->lattce_meam[a][b]);

              Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[a][a][b],

                       this->Cmax_meam[a][a][b], arat, scrn);

                       this->Cmax_meam[a][a][b], this->stheta_meam[a][b],  arat, scrn);

              this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn / (2 * Z1) * phiaa;

              Z2 = get_Zij2(this->lattce_meam[a][b], this->Cmin_meam[b][b][a],

                       this->Cmax_meam[b][b][a], arat, scrn2);

                       this->Cmax_meam[b][b][a], this->stheta_meam[a][b], arat, scrn2);



              this->phir[nv2][j] = this->phir[nv2][j] - Z2 * scrn2 / (2 * Z1) * phibb;



            } else if (this->lattce_meam[a][b] == L12) {
@@ -277,10 +273,7 @@ MEAM::compute_pair_meam(void) 
            }



          } else {

            for (n = 1; n <= 10; n++) {

              this->phir[nv2][j] =

                this->phir[nv2][j] + pow((-Z2 * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, b);

            }

            this->phir[nv2][j]+= phi_meam_series(scrn, Z1, Z2, a, b, r, arat);

          }

        }
@@ -325,11 +318,12 @@ MEAM::phi_meam(double r, int a, int b) 
  double rho02, rho12, rho22, rho32;

  double scalfac, phiaa, phibb;

  double Eu;

  double arat, scrn /*unused:,scrn2*/;

  double arat, scrn, scrn2;

  int Z12, errorflag;

  int n, Z1nn;

  int n, Z1nn, Z2nn;

  lattice_t latta /*unused:,lattb*/;

  double rho_bkgd1, rho_bkgd2;

  double b11s, b22s;



  double phi_m = 0.0;
@@ -403,14 +397,14 @@ MEAM::phi_meam(double r, int a, int b) 
      if (this->ibar_meam[a] <= 0)

        G1 = 1.0;

      else {

        get_shpfcn(this->lattce_meam[a][a], s1);

        get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], s1);

        Gam1 = (s1[0] * t11av + s1[1] * t21av + s1[2] * t31av) / (Z1 * Z1);

        G1 = G_gam(Gam1, this->ibar_meam[a], errorflag);

      }

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

        G2 = 1.0;

      else {

        get_shpfcn(this->lattce_meam[b][b], s2);

        get_shpfcn(this->lattce_meam[b][b], this->stheta_meam[b][b], this->ctheta_meam[b][b],  s2);

        Gam2 = (s2[0] * t12av + s2[1] * t22av + s2[2] * t32av) / (Z2 * Z2);

        G2 = G_gam(Gam2, this->ibar_meam[b], errorflag);

      }
@@ -470,20 +464,46 @@ MEAM::phi_meam(double r, int a, int b) 
  } else if (this->lattce_meam[a][b] == L12) {

    phiaa = phi_meam(r, a, a);

    //       account for second neighbor a-a potential here...

    Z1nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],

             this->Cmax_meam[a][a][a], arat, scrn);

    if (scrn > 0.0) {

      for (n = 1; n <= 10; n++) {

        phiaa = phiaa + pow((-Z1nn * scrn / Z1), n) * phi_meam(r * pow(arat, n), a, a);

    Z1nn = get_Zij(this->lattce_meam[a][a]);

    Z2nn = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],

             this->Cmax_meam[a][a][a], this->stheta_meam[a][b], arat, scrn);





    phiaa += phi_meam_series(scrn, Z1nn, Z2nn, a, a, r, arat);

    phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;



  } else if (this->lattce_meam[a][b] == CH4) {

    phi_m = (5 * Eu - F1 - 4*F2)/4;



  } else if (this->lattce_meam[a][b] == ZIG){

      if (a==b){

        phi_m = (2 * Eu - F1 - F2) / Z12;

      } else{

        Z1 = get_Zij(this->lattce_meam[a][b]);

        Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);

        b11s = -Z2/Z1*scrn;

        Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], scrn2);

        b22s = -Z2/Z1*scrn2;



        phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);

        phibb = phi_meam(2.0*this->stheta_meam[a][b]*r, b, b);

        phi_m = (2.0*Eu - F1 - F2 + phiaa*b11s + phibb*b22s) / Z12;

      }



  } else if (this->lattce_meam[a][b] == TRI) {

      if (a==b){

        phi_m = (3.0*Eu - 2.0*F1 - F2) / Z12;

     } else {

        Z1 = get_Zij(this->lattce_meam[a][b]);

        Z2 = get_Zij2_b2nn(this->lattce_meam[a][b], this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], scrn);

        b11s = -Z2/Z1*scrn;

        phiaa = phi_meam(2.0*this->stheta_meam[a][b]*r, a, a);

        phi_m = (3.0*Eu - 2.0*F1 - F2 + phiaa*b11s) / Z12;

      }

    phi_m = Eu / 3.0 - F1 / 4.0 - F2 / 12.0 - phiaa;



  } else {

    //

    // potential is computed from Rose function and embedding energy

    phi_m = (2 * Eu - F1 - F2) / Z12;

    //

  }



  // if r = 0, just return 0
@@ -494,6 +514,31 @@ MEAM::phi_meam(double r, int a, int b) 
  return phi_m;

}



//----------------------------------------------------------------------c

// Compute 2NN series terms for phi

//   To avoid nan values of phir due to rapid decrease of b2nn^n or/and

//   argument of phi_meam, i.e. r*arat^n, in some cases (3NN dia with low Cmin value)

//

const double

MEAM::phi_meam_series(const double scrn, const int Z1, const int Z2, const int a, const int b, const double r, const double arat)

{

  double phi_sum = 0.0;

  double b2nn, phi_val;

  if (scrn > 0.0) {

    b2nn = -Z2*scrn/Z1;

    for (int n = 1; n <= 10; n++) {

      phi_val = MathSpecial::powint(b2nn,n) * phi_meam(r * MathSpecial::powint(arat, n), a, b);

      if (iszero(phi_val)) {

        // once either term becomes zero at some point, all folliwng will also be zero

        // necessary to avoid numerical error (nan or infty) due to exponential decay in phi_meam

        break;

      }

      phi_sum += phi_val;

    }

  }

  return phi_sum;

}



//----------------------------------------------------------------------c

// Compute background density for reference structure of each element

void
@@ -509,7 +554,7 @@ MEAM::compute_reference_density(void) 
    if (this->ibar_meam[a] <= 0)

      Gbar = 1.0;

    else {

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

      get_shpfcn(this->lattce_meam[a][a], this->stheta_meam[a][a], this->ctheta_meam[a][a], shp);

      gam = (this->t1_meam[a] * shp[0] + this->t2_meam[a] * shp[1] + this->t3_meam[a] * shp[2]) / (Z * Z);

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

    }
@@ -524,7 +569,7 @@ MEAM::compute_reference_density(void) 
    //     screening)

    if (this->nn2_meam[a][a] == 1) {

      Z2 = get_Zij2(this->lattce_meam[a][a], this->Cmin_meam[a][a][a],

               this->Cmax_meam[a][a][a], arat, scrn);

               this->Cmax_meam[a][a][a], this->stheta_meam[a][a], arat, scrn);

      rho0_2nn = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * (arat - 1));

      rho0 = rho0 + Z2 * rho0_2nn * scrn;

    }
@@ -554,10 +599,15 @@ MEAM::get_tavref(double* t11av, double* t21av, double* t31av, double* t12av, dou 
    case FCC:

    case BCC:

    case DIA:

    case DIA3:

    case HCP:

    case B1:

    case DIM:

    case B2:

    case CH4:

    case LIN:

    case ZIG:

    case TRI:

      //     all neighbors are of the opposite type

      *t11av = t12;

      *t21av = t22;
@@ -603,7 +653,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* 
  double a1, a2;

  double s[3];

  lattice_t lat;

  int Zij2nn;

  int Zij,Zij2nn;

  double rhoa01nn, rhoa02nn;

  double rhoa01, rhoa11, rhoa21, rhoa31;

  double rhoa02, rhoa12, rhoa22, rhoa32;
@@ -624,6 +674,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* 


  lat = this->lattce_meam[a][b];



  Zij = get_Zij(lat);



  *rho11 = 0.0;

  *rho21 = 0.0;

  *rho31 = 0.0;
@@ -645,6 +697,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* 
      *rho02 = 6.0 * rhoa01;

      break;

    case DIA:

    case DIA3:

      *rho01 = 4.0 * rhoa02;

      *rho02 = 4.0 * rhoa01;

      *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
@@ -657,7 +710,7 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* 
      *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;

      break;

    case DIM:

      get_shpfcn(DIM, s);

      get_shpfcn(DIM, 0, 0, s);

      *rho01 = rhoa02;

      *rho02 = rhoa01;

      *rho11 = s[0] * rhoa12 * rhoa12;
@@ -692,13 +745,71 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* 
      *rho01 = 8.0 * rhoa02;

      *rho02 = 8.0 * rhoa01;

      break;

    case CH4:

      *rho01 = 4.0 * rhoa02; //in assumption that 'a' represent carbon

      *rho02 = rhoa01;	//in assumption that 'b' represent hydrogen



      get_shpfcn(DIM, 0, 0, s);	//H

      *rho12 = s[0] * rhoa11 * rhoa11;

      *rho22 = s[1] * rhoa21 * rhoa21;

      *rho32 = s[2] * rhoa31 * rhoa31;



      get_shpfcn(CH4, 0, 0, s); //C

      *rho11 = s[0] * rhoa12 * rhoa12;

      *rho21 = s[1] * rhoa22 * rhoa22;

      *rho31 = s[2] * rhoa32 * rhoa32;

      break;

    case LIN:

      *rho01 = rhoa02*Zij;

      *rho02 = rhoa01*Zij;



      get_shpfcn(LIN, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);

      *rho12 = s[0] * rhoa11 * rhoa11;

      *rho22 = s[1] * rhoa21 * rhoa21;

      *rho32 = s[2] * rhoa31 * rhoa31;

      *rho11 = s[0] * rhoa12 * rhoa12;

      *rho21 = s[1] * rhoa22 * rhoa22;

      *rho31 = s[2] * rhoa32 * rhoa32;

      break;

    case ZIG:

      *rho01 = rhoa02*Zij;

      *rho02 = rhoa01*Zij;



      get_shpfcn(ZIG, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);

      *rho12 = s[0] * rhoa11 * rhoa11;

      *rho22 = s[1] * rhoa21 * rhoa21;

      *rho32 = s[2] * rhoa31 * rhoa31;

      *rho11 = s[0] * rhoa12 * rhoa12;

      *rho21 = s[1] * rhoa22 * rhoa22;

      *rho31 = s[2] * rhoa32 * rhoa32;

      break;

    case TRI:

      *rho01 = rhoa02;

      *rho02 = rhoa01*Zij;



      get_shpfcn(TRI, this->stheta_meam[a][b], this->ctheta_meam[a][b], s);

      *rho12 = s[0] * rhoa11 * rhoa11;

      *rho22 = s[1] * rhoa21 * rhoa21;

      *rho32 = s[2] * rhoa31 * rhoa31;

      s[0] = 1.0;

      s[1] = 2.0/3.0;

      s[2] = 1.0 - 0.6*s[0];



      *rho11 = s[0] * rhoa12 * rhoa12;

      *rho21 = s[1] * rhoa22 * rhoa22;

      *rho31 = s[2] * rhoa32 * rhoa32;

      break;





    // default:

    //        call error('Lattice not defined in get_densref.')

  }



  if (this->nn2_meam[a][b] == 1) {



    Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b], arat, scrn);



    Zij2nn = get_Zij2(lat, this->Cmin_meam[a][a][b], this->Cmax_meam[a][a][b],

                      this->stheta_meam[a][b], arat, scrn);



    a1 = arat * r / this->re_meam[a][a] - 1.0;

    a2 = arat * r / this->re_meam[b][b] - 1.0;
@@ -725,7 +836,8 @@ MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, double* 
      *rho01 = *rho01 + Zij2nn * scrn * rhoa01nn;



      //     Assume Zij2nn and arat don't depend on order, but scrn might

      Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a], arat, scrn);

      Zij2nn = get_Zij2(lat, this->Cmin_meam[b][b][a], this->Cmax_meam[b][b][a],

                        this->stheta_meam[a][b], arat, scrn);

      *rho02 = *rho02 + Zij2nn * scrn * rhoa02nn;

    }

  }

CRA Git | Maintained and supported by SUSTech CRA and CCSE