tools/eam_generate/Cu_Mishin1.c

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <stdlib.h>
   4 #include <unistd.h>
   5 #include <stdarg.h>
   6 #include <string.h>
   7
   8 #define rc  5.50679// A
   9 #define h   0.50037// A
  10 #define E1  2.01458*1e2//eV
  11 #define E2  6.59288*1e-3//eV
  12 #define r01 0.83591//A
  13 #define r02 4.46867//A
  14 #define a1  2.97758//1/A
  15 #define a2  1.54927//1/A
  16 #define d   0.86225*1e-2//A
  17 #define rs1 2.24//A
  18 #define rs2 1.8//A
  19 #define rs3 1.2//A
  20 #define S1  4.//eV/A^4
  21 #define S2  40.//eV/A^4
  22 #define S3  1.15*1e3//eV/A^4
  23 #define a   3.80362
  24 #define r03 -2.19885//A
  25 #define r04 -2.61984*1e2//A
  26 #define b1  0.17394//A^-2
  27 #define b2  5.35661*1e2//1/A
  28 #define F0  -2.28235//eV
  29 #define F2  1.35535//eV
  30 #define q1  -1.27775//eV
  31 #define q2  -0.86074//eV
  32 #define q3  1.78804//eV
  33 #define q4  2.97571//eV
  34 #define Q1  0.4
  35 #define Q2  0.3
  36
  37 #define EOK 0
  38 #define ERROR -1
  39
  40
  41 double M (double r, double r0, double alpha) {
  42    return exp (-2.*alpha*(r-r0)) - 2.*exp (-alpha*(r-r0));
  43 }
  44
  45 double psi (double x) {
  46    if (x >= 0.) return 0.;
  47    else return pow (x, 4.) / (1.+pow (x, 4.) );
  48 }
  49
  50 double H (double x) {
  51    if (x >= 0) return 1.;
  52    else return 0.;
  53 }
  54
  55 double V (double r) {
  56    return ( E1*M(r,r01,a1) + E2*M(r,r02,a2) + d ) * psi ( (r-rc)/h ) -
  57       H (rs1-r)*S1*pow (rs1-r, 4.) -
  58       H (rs2-r)*S2*pow (rs2-r, 4.) -
  59       H (rs3-r)*S3*pow (rs3-r, 4.);
  60 }
  61
  62 double rho (double r) {
  63    return ( a*exp (-b1*pow (r-r03, 2.)) + exp (-b2*(r-r04)) ) * psi ( (r-rc)/h );
  64 }
  65
  66 double F (double rho_) {
  67    if (rho_ < 1.)
  68       return F0 + .5*F2*pow (rho_-1., 2.)
  69          + q1*pow (rho_-1.,3.)
  70          + q2*pow (rho_-1.,4.)
  71          + q3*pow (rho_-1.,5.)
  72          + q4*pow (rho_-1.,6.);
  73    else if (rho_ > 1.)
  74       return ( F0 + .5*F2*pow (rho_-1., 2.) + q1*pow (rho_-1., 3.) + Q1*pow (rho_-1., 4.)) /
  75          ( 1. + Q2*pow (rho_-1., 3.) );
  76    else exit (ERROR);
  77
  78 }
  79
  80 int main (void) {
  81    int Nr = 10001;
  82    double rmax = 9.;
  83    double dr = rmax/(double)Nr;
  84    int Nrho = 10001;
  85    double rhomax = rho (0.);
  86    double drho = rhomax/(double)Nrho;
  87
  88    int atomic_number = 1;
  89    double mass = 63.55;
  90    double lattice_constant = 3.615;
  91    char lattice_type[] = "FCC";
  92
  93    int i;
  94
  95    char LAMMPSFilename[] = "Cu_Mishin1.eam";
  96    FILE *LAMMPSFile = fopen (LAMMPSFilename, "w");
  97    if (!LAMMPSFile) exit (ERROR);
  98
  99    // Header for setfl format
 100    fprintf (LAMMPSFile, \
 101          "#-> LAMMPS Potential File in DYNAMO 86 setfl Format <-#\n"\
 102          "# Mishin Cu EAM1 PRB(2001)63:224106\n"\
 103          "# Implemented by G. Ziegenhain (2007) gerolf@ziegenhain.com \n"\
 104          "%d Cu\n"\
 105          "%d %20.20f %d %20.20f %20.20f\n"\
 106          "%d %20.20f %20.20f %s\n",
 107          atomic_number,
 108          Nrho, drho, Nr, dr, rc,
 109          atomic_number, mass, lattice_constant, lattice_type);
 110
 111    // Embedding function
 112    for (i = 0; i < Nrho; i++)
 113       fprintf (LAMMPSFile, "%20.20f\n", F ((double)i*drho));
 114    // Density function
 115    for (i = 0; i < Nr; i++)
 116       fprintf (LAMMPSFile, "%20.20f\n", rho ((double)i*dr));
 117    // Pair potential
 118    for (i = 0; i < Nr; i++)
 119       fprintf (LAMMPSFile, "%20.20f\n", V ((double)i*dr) * (double)i*dr);
 120
 121    fclose (LAMMPSFile);
 122    return (EOK);
 123 }