From e9b11e3673a383f1e52e02b10caf570b9f3446b6 Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Sat, 12 Feb 2022 12:42:54 +0100 Subject: [PATCH] Added Python version of Zhou04_create_v2.f: create_eam.py and eamDatabase.py in tools/eam_database --- tools/eam_database/create_eam.py | 173 +++++++++ tools/eam_database/eamDatabase.py | 610 ++++++++++++++++++++++++++++++ 2 files changed, 783 insertions(+) create mode 100644 tools/eam_database/create_eam.py create mode 100644 tools/eam_database/eamDatabase.py diff --git a/tools/eam_database/create_eam.py b/tools/eam_database/create_eam.py new file mode 100644 index 0000000000..975774466d --- /dev/null +++ b/tools/eam_database/create_eam.py @@ -0,0 +1,173 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +Python version of the code Zhou04_create_v2.f +original author: X. W. Zhou, xzhou@sandia.gov +based on updates by: Lucas Hale lucas.hale@nist.gov +written by: Germain Clavier g.m.g.c.clavier@tue.nl +This script requires the numpy library. +""" + +import sys +import argparse as ap +import numpy as np +from eamDatabase import Database + + +def prof(at, r): + atom = Database[at] + f = np.zeros(r.shape) + f = atom.fe * np.exp(-atom.beta1 * (r[r >= 0.5] / atom.re - 1.0)) + f = f / (1.0 + (r / atom.re - atom.ramda1) ** 20) + return f + + +def pair(at1, at2, r): + atom = Database[at1] + psi1 = atom.A * np.exp(-atom.alpha * (r / atom.re - 1.0)) + psi1 /= 1.0 + (r / atom.re - atom.cai) ** 20 + psi2 = atom.B * np.exp(-atom.beta * (r / atom.re - 1.0)) + psi2 /= 1.0 + (r / atom.re - atom.ramda) ** 20 + if at1 == at2: + psi = psi1 - psi2 + else: + psia = psi1 - psi2 + atom2 = Database[at2] + psi1 = atom2.A * np.exp(-atom2.alpha * (r / atom2.re - 1.0)) + psi1 /= 1.0 + (r / atom2.re - atom2.cai) ** 20 + psi2 = atom2.B * np.exp(-atom2.beta * (r / atom2.re - 1.0)) + psi2 /= 1.0 + (r / atom2.re - atom2.ramda) ** 20 + psib = psi1 - psi2 + prof1 = prof(at1, r) + prof2 = prof(at2, r) + psi = 0.5 * (prof2 / prof1 * psia + prof1 / prof2 * psib) + return psi + + +def embed(at, rho): + atom = Database[at] + Fm33 = np.zeros(rho.shape) + Fm33[rho < atom.rhoe] = atom.Fm3 + Fm33[rho >= atom.rhoe] = atom.Fm4 + emb = np.zeros(rho.shape) + for i, r in enumerate(rho): + if r == 0: + emb[i] = 0 + elif r < atom.rhoin: + dr = r / atom.rhoin - 1 + emb[i] = atom.Fi0 + atom.Fi1 * dr + atom.Fi2 * dr**2 + atom.Fi3 * dr**3 + elif r < atom.rhoout: + dr = r / atom.rhoe - 1 + emb[i] = atom.Fm0 + atom.Fm1 * dr + atom.Fm2 * dr**2 + Fm33[i] * dr**3 + else: + dr = r / atom.rhos + emb[i] = atom.Fn * (1.0 - atom.fnn * np.log(dr)) * dr**atom.fnn + return emb + + +def write_file(attypes, filename, Fr, rhor, z2r, nrho, drho, nr, dr, rc): + struc = "fcc" + with open(filename, "w") as f: + f.write("DATE: 2022-02-11") + f.write(" CONTRIBUTORS: Xiaowang Zhou xzhou@sandia.gov and") + f.write(" Lucas Hale lucas.hale@nist.gov") + f.write(" Germain Clavier g.m.g.c.clavier@tue.nl/germain.clavier@gmail.com\n") + f.write(" CITATION: X. W. Zhou, R. A. Johnson, H. N. G. Wadley, Phys. Rev. B, 69, 144113(2004)\n") + f.write("Generated from create_eam.py, Python version of Zhou04_create_v2.f\n") + f.write("{:<5d} {:<24}\n".format(len(attypes), " ".join(attypes))) + f.write( + "{:<5d} {:<24.16e} {:<5d} {:<24.16e} {:<24.16e}\n".format( + nrho, drho, nr, dr, rc + ) + ) + for at in attypes: + atom = Database[at] + f.write( + "{:>5d} {:>15.5f} {:>15.5f} {:>8}\n".format( + atom.ielement, atom.amass, atom.blat, struc + ) + ) + for i, fr in enumerate(Fr[at]): + f.write(" {:>24.16E}".format(fr)) + if not (i + 1) % 5: + f.write("\n") + for i, rho in enumerate(rhor[at]): + f.write(" {:>24.16E}".format(rho)) + if not (i + 1) % 5: + f.write("\n") + for n1 in range(len(attypes)): + for n2 in range(n1 + 1): + for i, z in enumerate(z2r[n1, n2]): + f.write(" {:>24.16E}".format(z)) + if not (i + 1) % 5: + f.write("\n") + + return + + +def main(): + + parser = ap.ArgumentParser(description="Script to make EAM alloy file inputs.") + parser.add_argument("-n", "--names", dest="names", nargs="+", help="Atom names.") + parser.add_argument("-nr", dest="nr", type=int, default=2000, help="Number of point in r space [default 2000].") + parser.add_argument("-nrho", dest="nrho", type=int, default=2000, help="Number of point in rho space [default 2000].") + args = parser.parse_args() + atnames = args.names + nr = args.nr + nrho = args.nrho + + for n in atnames: + try: + Database[n] + except KeyError: + output = "Atom {} not found.\n".format(n) + valid = "Valid inputs are: {}".format(" ".join(Database.keys())) + sys.exit("".join([output, valid])) + + ntypes = len(atnames) + outfilename = "".join([*atnames, ".set"]) + rhor = {} + Fr = {} + + alatmax = max([Database[at].blat for at in atnames]) + rhoemax = max([Database[at].rhoe for at in atnames]) + rc = np.sqrt(10.0) / 2.0 * alatmax + rst = 0.5 + r = np.linspace(0.0, rc, num=nr, dtype=np.double) + dr = r[1] - r[0] + r[r < rst] = rst + z2r = np.zeros([ntypes, ntypes, nr]) + fmax = -np.inf + rhomax = -np.inf + + for i, n1 in enumerate(atnames): + for j, n2 in enumerate(atnames): + if j > i: + continue + elif i == j: + rhor[n1] = prof(n1, r) + rhomax = np.max(rhor[n1]) if rhomax < np.max(rhor[n1]) else rhomax + z2r[i, j, :] = r * pair(n1, n2, r) + else: + z2r[i, j, :] = r * pair(n1, n2, r) + z2r = np.where(z2r, z2r, z2r.transpose((1, 0, 2))) + if rhomax < 2.0 * rhoemax: + rhomax = 2.0 * rhoemax + if rhomax < 100.0: + rhomax = 100.0 + rho = np.linspace(0.0, rhomax, num=nrho, dtype=np.double) + drho = rho[1] - rho[0] + for i, n1 in enumerate(atnames): + Fr[n1] = embed(n1, rho) + + write_file(atnames, outfilename, Fr, rhor, z2r, nrho, drho, nr, dr, rc) + + return + + +if __name__ == "__main__": + try: + main() + except KeyboardInterrupt: + raise SystemExit("User interruption.") diff --git a/tools/eam_database/eamDatabase.py b/tools/eam_database/eamDatabase.py new file mode 100644 index 0000000000..26d22d2f3a --- /dev/null +++ b/tools/eam_database/eamDatabase.py @@ -0,0 +1,610 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +Python version of the code Zhou04_create_v2.f +original author: X. W. Zhou, xzhou@sandia.gov +based on updates by: Lucas Hale lucas.hale@nist.gov +written by: Germain Clavier g.m.g.c.clavier@tue.nl + +This file contains atom attributes for the EAM model and alloy combination. The +original file is EAM_code. It is designed to be used with create_eam.py script. +To add new contribution, just add new AtType instances. +""" + +import math + +class AtType: + def __init__( + self, + name, + re, + fe, + rhoe, + rhos, + alpha, + beta, + A, + B, + cai, + ramda, + Fi0, + Fi1, + Fi2, + Fi3, + Fm0, + Fm1, + Fm2, + Fm3, + fnn, + Fn, + ielement, + amass, + Fm4, + beta1, + ramda1, + rhol, + rhoh, + ): + self.name = name + self.re = re + self.fe = fe + self.rhoe = rhoe + self.rhos = rhos + self.alpha = alpha + self.beta = beta + self.A = A + self.B = B + self.cai = cai + self.ramda = ramda + self.Fi0 = Fi0 + self.Fi1 = Fi1 + self.Fi2 = Fi2 + self.Fi3 = Fi3 + self.Fm0 = Fm0 + self.Fm1 = Fm1 + self.Fm2 = Fm2 + self.Fm3 = Fm3 + self.fnn = fnn + self.Fn = Fn + self.ielement = ielement + self.amass = amass + self.Fm4 = Fm4 + self.beta1 = beta1 + self.ramda1 = ramda1 + self.rhol = rhol + self.rhoh = rhoh + self.blat = math.sqrt(2.0) * self.re + self.rhoin = self.rhol * self.rhoe + self.rhoout = self.rhoh * self.rhoe + + def __repr__(self): + output = """{}: + re = {}; fe = {} + rhoe = {}; rhos = {}; + alpha = {}; beta = {}; + A = {}; B = {} + cai = {}; ramda = {} + Fi0 = {}; Fi1 = {}; Fi2 = {}; Fi3 = {} + Fm0 = {}; Fm1 = {}; Fm2 = {}; Fm3 = {}; Fm4 = {} + fnn = {}; Fn = {} + ielement = {}; amass = {} + beta1 = {}; ramda1 = {} + rhol = {}; rhoh = {} + blat = {}; rhoin = {}; rhoout = {}""" + return output.format( + self.name, + self.re, + self.fe, + self.rhoe, + self.rhos, + self.alpha, + self.beta, + self.A, + self.B, + self.cai, + self.ramda, + self.Fi0, + self.Fi1, + self.Fi2, + self.Fi3, + self.Fm0, + self.Fm1, + self.Fm2, + self.Fm3, + self.Fm4, + self.fnn, + self.Fn, + self.ielement, + self.amass, + self.beta1, + self.ramda1, + self.rhol, + self.rhoh, + self.blat, + self.rhoin, + self.rhoout, + ) + + +Database = {} +Database["Cu"] = AtType( + "Cu", # Name + 2.556162, # re + 1.554485, # fe + 21.175871, # rhoe + 21.175395, # rhos + 8.127620, # alpha + 4.334731, # beta + 0.396620, # A + 0.548085, # B + 0.308782, # cai + 0.756515, # ramda + -2.170269, # Fi0 + -0.263788, # Fi1 + 1.088878, # Fi2 + -0.817603, # Fi3 + -2.19, # Fm0 + 0.00, # Fm1 + 0.561830, # Fm2 + -2.100595, # Fm3 + 0.310490, # fnn + -2.186568, # Fn + 29, # ielement + 63.546, # amass + -2.100595, # Fm4 + 4.334731, # beta1 + 0.756515, # ramda1 + 0.85, # rhol + 1.15, # rhoh +) +Database["Ag"] = AtType( + "Ag", + 2.891814, + 1.106232, + 14.604100, + 14.604144, + 9.132010, + 4.870405, + 0.277758, + 0.419611, + 0.339710, + 0.750758, + -1.729364, + -0.255882, + 0.912050, + -0.561432, + -1.75, + 0.00, + 0.744561, + -1.150650, + 0.783924, + -1.748423, + 47, + 107.8682, + -1.150650, + 4.870405, + 0.750758, + 0.85, + 1.15, +) +Database["Au"] = AtType( + "Au", + 2.885034, + 1.529021, + 19.991632, + 19.991509, + 9.516052, + 5.075228, + 0.229762, + 0.356666, + 0.356570, + 0.748798, + -2.937772, + -0.500288, + 1.601954, + -0.835530, + -2.98, + 0.00, + 1.706587, + -1.134778, + 1.021095, + -2.978815, + 79, + 196.96654, + -1.134778, + 5.075228, + 0.748798, + 0.85, + 1.15, +) +Database["Ni"] = AtType( + "Ni", + 2.488746, + 2.007018, + 27.562015, + 27.562031, + 8.383453, + 4.471175, + 0.429046, + 0.633531, + 0.443599, + 0.820658, + -2.693513, + -0.076445, + 0.241442, + -2.375626, + -2.70, + 0.00, + 0.265390, + -0.152856, + 0.445470, + -2.7, + 28, + 58.6934, + -0.152856, + 4.471175, + 0.820658, + 0.85, + 1.15, +) +Database["Pd"] = AtType( + "Pd", + 2.750897, + 1.595417, + 21.335246, + 21.940073, + 8.697397, + 4.638612, + 0.406763, + 0.598880, + 0.397263, + 0.754799, + -2.321006, + -0.473983, + 1.615343, + -0.231681, + -2.36, + 0.00, + 1.481742, + -1.675615, + 1.130000, + -2.352753, + 46, + 106.42, + -1.675615, + 4.638612, + 0.754799, + 0.85, + 1.15, +) +Database["Pt"] = AtType( + "Pt", + 2.771916, + 2.336509, + 33.367564, + 35.205357, + 7.105782, + 3.789750, + 0.556398, + 0.696037, + 0.385255, + 0.770510, + -1.455568, + -2.149952, + 0.528491, + 1.222875, + -4.17, + 0.00, + 3.010561, + -2.420128, + 1.450000, + -4.145597, + 78, + 195.08, + -2.420128, + 3.789750, + 0.770510, + 0.25, + 1.15, +) +Database["Al"] = AtType( + "Al", + 2.863924, + 1.403115, + 20.418205, + 23.195740, + 6.613165, + 3.527021, + 0.314873, + 0.365551, + 0.379846, + 0.759692, + -2.807602, + -0.301435, + 1.258562, + -1.247604, + -2.83, + 0.00, + 0.622245, + -2.488244, + 0.785902, + -2.824528, + 13, + 26.981539, + -2.488244, + 3.527021, + 0.759692, + 0.85, + 1.15, +) +Database["Pb"] = AtType( + "Pb", + 3.499723, + 0.647872, + 8.450154, + 8.450063, + 9.121799, + 5.212457, + 0.161219, + 0.236884, + 0.250805, + 0.764955, + -1.422370, + -0.210107, + 0.682886, + -0.529378, + -1.44, + 0.00, + 0.702726, + -0.538766, + 0.935380, + -1.439436, + 82, + 207.2, + -0.538766, + 5.212457, + 0.764955, + 0.85, + 1.15, +) +Database["Fe"] = AtType( + "Fe", + 2.481987, + 1.885957, + 20.041463, + 20.041463, + 9.818270, + 5.236411, + 0.392811, + 0.646243, + 0.170306, + 0.340613, + -2.534992, + -0.059605, + 0.193065, + -2.282322, + -2.54, + 0.00, + 0.200269, + -0.148770, + 0.391750, + -2.539945, + 26, + 55.847, + -0.148770, + 5.236411, + 0.340613, + 0.85, + 1.15, +) +Database["Mo"] = AtType( + "Mo", + 2.728100, + 2.723710, + 29.354065, + 29.354065, + 8.393531, + 4.476550, + 0.708787, + 1.120373, + 0.137640, + 0.275280, + -3.692913, + -0.178812, + 0.380450, + -3.133650, + -3.71, + 0.00, + 0.875874, + 0.776222, + 0.790879, + -3.712093, + 42, + 95.94, + 0.776222, + 4.476550, + 0.275280, + 0.85, + 1.15, +) +Database["Ta"] = AtType( + "Ta", + 2.860082, + 3.086341, + 33.787168, + 33.787168, + 8.489528, + 4.527748, + 0.611679, + 1.032101, + 0.176977, + 0.353954, + -5.103845, + -0.405524, + 1.112997, + -3.585325, + -5.14, + 0.00, + 1.640098, + 0.221375, + 0.848843, + -5.141526, + 73, + 180.9479, + 0.221375, + 4.527748, + 0.353954, + 0.85, + 1.15, +) +Database["W"] = AtType( + "W", + 2.740840, + 3.487340, + 37.234847, + 37.234847, + 8.900114, + 4.746728, + 0.882435, + 1.394592, + 0.139209, + 0.278417, + -4.946281, + -0.148818, + 0.365057, + -4.432406, + -4.96, + 0.00, + 0.661935, + 0.348147, + 0.582714, + -4.961306, + 74, + 183.84, + 0.348147, + 4.746728, + 0.278417, + 0.85, + 1.15, +) +Database["Mg"] = AtType( + "Mg", + 3.196291, + 0.544323, + 7.132600, + 7.132600, + 10.228708, + 5.455311, + 0.137518, + 0.225930, + 0.5, + 1.0, + -0.896473, + -0.044291, + 0.162232, + -0.689950, + -0.90, + 0.00, + 0.122838, + -0.226010, + 0.431425, + -0.899702, + 12, + 24.305, + -0.226010, + 5.455311, + 1.0, + 0.85, + 1.15, +) +Database["Co"] = AtType( + "Co", + 2.505979, + 1.975299, + 27.206789, + 27.206789, + 8.679625, + 4.629134, + 0.421378, + 0.640107, + 0.5, + 1.0, + -2.541799, + -0.219415, + 0.733381, + -1.589003, + -2.56, + 0.00, + 0.705845, + -0.687140, + 0.694608, + -2.559307, + 27, + 58.9332, + -0.687140, + 4.629134, + 1.0, + 0.85, + 1.15, +) +Database["Ti"] = AtType( + "Ti", + 2.933872, + 1.863200, + 25.565138, + 25.565138, + 8.775431, + 4.680230, + 0.373601, + 0.570968, + 0.5, + 1.0, + -3.203773, + -0.198262, + 0.683779, + -2.321732, + -3.22, + 0.00, + 0.608587, + -0.750710, + 0.558572, + -3.219176, + 22, + 47.88, + -0.750710, + 4.680230, + 1.0, + 0.85, + 1.15, +) +Database["Zr"] = AtType( + "Zr", + 3.199978, + 2.230909, + 30.879991, + 30.879991, + 8.559190, + 4.564902, + 0.424667, + 0.640054, + 0.5, + 1.0, + -4.485793, + -0.293129, + 0.990148, + -3.202516, + -4.51, + 0.00, + 0.928602, + -0.981870, + 0.597133, + -4.509025, + 40, + 91.224, + -0.981870, + 4.564902, + 1.0, + 0.85, + 1.15, +)