Polarizable Potentials For Metals: The Density Readjusting Embedded Atom Method (DR-EAM)

TL;DR

This paper introduces DR-EAM, an extension of the embedded atom method that models electronic polarization in metals by dynamically adjusting valence electron densities, enabling more accurate simulations of metallic interfaces and responses to external fields.

Contribution

The paper presents DR-EAM, a novel modification of EAM that incorporates fluctuating valence densities to simulate electronic polarization effects in metals.

Findings

Successfully models polarization and image charge effects.

Predicts surface charging and shielding in electric fields.

Shows charge transfer in alloys affecting unit cell geometries.

Abstract

In simulations of metallic interfaces, a critical aspect of metallic behavior is missing from the some of the most widely used classical molecular dynamics force fields. We present a modification of the embedded atom method (EAM) which allows for electronic polarization of the metal by treating the valence density around each atom as a fluctuating dynamical quantity. The densities are represented by a set of additional fluctuating variables (and their conjugate momenta) which are propagated along with the nuclear coordinates. This ``density readjusting EAM'' (DR-EAM) preserves nearly all of the useful qualities of traditional EAM, including bulk elastic properties and surface energies. However, it also allows valence electron density to migrate through the metal in response to external perturbations. We show that DR-EAM can successfully model polarization in response to external charges, capturing the image charge effect in atomistic simulations. DR-EAM also captures some of the behavior of metals in the presence of uniform electric fields, predicting surface charging and shielding internal to the metal. We further show that it predicts charge transfer between the constituent atoms in alloys, leading to novel predictions about unit cell geometries in layered L$1_0$ structures.

Tables

Table 1. Table 1: DR-EAM uses a count of the effective number of charge carriers per atom ( N eff subscript 𝑁 eff N_{\mathrm{eff}} ) as a valency count in Eq. ( 2 ) and the chemical hardness, K Slater subscript 𝐾 Slater K_{\mathrm{Slater}} , defined in Eq. ( 18 ), to parametrize the harmonic term ( a 2 ) subscript 𝑎 2 (a_{2}) in the self potential. NB: R H subscript 𝑅 𝐻 R_{H} is typically reported in units of cm 3 ​ C − 1 superscript cm 3 superscript C 1 \mathrm{cm}^{3}\mathrm{~{}C}^{-1}

	$\frac{m}{\rho \left|R_H\right|}\cdot \frac{1}{N_{A}e}$	Neff	K
			$K_{\mathrm{Slater}}$	$K_{\mathrm{empirical}}$
Cu	1.42	0.57	7.63	7.87
Ag	1.21	0.48	7.58	8.10
Au	1.48	0.59	7.92	9.38
Ni	1.12	0.45	8.08	7.78
Pd	1.21	0.48	7.59	9.13
Pt	4.10	1.64	6.15	12.34
Al	3.05	1.22	5.54	11.01
Pb1	4.57	1.83	6.70	12.50
Fe	3.00	1.20	9.49	9.88
Mo	0.77	0.31	7.38	5.62
Ta	1.22	0.49	7.12	7.51
W	0.86	0.34	7.79	6.33
Mg	1.68	0.67	7.68	8.31
Co	0.82	0.33	8.31	6.86
Ti2	1.15	0.46	7.18	7.15
Zr	1.12	0.45	6.24	5.85

Table 2. Table 2: Coefficients for the self potential [Eq. ( 5 )] used in DR-EAM. All parameters have units that give the self potential in eV.

Element	$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$a_6$
Cu	10.75	7.63	11.12	65.63	-72.22	21.88
Ag	10.90	7.58	14.04	70.03	-82.35	25.21
Au	13.89	7.92	13.19	68.17	-88.02	28.85
Ni	10.80	8.08	5.62	58.45	-58.13	16.77
Pd	12.22	7.59	2.36	84.42	-89.92	26.68
Pt	12.61	6.15	14.07	67.32	-86.75	28.28
Al	9.33	5.54	19.18	81.65	-136.26	57.17
Pb	12.86	6.70	-16.03	79.12	-59.70	13.84
Fe	10.41	9.49	1.14	64.45	-69.63	22.15
Mo	12.30	7.38	-5.77	70.98	-68.99	19.90
Ta	8.80	7.12	13.48	67.67	-92.80	31.22
W	13.01	7.79	-4.39	74.83	-78.12	23.51
Mg	7.86	7.68	52.67	97.79	-338.60	211.34
Co	10.66	8.31	0.67	51.63	-41.62	9.92
Ti	8.99	7.18	-1.57	58.59	-56.70	16.41
Zr	7.96	6.42	5.11	42.34	-44.13	12.65

Table 3. Table 3: Comparison of Mechanical Properties for Bulk Metals. For each element, the top row includes the cohesive energy per atom ( E coh subscript 𝐸 coh E_{\mathrm{coh}} ), vacancy formation energy ( E v f superscript subscript 𝐸 𝑣 𝑓 E_{v}^{f} ), bulk modulus ( B 𝐵 B ), shear modulus ( G 𝐺 G ), and Poisson’s ratio ( σ 𝜎 \sigma ). These values (identical for DR-EAM and the unmodified EAM model within numerical precision) are compared with a second line from DFT calculations (when available). A third line with experimental values is also presented. McKee  et al. ( 1972 ); Medasani  et al. ( 2015 ); Ong  et al. ( 2015 ); Kim and Buyers ( 1978 ); Tzanetakis  et al. ( 1976 ); Janthon  et al. ( 2014 ); Kittel ( 2005 ); Simmons and Wang ( 1971 )

	Ecoh (eV)	E${}_{}{}^f{}_v{}^{}$ (eV)	B (GPa)	G (GPa)	$\sigma$			Ecoh (eV)	E${}_{}{}^f{}_v{}^{}$ (eV)	B (GPa)	G (GPa)	$\sigma$
Cu	3.65	1.42	125.74	64.48	0.28		Al	3.74	0.74	52.76	17.52	0.28
	4.09	-	145	57	0.34			3.74	-	83	25	0.37
	3.43	1.28	140	59.3	0.31			3.39	0.68	76	29.0	0.35
Ag	2.89	1.12	105.77	41.29	0.32		Pb	2.01	0.60	47.65	13.10	0.34
	2.83	-	88	28	0.38			3.71	-	37	14	0.35
	2.95	1.11	100	37.5	0.34			2.03	0.50	46	15.7	0.37
Au	3.94	0.92	169.65	38.39	0.39		Fe	4.34	1.65	117.78	58.08	0.39
	3.27	-	137	19	0.45			8.45	-	182	78	0.32
	3.81	0.89	180	33.7	0.41			4.28	1.40	170	94.1	0.27
Ni	4.51	1.78	146.23	74.41	0.28		Mo	7.07	2.76	179.48	60.70	0.31
	5.78	-	198	102	0.29			10.86	-	262	127	0.30
	4.44	1.79	180	101.1	0.27			6.82	3.0	230	130.4	0.28
Pd	4.00	1.56	190.36	62.04	0.35		Ta	8.37	2.80	170.83	39.06	0.28
	5.18	-	160	50	0.38			11.85	-	194	63	0.35
	3.89	1.85	180	54.3	0.37			8.10	-	200	74.0	0.33
Pt	5.98	1.34	252.62	56.49	0.37		W	8.95	3.62	191.97	82.61	0.31
	6.05	-	247	49	0.41			12.95	-	304	147	0.29
	5.84	1.32	230	67.3	0.39			8.90	4.0	310	163.4	0.27
Mg	1.54	0.65	34.83	17.19	0.28		Co	4.40	1.83	190.54	79.99	0.31
	1.59	-	37	18	0.29			7.11	-	212	106	0.29
	1.51	0.79	35.4	19.4	0.27			4.39	1.34	191.4	91.1	0.29
Ti	4.54	0.77	214.06	109.5	0.28		Zr	6.33	1.95	89.28	31.92	0.34
	7.88	-	113	47	0.32			8.54	-	94	35	0.34
	4.85	1.55	110	50.5	0.30			6.29	-	83.3	42.5	0.31

Table 4. Table 4: Comparison of surface energies (all reported in J m -2 ), and average atomic charge on the surface layers ⟨ q ⟩ delimited-⟨⟩ 𝑞 \langle q\rangle (in e). Experimental surface energies are averaged over many exposed facets. Vitos  et al. ( 1998 ) In the first data column, the DR-EAM surface energy is provided. When nonpolarizable EAM predicts a different surface energy, this value is given in parentheses. Values of ⟨ q ⟩ delimited-⟨⟩ 𝑞 \langle q\rangle smaller than 10 − 4 superscript 10 4 10^{-4} e are indicated with an asterisk.

	Surface	DR-EAM (EAM)	DFT Ong et al. (2015)	Experiment Vitos et al. (1998)	$⟨q⟩$			Surface	DR-EAM (EAM)	DFT Ong et al. (2015)	Experiment Vitos et al. (1998)	$⟨q⟩$
	(111)	1.599	1.31		*			(111)	0.352	0.25		*
Cu	(100)	1.655	1.47	1.790	*		Pb	(100)	0.395	0.28	0.593	*
	(110)	1.852	1.56		*			(110)	0.428	0.33		*
	(111)	0.947	0.77		*			(111)	1.988	2.73		*
Ag	(100)	1.025	0.81	1.246	*		Fe	(100)	1.441	2.50	2.417	*
	(110)	1.119 (1.115)	0.87		0.001			(110)	1.350	2.45		*
	(111)	0.923	0.74		*			(111)	2.908	2.96		*
Au	(100)	1.043	0.86	1.506	*		Mo	(100)	2.531	3.18	2.907	*
	(110)	1.124	0.83		*			(110)	2.351	2.80		*
	(111)	1.949	1.92		*			(111)	2.463	2.70		*
Ni	(100)	2.041	2.21	2.380	*		Ta	(100)	2.342	2.47	2.902	*
	(110)	2.231	2.29		*			(110)	1.954	2.34		*
	(111)	1.601	1.34		*			(111)	3.028	3.47		*
Pd	(100)	1.728	1.53	2.003	*		W	(100)	2.967	3.95	3.265	*
	(110)	1.890	1.57		*			(110)	2.826	3.23		*
	(111)	1.978	1.48		*			(111)	0.914	0.80		*
Pt	(100)	2.384	1.84	2.489	*		Al	(100)	0.891 (0.883)	0.92	1.143	-0.002
	(110)	2.150	1.68		*			(110)	1.011 (0.986)	0.98		0.001
Mg	(0001)	0.377	0.54	0.760	*		Zr	(0001)	1.278	1.97	1.909	*
Co	(0001)	2.028	2.11	2.522	*		Ti	(0001)	5.577	1.61	1.989	*

Table 5. Table 5: The lattice constants ( a 𝑎 a ), bulk modulus ( B ) 𝐵 (B) , shear modulus ( G ) 𝐺 (G) , Poisson’s ratio ( σ ) 𝜎 (\sigma) , energy per atom ( E a ) subscript 𝐸 𝑎 (E_{a}) and charge on atom B at surface ⟨ q B ⟩ delimited-⟨⟩ subscript 𝑞 B \langle q_{\mathrm{B}}\rangle for L 1 2 subscript 1 2 1_{2} ( \ce AB3) alloys. Note: ⟨ q A ⟩ = − 3 ​ ⟨ q B ⟩ delimited-⟨⟩ subscript 𝑞 A 3 delimited-⟨⟩ subscript 𝑞 B \langle q_{\mathrm{A}}\rangle=-3\langle q_{\mathrm{B}}\rangle for these alloys. For each alloy, the top row is the value calculated using DR-EAM, the middle row is from nonpolarizable EAM, and the lower values are derived from plane-wave DFT calculations. Ong  et al. ( 2015 )

Alloy	$a$ ($\AA$)	$B$ (GPa)	$G$ (GPa)	$\sigma$	$E_a$ (eV)	$⟨q_B⟩$ (e)		Alloy	$a$ ($\AA$)	$B$ (GPa)	$G$ (GPa)	$\sigma$	$E_a$ (eV)	$⟨q_B⟩$ (e)
NiAg3	3.99	111.70	39.01	0.34	-3.20	-0.002		PbPd3	4.17	134.67	36.20	0.37	-3.22	0.014
	3.99	111.70	39.01	0.34	-3.20	-			4.17	134.66	36.19	0.37	-3.21	-
	4.00	64	27	0.33	-3.37	-			4.12	132	47	0.35	-5.11	-
CuAu3	3.99	147.71	37.95	0.37	-4.09	-0.047		AgAu3	4.09	146.81	35.77	0.38	-3.78	-0.043
	3.99	147.57	37.88	0.38	-3.98	-			4.09	146.70	35.73	0.38	-3.72	-
	4.05	143	31	0.41	-3.49	-			4.17	126	23	0.44	-3.19	-
AuCu3	3.76	122.06	50.17	0.31	-3.94	0.057		NiAu3	3.98	153.92	37.58	0.38	-4.19	-0.049
	3.76	122.04	50.01	0.31	-3.86	-			3.98	153.72	37.50	0.38	-4.13	-
	3.77	142	52	0.35	-3.92	-			4.02	147	22	0.43	-3.80	-
PdCu3	3.73	132.73	42.50	0.34	-3.61	0.030		PtCu3	3.64	157.64	59.55	0.33	-4.46	0.051
	3.72	132.73	54.65	0.31	-3.79	-			3.64	157.22	59.42	0.33	-4.42	-
	3.71	146	53	0.36	-4.47	-			3.72	159	59	0.35	-4.72	-
PtFe3	3.55	236.48	94.07	0.32	-5.17	0.054		FePd3	3.82	175.88	57.43	0.35	-4.16	-0.031
	3.56	235.53	93.52	0.32	-5.11	-			3.82	175.78	57.42	0.35	-4.13	-
	3.75	180	50	0.40	-7.94	-			3.89	174	76	0.32	-6.11	-
AlNi3	3.66	111.60	48.83	0.30	-4.37	0.014
	3.65	111.56	48.74	0.30	-4.36	-
	3.56	117	96	0.29	-5.70	-

Table 6. Table 6: The lattice constants ( a 𝑎 a and c 𝑐 c ), bulk modulus ( B ) 𝐵 (B) , shear modulus ( G ) 𝐺 (G) , Poisson’s ratio ( σ ) 𝜎 (\sigma) , energy per atom ( E a ) subscript 𝐸 𝑎 (E_{a}) and charge on B atoms ⟨ q B ⟩ delimited-⟨⟩ subscript 𝑞 B \langle q_{\mathrm{B}}\rangle for L 1 0 subscript 1 0 1_{0} ( \ce AB) alloys. Note: ⟨ q A ⟩ = − ⟨ q B ⟩ delimited-⟨⟩ subscript 𝑞 A delimited-⟨⟩ subscript 𝑞 B \langle q_{\mathrm{A}}\rangle=-\langle q_{\mathrm{B}}\rangle for these alloys. For each alloy, the top row is the value calculated using DR-EAM, the second row is from nonpolarizable EAM, the third row are derived from plane-wave DFT calculations in Ref. Ong  et al. , 2015 and, the bottom row are the experimental values in Ref. Pearson, 1958 (when available).

Alloy	$c$ ($\AA$)	$a$ ($\AA$)	$c/a$	$B$ (GPa)	$G$ (GPa)	$\sigma$	$E_a$ (eV)	$⟨q_B⟩$ (e)		Alloy	$c$ ($\AA$)	$a$ ($\AA$)	$c/a$	$B$ (GPa)	$G$ (GPa)	$\sigma$	$E_a$ (eV)	$⟨q_B⟩$ (e)
AgTi	4.26	4.09	1.04	117.33	40.74	0.34	-3.95	0.090		FePd	3.62	3.77	0.95	179.37	67.84	0.33	-4.27	-0.070
	4.29	4.08	1.05	117.22	40.58	0.34	-3.70	-			3.62	3.77	0.96	178.99	67.65	0.33	-4.23	-
	4.17	4.085	1.02	111	37	0.36	-5.43	-			3.78	3.84	0.98	163	82	0.3	-6.88	-
	4.09	4.09	0.993	-	-	-	-	-			3.72	3.85	0.96	-	-	-	-	-
AlMg	4.29	4.24	1.01	73.19	35.30	0.29	-2.61	0.050		FePt	3.49	3.77	0.92	236.42	75.57	0.35	-5.66	-0.104
	4.29	4.24	1.01	73.47	35.50	0.29	-2.58	-			3.53	3.75	0.94	234.42	75.28	0.35	-5.60	-
	4.23	3.02	1.40	58	33	-	-2.68	-			3.72	3.76	1.01	201	98	0.3	-7.5	-
	-	-	-	-	-	-	-	-			3.71	3.83	0.96	-	-	-	-	-
AlTi	4.27	4.10	1.04	81.09	29.98	0.33	-4.29	0.020		MgNi	4.14	3.99	1.03	101.53	38.89	0.33	-2.54	-0.079
	4.27	4.10	1.04	81.17	29.91	0.33	-4.28	-			4.16	3.99	1.04	101.12	38.63	0.33	-2.47	-
	4.08	3.99	1.02	115	76	0.25	-6.22	-			3.16	2.99	1.05	150	55	0.35	-3.84	-
	4.06	3.99	1.02	-	-	-	-	-			-	-	-	-	-	-	-	-
AuCu	3.52	4.04	0.87	136.21	43.23	0.35	-4.06	0.119		MgZr	4.85	4.22	1.14	105.69	41.26	0.32	-4.03	-0.004
	3.66	3.97	0.92	135.45	42.84	0.35	-3.95	-			4.85	4.22	1.14	105.69	41.26	0.32	-4.03	-
	3.66	4.05	0.90	135	37	0.4	-3.73	-			4.46	3.16	1.41	-	-	-	-4.99	-
	3.66	3.92	0.92	-	-	-	-	-			-	-	-	-	-	-	-	-
CoPt	3.62	3.77	0.95	213.43	61.59	0.36	-5.52	-0.103		NiPt	3.59	3.77	0.95	202.64	59.68	0.36	-5.56	-0.094
	3.64	3.77	0.96	212.85	64.82	0.36	-5.46	-			3.43	3.86	0.88	200.37	58.23	0.36	-5.51	-
	3.72	3.82	0.97	216	114	0.29	-6.67	-			3.63	3.85	0.94	214	95	0.31	-6.02	-
	3.70	3.81	0.97	-	-	-	-	-			3.58	3.81	0.93	-	-	-	-	-
FeNi	3.51	3.51	1.00	162.43	83.79	0.27	-4.51	-0.017		PdTi	3.97	4.09	0.97	143.42	45.36	0.35	-4.65	0.119
	3.53	3.51	1.00	162.39	83.76	0.27	-4.51	-			4.03	4.06	0.99	142.77	44.99	0.35	-4.54	-
	3.58	3.55	1.00	187	103	0.29	-7.19	-			3.87	2.84	1.46	-	-	-	-7.07	-
	3.61	3.56	0.98	-	-	-	-	-			-	-	-	-	-	-	-	-

Table 7. Table 7: Parameters of the effective image charge potential [Eq. ( 30 )], s 𝑠 s and δ 𝛿 \delta (Å), for three different metals (Cu, Au, Pt) with exposed (100) and (111) facets. A point charge representing a chloride ion is brought close to the surface, and the effective potential is computed. For a (flat) infinite planar conductor, we would expect s → 2 → 𝑠 2 s\rightarrow 2 and δ → 0 → 𝛿 0 \delta\rightarrow 0 . This data suggests that the atomistic polarizable model in DR-EAM is capturing most of the classical image effect.

Metals	Facet	Atop		Hollow		Bridge
		$s$	$\delta$	$s$	$\delta$	$s$	$\delta$
Copper	100	1.86	0.84	1.76	1.08	1.80	0.98
	111	1.84	0.72	1.80	0.80	1.74	0.98
Gold	100	1.78	1.25	1.78	1.24	1.70	1.49
	111	1.84	1.07	1.79	1.22	1.79	1.22
Platinum	100	2.19	0.36	1.51	1.69	1.64	1.31
	111	2.18	0.34	1.59	1.36	1.62	1.26

Equations

$$V\left(\{\mathbf{r}\},\{q\}\right)=\sum _i{F}_i\left[{\bar{\rho }}_i\right]+\frac{1}{2}\sum _i\sum _{j\ne i}{\varphi }_{ij}(r_{ij},q_i,q_j)+\frac{1}{2}\sum _i\sum _{j\ne i}{q}_i{q}_{j}J\left(r_{ij}\right)+\sum _i{V}_{\mathrm{self}}\left(q_i\right).$$

$${\bar{\rho }}_i=\sum _{j\ne i}\left(1-\frac{q_j}{N_j}\right)f_j\left(r_{ij}\right),$$

$${\varphi }_{ij}\left(r,q_i,q_j\right)=\frac{1}{2}\left[\frac{f_j(r)\left(1-\frac{q_j}{N_j}\right)}{f_i(r)\left(1-\frac{q_i}{N_i}\right)}{\varphi }_{ii}(r)+\frac{f_i(r)\left(1-\frac{q_i}{N_i}\right)}{f_j(r)\left(1-\frac{q_j}{N_j}\right)}{\varphi }_{jj}(r)\right],$$

$$\begin{split}J(r)=\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r\right)}{r}-\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}}\\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_{\mathrm{c}}\right)}{R_{\mathrm{c}}^{2}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^{2}R_{\mathrm{c}}^{2}\right)}{R_{\mathrm{c}}}\right)\left(r-R_{\mathrm{c}}\right)\ \right]\quad r\leqslant R_{\textrm{c}},\end{split}$$

$$V_{\mathrm{self}}(q)=\sum _{n=1}^6{a}_n{q}^n.$$

$${\chi }_i={\left(\frac{\partial V_{\mathrm{self}}}{\partial q}\right)}_{q=0}=\frac{IP-EA}{2},$$

$$J_{ii}={\left(\frac{{\partial }^2{V}_{\mathrm{self}}}{\partial q^2}\right)}_{q=0}=IP+EA,$$

$${\chi }_{\mathrm{Pauling}}\approx 0.5106{\chi }_{\mathrm{Mulliken}}-0.2963,$$

$$a_1=2.5{\chi }_{\mathrm{Mulliken}}=4.8962\left({\chi }_{\mathrm{Pauling}}+0.2963\right).$$

$$K=\frac{1}{2}{\left(\frac{{\partial }^2{V}_{\mathrm{self}}}{\partial q^2}\right)}_{q=0}=\frac{1}{2}{J}_{ii},$$

$$K_{\text{empricial}}=\frac{1}{4\pi {ϵ}_o}{\left(\frac{N}{\alpha }\right)}^{\frac{1}{3}},$$

$$N=\frac{m}{\rho |R_H|}\cdot \frac{1}{N_A\mathrm{e}},$$

$$N_{\mathrm{eff}}=0.4\frac{m}{\rho |R_H|}\cdot \frac{1}{N_A\mathrm{e}}.$$

$$f(r)=\frac{f_e\exp \left[-\beta \left(\frac{r}{r_e}-1\right)\right]}{1+{\left(\frac{r}{r_e}-\lambda \right)}^{20}},$$

$$K=\frac{1}{2}{J}_{ii}=\frac{1}{2}\int \int d^3\mathbf{r}{d}^3{\mathbf{r}}^{\prime }{n}_i(\mathbf{r})\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{n}_i({\mathbf{r}}^{\prime }),$$

$$n(\mathbf{r})=\frac{{\tau }_{\alpha }^3}{8\pi }{e}^{-{\tau }_{\alpha }|\mathbf{r}-\mathbf{R}|},$$

$$U_{\alpha }=\frac{1}{4\pi {ϵ}_0}\frac{5}{16}{\tau }_{\alpha },$$

$$K_{\mathrm{Slater}}=\frac{1}{4\pi {ϵ}_0}\frac{5}{16}\frac{\beta }{r_e}.$$

$$f_j(r,q_j)\to \left(1-\frac{q_j}{N_j}\right)f_j(r),$$

$${\varphi }_{ij}(r)=\frac{1}{2}\left[\frac{f_j(r)}{f_i(r)}{\varphi }_{ii}(r)+\frac{f_i(r)}{f_j(r)}{\varphi }_{jj}(r)\right].$$

$$g(\left\{q\right\})=\sum _i^N{q}_i-Q,$$

$$\mathcal{L}=\sum _{i=1}^N\left[\frac{1}{2}{m}_i{\dot{\mathbf{r}}}_i^2+\frac{1}{2}{M}_q{\dot{q}}_i^2\right]-V\left(\left\{\mathbf{r}\right\},\left\{q\right\}\right)-\lambda \left(\sum _{i=1}^N{q}_i-Q\right),$$

$$m_i{\ddot{\mathbf{r}}}_i={\text{F}}_i=\frac{\partial \mathcal{L}}{\partial {\text{r}}_i}=-\frac{\partial V}{\partial {\mathbf{r}}_i},$$

$$M_q{\ddot{q}}_i=f_i=\frac{\partial \mathcal{L}}{\partial q_i}=-\frac{\partial V}{\partial q_i}-\lambda ,$$

$$\lambda =-\frac{1}{N}\sum _{i=1}^N\frac{\partial V}{\partial q_i}.$$

$${\mathbf{r}}_1$$

$${\mathbf{r}}_2$$

$${\mathbf{r}}_1$$

$${\mathbf{r}}_2$$

$${\mathbf{r}}_3$$

$${\mathbf{r}}_4$$

$$V_{\mathrm{ext}}=-\sum _{i=1}^N\mathrm{e}{q}_i{\mathbf{r}}_i\cdot \mathbf{E}.$$

Full text

I Introduction

The metallic models normally used in molecular dynamics simulations of bulk metals (EAM,Daw and Baskes (1984); Foiles et al. (1986); Johnson (1989); Daw (1989); Plimpton and Hendrickson (1993); Voter (1995); Lu and Szpunar (1997); Alemany et al. (1998) Finnis-Sinclair,Finnis and Sinclair (1984); Sutton and Chen (1990) MEAM,Baskes and Johnson (1994); Lee and Baskes (2000); Thijsse (2002); Timonova and Thijsse (2011); Lee et al. (2001); van Beurden et al. (2002) and Quantum Sutton-ChenKimura et al. (1998); Qi et al. (1999)) have all been widely used by the materials simulation community for work on bulk and nanoparticle properties,Chui and Chan (2003); Wang et al. (2005); Medasani et al. (2007) melting,Belonoshko et al. (2000); Sankaranarayanan et al. (2006, 2005) fracture,Shastry and Farkas (1996, 1998) crack propagation,Becquart et al. (1993) and alloying dynamics.Shibata et al. (2002) One of the strengths common to all of the methods is the relatively large library of metals for which these potentials have been parameterized.Foiles et al. (1986); Johnson (1988); Rifkin et al. (1992); Mishin et al. (1999, 2001, 2002); Zope and Mishin (2003); Mishin et al. (2005); Zhou et al. (2001, 2004a) However, none of these models allow the metal atoms to polarize, so they neglect a vital interaction with ionic or polar groups.

Streitz and Mintmire developed an electrostatic+EAM (ES+) potential energy model for aluminum oxide surfaces,Streitz and Mintmire (1994) which combined a variable-charge electrostatic approach with the Finnis-Sinclair variant of EAM Finnis and Sinclair (1984). The ES+ potential was originally parameterized only for Aluminum Oxide phases, but has recently been adapted to other materials. Charge-equalization approaches based on the extended Lagrangian charge equalization method pioneered by Rick, Stuart, and Berne,Rick et al. (1994, 1995); Rick and Stuart (2002) have also been utilized to parametrize similar EAM+Charge Transfer models for Cu/\ceCuO interfaces,Devine et al. (2011) titanium/\ceTiO2 interfaces,Ogata et al. (1999) oxides of the ternary Al-Ni-Fe alloys,Jeon et al. (2011) and binary Al-Zr oxides.Zhou et al. (2004b) Another charge transfer ionic-embedded atom method potential was parametrized by Zhou et al. to study the growth of \ceO–\ceAl–\ceNi–\ceCo–\ceFe system.Zhou and Wadley (2005)

Our primary interest is the interactions between metals and nonreactive, but polar, species (e.g., water, carbon monoxide, ions, etc.). These interactions have a large role in the interfacial thermal conductance, the ordering of water on metal surfaces, surface restructuring under gas overpressure, surface friction, and slip/stick hydrodynamics. A charge or multipole that comes close to a conducting surface creates a disturbance in the valence electron density in the conductor. The charge density is altered only at the surface of the conductor, but the effective interactions are often treated using image charges or multipoles which are located inside the metal. To accurately capture this effect, EAM must be modified to handle perturbations in the valence electron density due to the presence of external species. The goal is an atomistic model of a metal which can accurately reproduce the image charge effects exhibited by real metal interfaces. Although the ES+ and related approaches that add fluctuating charges to EAM may exhibit some of these features, they have been tuned to bulk-like properties for fully or partially oxidized metals, and not for surface interactions and screening of metal interfaces. There are two significant assumptions made by the ES+ model:

The fluctuating charges interact primarily via Coulomb integrals between Slater-type orbitals centered on atomic sites. 2. 2.

The self-energy for modifying a charge on a site is essentially a parabolic function which depends on the ionization potential (IP) and electron affinity (EA) of the neutral species.

We present here a model in which the core/valence distinction of EAM-like models is left in place, and the fluctuating charges alter the valence densities, which are relatively diffuse in EAM-based models. We also parametrize the “self” potential as a sixth-order polynomial, using information from thermodynamically-derived models for charge transfer (i.e. Pauling electronegativities), experimental measures of charge mobility in metals (Hall coefficients), and experimental ionization data beyond the $+1$ and $-1$ oxidation states.

II Methodology

For a collection of atoms with instantaneous positions, $\{\mathbf{r}\}$, and partial charges, $\{q\}$, the configurational potential energy in DR-EAM is given by

[TABLE]

Here, $F_{i}\left[\bar{\rho}_{i}\right]$ is an energy functional for embedding atom $i$ in a total valence density, $\bar{\rho}_{i}$, located at $\mathbf{r}_{i}$, the position of the pseudo-atom $i$ (nucleus + core electrons). $\phi_{ij}$ is a pair potential that represents the (mostly) repulsive overlap of the two pseudo-atom cores at a distance $r_{ij}=\left|\mathbf{r}_{j}-\mathbf{r}_{i}\right|$. $J\left(r_{ij}\right)$ is the Coulomb integral that accounts for electrostatic contributions from the fluctuations in the valence charge density distributions, and $V_{\mathrm{self}}$ is the energy for under- or over-charging each atom.

The instantaneous electron density due to the valence electrons from all the other atomic sites is computed at the location of each atom. For atom $i$,

[TABLE]

where $f_{j}\left(r\right)$ is the radial dependence of the valence density of neutral atom $j$, and $q_{j}$ is a dynamic charge variable that governs the instantaneous fluctuations in the valence density. $N_{j}$ is a “valency count” for atom $j$ that is determined by the number of free charge carriers in the bulk metal. Changes in the partial charge value allow for conduction band electrons $(q<0)$ or holes $(q>0)$ to migrate into a spatially localized cloud surrounding each atom.

The pair potential, $\phi_{ij}$, in DR-EAM also depends on the instantaneous valence densities at sites $i$ and $j$.

[TABLE]

where $\phi_{ii}$ is the pair interaction of two $i$ atoms in the pure bulk metal.

A treatment of electrostatic interactions is required to account for local perturbations to the background electron density. Here, we have adopted the damped shifted force (DSF) kernelFennell and Gezelter (2006),

[TABLE]

which has energies and forces that go smoothly to zero approaching a cutoff value, $R_{\textrm{c}}$. The damping parameter, $\alpha$, describes the effective screening length of the charge, essentially treating density perturbations as Gaussians of width $\alpha^{-1}$. Other electrostatic models, including analytical integration of the Slater Coulomb integrals, have also been adopted by other fluctuating charge approaches to polarization.Zhou and Wadley (2005); Rick et al. (1994)

II.1 The self potential

The self potential in DR-EAM accounts for the energetic penalty for over-charging or under-charging a neutral atom, and is modeled with a polynomial,Iczkowski and Margrave (1961)

[TABLE]

The parameters have been tuned using a range of electron affinities and ionization potentials for commonly exhibited oxidation states in bulk materials,Hotop et al. (1973); Bilodeau et al. (1998); Andersen et al. (1999); Scheer et al. (1998a); Chen and Ning (2016a); Bilodeau et al. (1999); Scheer et al. (1998b); Chen and Ning (2016b); Chen et al. (2016); Feigerle et al. (1981); Lindahl et al. (2010); Bratsch and Lagowski (1986); Ilin et al. (1987); Fu et al. (2017) as well as the atomic polarizabilities of the neutral metals. See Figure. 1 for an example of the self potential in both the DR-EAM and harmonic models.

The standard electronegativity equalization models, e.g., charge equilibration (QEq)Rappé and Goddard III (1991) and fluctuating charge (fluc-q)Rick et al. (1994), treat inter-molecular fluctuating charge interactions with a Coulombic potential, but use Slater Coulomb integrals for intra-molecular interactions, and include harmonic self potentials. The two parameters that are usually used to describe the harmonic self energy for atom $i$ are the Mulliken electronegativity, $\chi_{i}$ and $J_{ii}$ (twice the chemical hardness),

[TABLE]

[TABLE]

where energies at fixed charge states, ($q=0,\pm 1$), i.e., the ionization potential (IP) and electron affinity (EA), are used to set both of these parameters.Rappé and Goddard III (1991) However, there is now significant evidence that the charge states in the condensed phase are better represented by partial charges that are significantly smaller. Ab initio calculations of ionic liquids show the charges of ions in ionic liquids is typically between 0.6 to 0.9 in the units of electron charge.Zhang and Maginn (2012); Bhargava and Balasubramanian (2007); Schmidt et al. (2010); Wendler et al. (2011); Dommert et al. (2012) Also, a molecular dynamics study of charge transfer in 1,3-dimethylimidazolium bis(trifluoromethylsulfonyl) imide shows the average charge on ions is $\pm$ 0.84e.Ishizuka and Matubayasi (2016) Pluhařová et al. carried out an ab initio molecular dynamics simulation of fluoride and lithium ion pairing in water, and found the fractional charges in fluoride (-0.79 to -0.82e) and lithium (0.89 to 0.90e).Pluhařová et al. (2013)

More relevant to EAM-like models, Seriani et al. Seriani et al. (2007) found that in various form of Platinum Oxides, the calculated Bader charges on Pt atoms (in units of electron charge) are: 1.53 for $\alpha$-\cePtO2, 1.62 for $\beta$-\cePtO2, 1.13 for \cePt3O4 and 0.86 for \cePtO yielding an effective mapping of 0.4e for each successive oxidation state. For oxygen atoms, the charges remain roughly constant: -0.76 for $\alpha$-\cePtO2, -0.81 for $\beta$-\cePtO2, -0.85 for \cePt3O4 and -0.86 for \cePtO (all corresponding to a -2 oxidation state for oxygen). Wang et al. found the Bader charges in metal atoms for oxides of \ceTiO2, \ceMnO2, \ceCoO2, \ceNiO2, and \ceZrO2 to be 2.01e, 1.58e, 1.32e, 1.29e, 2.11e, respectively,Wang et al. (2014) providing an average charge of 0.415e for each oxidation state. For this reason, DR-EAM adopts an effective scaling of 0.4e for the partial charge on each successive oxidation state.

One additional consideration is the reliance on IP and EA values when other measures of electronegativity may be more useful. There are a number of approximate linear transformations between the Mulliken electronegativity and the Pauling electronegativity scale,Huheey (1978) which is widely used to predict polarization and charge transfer in the chemical literature. Because we seek to model polarization and charge transfer between a range of elements, we begin with the Pauling-Allred scale,Allred (1961) which is derived from thermodynamic data, and fit a linear relationship to the Mulliken scale (in eV) for all metals that have well-characterized electron affinity values. This relationship,

[TABLE]

is then used to set the first derivative of the self potential [$a_{1}$ in Eq. (5)]. Because we have adopted fractional partial charges (0.4e) to represent each oxidation state, the Mulliken electronegativity is related to the first-order coefficient ($a_{1}$) in the self potential by a factor of 2.5,

[TABLE]

For some metals (e.g., Mg), there are only approximate electron affinities available, so in order to obtain the first order coefficients in the self potential, we use tabulated Pauling electronegativity values to set $a_{1}$ for all of the parameterized elements.

Traditionally, the chemical hardness has been estimated via Coulomb integrals for Slater orbitals centered on the atomic sites.Rick et al. (1994) Because the total charge potential in standard fluc-q models is harmonic, the fluctuating charges then have a single low-energy solution, yielding unique sets of charges that satisfy electronegativity equalization at each atomic configuration.

In addition to overestimating of the partial charges carried by polarized sites in condensed phase simulations, we find that standard harmonic models of the self potential also overestimate the bulk polarizability (particularly for the coinage metals). One of the parameters $(a_{2})$ in the self energy (Eq. (5)) can be better estimated using either empirically derived atomic polarizability data, or via a Coulomb integral of only the valence charge density. The chemical hardness,

[TABLE]

has been correlated with the atomic polarizablity $\alpha$ via an empirical relationship,Nagle (1990)

[TABLE]

where $N$ is the number of charge carriers per atom. The number of charge carriers per atom can be experimentally determined using the Hall effect. The Hall effect is the production of a potential difference on a current-carrying conductor when a magnetic field is applied perpendicular to the current. The Hall coefficient, $R_{H}$, a constant relating the magnetic field, current, and the potential produced, can be used to find the number of charge carriers. In DR-EAM, charge carrier densities are determined using Hall coefficients ($R_{H}$),

[TABLE]

where $\rho$, and $m$ are density of the metal, and the atomic mass, respectively, and $N_{A}$ is Avogadro’s number. Each of the charge states has been scaled by 0.4, so, the number of effective charge carriers for atom $i$ is

[TABLE]

In DR-EAM, $N_{\mathrm{eff}}$, is also taken to be the effective valence count in Eq.(2).

Fluctuating charge self potentials based on IP and EA can also overestimate the polarizablity of metals, sometimes badly enough to trigger a polarization catastrophe. This happens because the Coulombic interaction is bilinear in the charge degrees of freedom $(q_{i}\times q_{j})$ so charges with opposing signs will tend to amplify their differences without a fast-rising self-potential to bound the charges. If the curvature of a harmonic self-potential is too small, the Coulombic interaction will dominate, leading to a polarization catastrophe.

Atomic polarizablity depends on the chemical hardness of the atom, which is related to the coefficient of second order in the self potential. In DR-EAM, the self potential is determined by imposing the constraint that the coefficient of second order is equal to the hardness of atom. The scaling of charge states and the use of a higher order polynomial helps solve the overpolarization issues. A fourth or sixth-order polynomial with a positive coefficient in the highest order term will eventually rise faster than the Coulombic energy falls. The polarization catastrophe was solved by Zhou et al. in their modified charge transfer-EAM method using upper and lower bounds on the charge values that oxygen and metals can take.Zhou et al. (2004b); Zhou and Wadley (2005) We find that scaled charges and a higher order polynomial are sufficient to prevent the polarization catastrophe in DR-EAM.

Zhou et al.Zhou et al. (2004a) proposed EAM parameters and functionals for 16 metals. In this paper, we are implementing DR-EAM based on these parameters. The valence electron density in the Zhou et al. parametrization takes the form

[TABLE]

where the denominator acts as a cutoff function that takes $f(r)$ smoothly to zero at $r/r_{e}$ values ranging from 2.5 to 3.5. The valence electron density without the denominator has form of a Slater function.

In the local density approximation (LDA), there is another way of estimating the chemical hardness from the electronic density of isolated atoms.Elstner et al. (1998)

[TABLE]

where $n_{i}(\textbf{r})$ is the electronic density of the isolated atom. Elstner et al.Elstner et al. (1998) used normalized spherical charge densities,

[TABLE]

where $\mathbf{R}$ is the location of the atom, to derive the chemical hardness of a spin-unpolarized atom or Hubbard parameters $U_{\alpha}$ in terms of $\tau_{\alpha}$.

[TABLE]

We have adopted the Elstner et al. result along with the Slater-form of the EAM densities in Eq. (14). We find that by setting $\tau_{\alpha}=\beta/r_{e}$, a chemical hardness appropriate for use the self potential can come directly from the valence density functions,

[TABLE]

The valency count, $N_{\mathrm{eff}}$, and the hardness, $K$, values that were used to parametrize DR-EAM are provided in Table 1.

The second-order coefficient ($a_{2}$) is fixed to value of chemical hardness ($K$) obtained from Eq. (18). All remaining coefficients are determined by fitting gas phase electron affinities and ionization potentials where the charges have been scaled to match the 0.4e per oxidation state in the condensed phase. While performing the fits, the coefficients of the highest even order in the polynomial (typically $q^{6}$) were also constrained to be positive. Self-potential parameters are given in Table 2, and plots of the DR-EAM self-potentials are given in Figure. 2.

II.2 Differences from standard EAM

The embedding functional and pair potential are not unique in EAM, and it is possible to parametrize reasonable metal potentials by transferring some weight between the embedding functional and the effective pair potential. In DR-EAM, we have adopted the same form and parametrization for the embedding functional in DR-EAM adopted by Zhou et al.Zhou et al. (2001); Wadley et al. (2001); Zhou et al. (2004a). Although the original EAM potentials have only positional degrees of freedom that alter the local valence density, DR-EAM adds an additional charge degree of freedom to each atomic site. The system evolves with the time-dependent positions and partial charges, readjusting the valence density contributed by each atom. The adjusted valence density in DR-EAM,

[TABLE]

where $q_{j}$ and $N_{j}$ are the partial charge and valency count for atom $j$.

The form of the pair potential is tied to the choice of embedding functional. It is possible to transform $F$ and $\phi$ so that the slope of the embedding function is zero at the equilibrium electron density. Mixed-element pair potentials $(\phi_{ij})$ can be defined many ways in terms of pair potentials for the two elements separately ($\phi_{ii}$ and $\phi_{jj}$). Arithmetic means and geometric means have both been used previously,Daw (1989) but when the slope of the embedding function is set to zero at the equilibrium electron density, the effective two body potentials $(\phi)$ are negative at some distances, and the geometric mean cannot be used. Johnson derived a mixed-element pair potential for alloys in terms of the valence densities,Johnson (1989)

[TABLE]

This form of the pair potential is independent of the arbitrary electron-density transformation. Since DR-EAM uses the readjusted valence density, the pair potential in DR-EAM is given by Eq. (3). The self energy, Vself, and the Coulombic interactions could also be included in a modified functional, but for simplicity, DR-EAM treats the energies due to alteration of net atomic charges as separate terms in the potential energy.

II.3 Charge conservation

Because the fluctuating charge variables represent physical charge densities, they are not independent variables, and a charge conservation constraint is required when propagating these degrees of freedom along with the nuclear coordinates. Charge conservation is implemented using the method of undetermined multipliers to enforce the constraints. The net charge constraint is written as

[TABLE]

where $Q$ is the (fixed) charge on the system in units of electron charge and $q_{i}$ is the fluctuating charge of the $\mathrm{i^{th}}$ atom. Each contiguous metallic region can be thought of as a single molecule, so it is also possible to constrain metallic regions separately using several undetermined multipliers. However, using multiple constraints prevents charge transfer between the separate blocks.

II.4 Dynamics

The minimum energy of the fluctuating charge system is first determined using steepest-descent minimization on the configurational energy of the system with frozen nuclear coordinates. Once an optimal set of charges is known, perturbations to these charges depend only on nuclear coordinates (which are relatively slowly moving degrees of freedom). For this reason, propagation of the fluctuating charge variables using extended-Lagrangian simulations has become widespread.

The extended Lagrangian in DR-EAM (with the charge constraint) is given by

[TABLE]

where $m_{i}$ is the mass of the $i^{\mathrm{th}}$ atom, $M_{q}$ is the fictitious mass assigned to the fluctuating charges and $\lambda$ is Lagrange multiplier for the charge constraint. The first term in the kinetic energy is the contribution from the motion of atoms themselves while the second term is a fictitious kinetic energy, the contribution from the charge velocities $\left\{\dot{q}\right\}$. $V\left(\left\{\mathbf{r}\right\},\left\{q\right\}\right)$ is the configurational energy for the extended system defined in Eq. (1).

Equations of motion for the fictitious charge $q$ and nuclear coordinates $\mathbf{r}$ can be obtained by calculating the forces on the dynamical charge variables and on the atoms themselves,

[TABLE]

[TABLE]

where $\lambda$ is the charge conservation constraint,

[TABLE]

Two temperatures, nuclear and electronic, are defined for the system. The nuclear temperature arises from the kinetic energy of atoms whereas the electronic temperature comes from the fictitious kinetic energy of the fluctuating charges. Perturbations in the fluctuating charge forces will be due primarily to the motion of surface-adsorbed molecules (e.g., water), so the timescale for charge fluctuations should be approximately the same as for nuclear motion. For this reason, $M_{q}$ is chosen to be large enough ($M_{q}=600\mathrm{~{}kcal~{}mol^{-1}~{}fs^{2}~{}e^{-2}}$) so that changes in electronic degrees of freedom can be integrated along with nuclear coordinates. In theory, if $M_{q}$ is chosen to be very small, it would be possible to simulate the collective electronic motion (e.g., plasmons), but the time steps for those simulations would by necessity be extremely small relative to molecular motion.

The charge and nuclear degrees of freedom are coupled by the DR-EAM potential energy in Eq. (1). Although it is possible to propagate the entire system in the microcanonical (NVE) ensemble, equipartition eventually brings the electronic temperature to the same value as the nuclear temperature. Electronegativity equalization methods assume a single low-energy solution to the charge degrees of freedom, i.e., an electronic temperature of zero, and there are many ways to achieve this. One could minimize the energy with respect to the charge degrees of freedom at every time step, but this would be prohibitively expensive. It is also possible to propagate the nuclear coordinates in the microcanonical (NVE) ensemble while simultaneously keeping the electronic degrees of freedom at a much lower temperature, $T_{e}<1$ K, with a Nosé-Hoover thermostat. In many of the tests and simulations described below, the electronic degrees of freedom were kept at a low temperature using Langevin dynamics, with a drag coefficient of $0.1\mathrm{~{}kcal~{}mol^{-1}~{}fs~{}e^{-2}}$. As in any fluctuating charge method, the coupling between the electronic and nuclear degrees of freedom will eventually transfer energy between the subsystems, so we expect that separately thermostatting the electronic coordinates will be required.

We note that the fluctuating charge forces modified by the charge conservation constraint are used in all steps of this method, including the initial minimization to find the optimal set of charges.

III Tests and Applications

To test the new method, we carried out simulations using DR-EAM where the choice of the underlying EAM functions and parameters were taken from Zhou et al.Zhou et al. (2004a) as a basis for further refinement. We have tested the new method on: (1) pure bulk metals, (2) common surfaces of the bulk metals, (3) ordered structures (L$1_{0}$ and L$1_{2}$), (4) a metal slab in a uniform electric field, and (4) metal slabs with external fixed charges approaching the surface. In all tests, electrostatic interactions were calculated with the damped shifted force (DSF) method with a damping parameter ($\alpha$) of 0.14 $\mathrm{\SIUnitSymbolAngstrom}$-1.

III.1 Bulk metals

The cohesive energy per atom, vacancy formation energy, bulk modulus, shear modulus, and Poisson’s ratio were computed for both DR-EAM and the unmodified EAM energy function, and these are provided in Table 3. Elastic stiffness tensors were calculated using the Energy versus Strain method of Yu et al.Yu et al. (2010); Golesorkhtabar et al. (2013) and elastic constants were computed using the same definitions that are utilized by the Materials Project.de Jong et al. (2015); Gaillac et al. (2016)

In all of the parameterized elements, the bulk metals do not polarize to any significant degree. There is negligible (nearly zero) charge transfer between the atoms in pure metals, so DR-EAM and EAM produce identical bulk properties within the limit of precision used in calculation.

III.2 Surface energies of bulk metals

In DR-EAM, the anisotropic environment around atoms at the surface of a metal encourages charge transfer into the near-surface atoms to allow the surface atoms to approach the equilibrium valence density ($\rho_{e}$). In traditional EAM, the only way for atoms to approach equilibrium densities is by contracting the surface layers closer to the underlying bulk.

To study this difference, the surface formation energies for (111), (110), and (100) surfaces are computed for FCC and BCC metals using both DR-EAM and EAM. For HCP metals, the basal plane (0001) was exposed to calculate the surface energy. This data is provided in Table 4. For most of the surfaces, DR-EAM and EAM again produce identical surface energies within the limit of precision used in the calculations. Notable exceptions are the Ag(110), Al(100), and Al(110) surfaces. Even for these surfaces, however, the surface polarization yields only a small charge, $\langle q\rangle<0.0025$, on the atoms in the subsurface layers. Surface atoms exhibit even smaller charges.

III.3 Ordered structures

Ordered alloys are composed of at least two elements which alternate in position in a regular pattern. Two ordered structures (L$1_{2}$ and L$1_{0}$) were used to study how DR-EAM predicts charge transfer and structural changes in alloys (see Figure. 3). L$1_{0}$ alloys have \ceAB composition and space group P4/mmm, which exhibits a layered structure. The basis vectors for L$1_{0}$ with two components, A and B, are

[TABLE]

where a and c are the lattice parameters. The alternation of A and B layers lies along the $c$ direction, while all atoms along the two directions with the $a$ lattice parameters are identical. This means that charge transfer between the two elements will result in modification of the $c/a$ ratio exhibited by the crystal.

L$1_{2}$ alloys have a \ceAB3 composition and space group Pm$\bar{\mathrm{3}}$m. The basis vectors for L$1_{2}$ with two components, A and B, are

[TABLE]

where a is the lattice parameter. Each of the A atoms in an L$1_{2}$ structure is surrounded by a symmetric set of B atoms, so charge transfer between A and B will result in isotropic volume contraction due to electrostatic interactions.

Ordered structures were formed by replacing appropriate atoms in a bulk FCC lattice of A atoms with B atoms. Following this replacement, the geometry of the periodic box was optimized, while allowing the fluctuating densities in DR-EAM to re-optimize simultaneously with each change in the geometry.

Bulk properties of the L$1_{0}$ and L$1_{2}$ ordered alloys (lattice constants, bulk modulus, shear modulus, Poisson’s ratio and energy per atom) are compared with DFT calculations and experimental values (when known) in Tables 5 and 6. Significant charge transfer was observed between the two components of the ordered structures. The direction and magnitude of the transfer is governed by the electronegativities ($a_{1}$) of the two elements. Even with charge transfer, for L12 structures, most of the lattice constants are nearly equal to the structures optimized under nonpolarizable EAM. This is likely due to cancellation of added A-B attraction with the simultaneous B-B repulsion in L$1_{2}$ structure. Even though the DR-EAM and EAM structures are similar, the bulk modulus, shear modulus, Possion’s ratio and energy per atom do exhibit small changes due to charge transfer. In the L$1_{0}$ structures (see Table 6), the DR-EAM optimized geometries are different from structures optimized without charge transfer. However, when compared with DFT structures or experimental lattice parameters, the two EAM methods are remarkably similar. Some notable exceptions (e.g., NiPt and AuCu) occur when charge transfer between the two metallic components is relatively large. Inter-metallic charge transfer in L$1_{0}$ alloys will usually increase in-layer Coulombic repulsion, while simultaneously increasing inter-layer attraction. This results in changes to the $c/a$ ratio that appear to depend largely on the extent of the charge transfer. However, very large charge transfers can also alter intermetallic pair potentials by changing the density ratios in Eq. 3, modifying the equilibrium geometry of the alloy.

III.4 Metal slabs in a uniform field

An atomic partial charge $q_{i}$ in an electric field $\mathbf{E}$ feels a physical force $\left(\mathbf{F}_{i}=eq_{i}\mathbf{E}\right)$ due to that field. To derive forces on the fluctuating charge variables due to the presence of an external field, we define a potential that depends on the location of each atom,

[TABLE]

This potential provides a fluctuating charge force due to the field,

[TABLE]

which depends on the position of the atom relative to the origin of the coordinate system. (A spatial derivative of $V_{\mathrm{ext}}$ also yields the correct physical force on the atom’s coordinates.)

We note that this uniform field is not realizable in any experiment, particularly when discussing a conducting slab. Uniform fields may be approximated in a region between charged plates, but interior to a conductor, the field will be shielded by the skin, and an external field should not interact with the atoms on the interior. However, it is useful to test the behavior of the new methodology on exposure to the field. If DR-EAM can approximate the correct behavior of a conductor, the response should be an effective cancellation of the external field inside a metal slab.

A local electric dipole density can be defined between consecutive atomic layers. The origin of a local coordinate system ($O_{n}$) is set between layers $n$ and $n+1$, and the local dipole density between these layers is measured with respect to this origin,

[TABLE]

where $\langle q_{n}\rangle$ is the average charge in the nth layer, $\mathbf{R}_{n}$ and $\mathbf{R}_{n+1}$ are the average positions of the two layers with respect to $O_{n}$, and $v$ is the volume between the two layers $(v=L_{x}L_{y}\Delta z=A\Delta z)$.

To test the shielding properties of DR-EAM, an electric field was applied along the $z$-axis of a slab comprising 18 atomic layers of each metal. Because the local origin is located half-way between each layer, the polarization density,

[TABLE]

depends only on the average charges on the two layers. The bound charge density inside the slab, $\rho_{b}=-\nabla\cdot\mathbf{P}$. Since our external field (and polarization density) both point along the $z$-axis, we can simply write the bound charge density, $\rho_{b}=-\partial P/\partial z$.

When an electric field is applied to a Pt slab with 18 layers, the resulting charge distribution in DR-EAM screens the electric field in the interior of the slab. The average charges in each layer are approximately linear in $z$ in the interior of the slab, which yields a nearly constant dipole density in the interior, and nearly zero bound charge density inside the slab. Figure. 4 shows the average charge in each layer as well as the local dipole density and screened electric field values at the surface and interior of the slab. These are the expected responses of a classical conductor in an external field. The electric field penetrates up to three layers deep from the atomic skin, but is fully screened beyond that.

III.5 Image charge effects

When a ion is placed directly above a conducting metal slab, the ion polarizes the slab and an oppositely charged image is induced on the surface of the slab directly beneath the ion. Classical treatments of charges at the surfaces of infinite planar conductors predict an effective interaction of the ion with its own image, $U_{\mathrm{eff}}(d)=-q^{2}\mathrm{e}^{2}/(4\pi\epsilon_{0}\,2d)$, where $d$ is the separation of the ion from the surface of the conductor, and $q$ is the charge on the ion (in units of electrons, $e$).

When computing electrostatic interactions between atoms in the DR-EAM model and external fixed charges, we have used the Damped Shifted Force (DSF) potential. The electrostatic kernel in DSF is given by Eq. (4), so it is reasonable to fit the effective potential between the unit point charge ($q=\pm 1$) and the polarizable DR-EAM surface with a similarly damped interaction model,

[TABLE]

where $\alpha$ is the value of the damping coefficient (in $\mathrm{\SIUnitSymbolAngstrom}$-1) used in the calculation. If the DR-EAM model captures the physics of classical image charges, the fitting parameter $s$ would be expected to be exactly two for infinite planar conducting surfaces, $\delta$ measures the offset of the surface of zero potential from the atomic layer at the top of the slab (see Figure. 5). For a perfectly flat classical conductor, we would expect $\delta=0$.

DR-EAM does exhibit image charge effects (shown graphically in Figure. 6). We placed a point charge models of chloride $(q=-1)$ ions above a range of metal surfaces, including Copper, Gold, and Platinum (111) and (100) surfaces. The ions are brought close to the surface above common sites (atop, hollow, and bridge) on these surfaces and the partial charges are allowed to find their lowest energy configurations. The surfaces polarize directly beneath the probe charge, and the distance dependence of the total interaction is used to fit $d$, $s$, and $\delta$ values in Eq. (30). For all of the surfaces presented in Table 7, we note that the scaling approaches the classical conductor value $s\rightarrow 2$, and $\delta$ sits approximately 0.5 layers beneath the atomic coordinates of the surface atoms.

IV CONCLUSIONS

DR-EAM is a polarizable force field for metals where each metal atom has an additional variable, or partial charge, that represents a fluctuation in the local valence density on that atom. The dynamics of the partial charges are calculated using an extended Lagrangian, and can aid in simulating polarization and charge transfer effects in systems that contain these metals. The polarization catastrophe in traditional electronegativity equalization models is solved using a sixth order polynomial for the self potential. The coefficients of this polynomial are tied to thermodynamically- and experimentally-derived data, notably the Pauling-Allred electronegativities, Hall coefficients, higher ionization energies, and bulk polarizabilities of these metals.

A number of important physical properties of metals are captured by DR-EAM. The partial charges distribute on the surface in response to external fields and produce screening of the electric field inside a metal slab. In addition, the surface charge response reproduces the classical image charge effect for point charges. These effects cannot be modeled using nonpolarizable versions of the embedded atom method. Additionally, all of the strengths of EAM (i.e. reasonably quantitative results for bulk elastic constants and surface energies) have been retained in DR-EAM. Relative to an unmodified EAM simulation, the extended Lagrangian in DR-EAM adds an additional 12.5% to the computational cost associated with a 2000 atom, 1 ns simulation. This is a relatively small additional cost given the new capability to simulate interactions with charged and polar molecules adsorbed on metal surfaces.

Acknowledgements.

Support for this project was provided by the National Science Foundation under Grant No. CHE-1663773. Computational time was provided by the Center for Research Computing (CRC) at the University of Notre Dame.