1 // EMTDefaultParameterProvider.cpp -- provides the default EMT parameters.
2
3 // Copyright (C) 2008-2011 Jakob Schiotz and Center for Individual
4 // Nanoparticle Functionality, Department of Physics, Technical
5 // University of Denmark. Email: schiotz@fysik.dtu.dk
6 //
7 // This file is part of Asap version 3.
8 // Asap is released under the GNU Lesser Public License (LGPL) version 3.
9 // However, the parts of Asap distributed within the OpenKIM project
10 // (including this file) are also released under the Common Development
11 // and Distribution License (CDDL) version 1.0.
12 //
13 // This program is free software: you can redistribute it and/or
14 // modify it under the terms of the GNU Lesser General Public License
15 // version 3 as published by the Free Software Foundation. Permission
16 // to use other versions of the GNU Lesser General Public License may
17 // granted by Jakob Schiotz or the head of department of the
18 // Department of Physics, Technical University of Denmark, as
19 // described in section 14 of the GNU General Public License.
20 //
21 // This program is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // and the GNU Lesser Public License along with this program. If not,
28 // see <http://www.gnu.org/licenses/>.
29
30
31 #include "EMTDefaultParameterProvider.h"
32 #include "Exception.h"
33 #include <string.h>
34 #include <cmath>
35 #include <iostream>
36
37 using namespace std;
38
39 const int EMTDefaultParameterProvider::shell0 = 3;
40 const int EMTDefaultParameterProvider::shell1 = 4;
41 const double EMTDefaultParameterProvider::bohr = 0.5291772; // Angstrom
42
43 typedef vector<emt_parameters *>::iterator param_p_iter;
44
EMTDefaultParameterProvider()45 EMTDefaultParameterProvider::EMTDefaultParameterProvider()
46 {
47 chi = NULL;
48 listcutofffactor = 1.04500185048;
49 }
50
~EMTDefaultParameterProvider()51 EMTDefaultParameterProvider::~EMTDefaultParameterProvider()
52 {
53 if (chi != NULL) delete chi;
54 for (param_p_iter ep = params.begin(); ep != params.end(); ++ep)
55 {
56 emt_parameters *e = *ep;
57 delete e;
58 }
59 }
60
Debug()61 void EMTDefaultParameterProvider::Debug()
62 {
63 cerr << "EMTDefaultParameterProvider debug information:" << endl;
64 cerr << "Length of params vector: " << params.size() << endl;
65 for (param_p_iter e = params.begin(); e != params.end(); ++e)
66 cerr << " " << (*e)->name << " Z=" << (*e)->Z << endl;
67 if (chi == 0) {
68 cerr << "Chi matrix: NOT ALLOCATED." << endl;
69 } else {
70 int n = params.size();
71 cerr << "Chi matrix: " << n << " x " << n << endl;
72 for (int i = 0; i < n; ++i)
73 for (int j = 0; j < n; ++j)
74 cerr << " chi[" << i << "][" << j << "] = " << (*chi)[i][j] << endl;
75 }
76 }
77
GetParameters(int element)78 const emt_parameters *EMTDefaultParameterProvider::GetParameters(int element)
79 {
80 emt_parameters *newpar;
81 // Check if we already have parameters for this element. If so:
82 // return them.
83 for (param_p_iter i=params.begin(); i != params.end(); ++i) {
84 if ((*i)->Z == element)
85 return *i;
86 }
87
88 // Get the parameters and add them to the vector.
89 newpar = GetNewParameters(element);
90 newpar->index = params.size(); // Count the elements, starting at zero.
91 params.push_back(newpar);
92 return newpar;
93 }
94
GetNewParameters(int element)95 emt_parameters *EMTDefaultParameterProvider::GetNewParameters(int element)
96 {
97 // This version returns the parameters using the units Angstrom
98 // and eV. To prevent typos the parameters are not converted in
99 // the if statement below, instead the conversion from bohr to
100 // Angstrom is done when the values are copied to the
101 // emt_parameters struct.
102
103 double E0, S0, n0, V0, eta2, kappa, lambda, mass;
104 double ls;
105 int Z;
106 string name;
107 emt_parameters *p;
108
109 if(element == 13){
110 name = "Al";
111 E0=-3.280; S0=3.000; V0=1.493; Z=13; eta2=1.240; kappa=2.000;
112 lambda=1.169; mass=26.98 ; n0=0.007 ; ls = 7.54871784;
113 }
114 else if(element == 29){
115 name = "Cu";
116 E0=-3.510; S0=2.67 ; V0=2.476; Z=29; eta2=1.652; kappa=2.74;
117 lambda=1.906; mass=63.54 ; n0=0.0091 ;
118 ls = 6.789382809;
119 }
120 else if(element == 47){
121 name = "Ag";
122 E0=-2.96 ; S0=3.01 ; V0=2.132; Z=47; eta2=1.652; kappa=2.790;
123 lambda=1.892; mass=107.87; n0=0.00547;
124 ls = 7.6790043;
125 }
126 else if(element == 79){
127 name = "Au";
128 E0=-3.80 ; S0=3.00 ; V0=2.321; Z=79; eta2=1.674; kappa=2.873;
129 lambda=2.182; mass=196.97; n0=0.00703;
130 ls = 7.66504117182;
131 }
132 else if(element == 28){
133 name = "Ni";
134 E0=-4.44 ; S0=2.60 ; V0=3.673; Z=28; eta2=1.669; kappa=2.757;
135 lambda=1.948; mass=58.71 ; n0=0.0103 ;
136 ls = 6.598896;
137 }
138 else if(element == 46){
139 name = "Pd";
140 E0=-3.90 ; S0=2.87 ; V0=2.773; Z=46; eta2=1.818; kappa=3.107;
141 lambda=2.155; mass=106.4 ; n0=0.00688;
142 ls = 7.330378;
143 }
144 else if(element == 78){
145 name = "Pt";
146 E0=-5.85 ; S0=2.90 ; V0=4.067; Z=78; eta2=1.812; kappa=3.145;
147 lambda=2.192; mass=195.09; n0=0.00802;
148 ls = 7.41119853;
149 }
150 else if(element == 12){
151 name = "Mg";
152 // these Mg parameters have been generated by N. Bailey using properties
153 // of Mg derived from DFT, as well as properties of Mg-Cu alloys (Cu
154 // parameters were also refit, but the defaults will be left as is)
155 // The fitting process used eV and Angstrom, so the reverse of the
156 // scaling below needs to be done here for parameters involving length
157 // units. The fitting was done in December 2002
158 E0=-1.487 ; S0=1.7664 / bohr ; V0=2.2298; Z=12; eta2=2.5411*bohr; kappa=4.435 * bohr;
159 lambda=3.2927 * bohr; mass=24.3050; n0=0.03554 * (bohr*bohr*bohr);
160 ls = 4.52004 / bohr;
161 }
162 else
163 {
164 throw AsapError("This element isn't defined in EMT.");
165 }
166
167 p = new emt_parameters;
168 p->e0 = E0;
169 p->seq = S0 * bohr;
170 p->neq = n0 / (bohr*bohr*bohr);
171 p->V0 = V0;
172 p->eta2 = eta2 / bohr;
173 p->kappa = kappa / bohr;
174 p->lambda = lambda / bohr;
175 p->mass = mass;
176 p->invmass = 1.0/mass;
177 p->gamma1 = 0.0; // These are calculated later!
178 p->gamma2 = 0.0;
179 p->Z = Z;
180 ASSERT(element == Z);
181 p->name = name;
182 p->lengthscale = ls / sqrt(2.0) * bohr;
183
184 return p;
185 }
186
CalcGammaEtc()187 void EMTDefaultParameterProvider::CalcGammaEtc()
188 {
189 // These can be overridden on an individual basis
190 calc_cutoff();
191 calc_gammas();
192 calc_chi();
193 }
194
GetMaxListCutoffDistance()195 double EMTDefaultParameterProvider::GetMaxListCutoffDistance() // Max value, useful before initialization.
196 {
197 maxseq = 3.01 * bohr; // Value for Ag
198 double maxcutoff = 0.5 * maxseq * Beta * (sqrt((double) shell0) +
199 sqrt((double) shell0 + 1));
200 return maxcutoff * listcutofffactor;
201 }
202
calc_cutoff()203 void EMTDefaultParameterProvider::calc_cutoff()
204 {
205 double r;
206
207 maxseq = 0.0;
208 for (param_p_iter i = params.begin(); i != params.end(); ++i)
209 if ((*i)->seq > maxseq)
210 maxseq = (*i)->seq;
211
212 cutoff = 0.5 * maxseq * Beta * (sqrt((double) shell0) +
213 sqrt((double) shell0 + 1));
214 // The cutoff slope is set (rather arbitrarily) to the same value
215 // used in the original ARTwork code. The idea of this
216 // expression is that the slope is such that the Fermi function
217 // has the value 1e-5 where the neighbor list is cut off.
218 r = cutoff * 2.0 * sqrt((double) shell0 + 1) / (sqrt((double) shell0) +
219 sqrt((double) shell0 + 1));
220 cutslope = log(9999.0) / (r - cutoff);
221 #if 0
222 cerr << "cutoff " << cutoff << endl;
223 cerr << "maxseq " << maxseq << endl;
224 cerr << "Beta " << Beta << endl;
225 cerr << "cutslope " << cutslope << endl;
226 cerr << "shell0 " << shell0 << endl;
227 cerr << "shell1 " << shell1 << endl;
228 #endif
229 }
230
calc_gammas()231 void EMTDefaultParameterProvider::calc_gammas()
232 {
233 static int shellpop[5] = {12, 6, 24, 12, 24}; // Population of the first 5 shells (fcc).
234 double w, d;
235
236 for (param_p_iter ep = params.begin(); ep != params.end(); ++ep)
237 {
238 emt_parameters *e = *ep;
239
240 e->gamma1 = e->gamma2 = 0.0;
241 // The neighborlist object does not return neighbors in th 4. shell!
242 // for (int is = 0; is < shell1 && is < 5; is++)
243 for (int is = 0; is < shell1 - 1 && is < 5; is++) // !!!!! xxxxx
244 {
245 d = sqrt((double) (is + 1)) * Beta * e->seq;
246 w = 1. / (1. + exp(cutslope * (d - cutoff)));
247 e->gamma1 += w * shellpop[is] * exp(-d * e->eta2);
248 e->gamma2 += w * shellpop[is] * exp(-d * e->kappa / Beta);
249 }
250 e->gamma1 /= shellpop[0] * exp(-Beta * e->seq * e->eta2);
251 e->gamma2 /= shellpop[0] * exp(-e->seq * e->kappa);
252 }
253 }
254
calc_chi()255 void EMTDefaultParameterProvider::calc_chi()
256 {
257 int n = params.size();
258
259 if (chi != NULL)
260 delete chi;
261 chi = new TinyDoubleMatrix(n,n);
262
263 for (int i = 0; i < n; ++i)
264 for (int j = 0; j < n; ++j) {
265 (*chi)[i][j] = params[j]->neq / params[i]->neq;
266 }
267 }
268
269