Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88483993
meam_setup_done.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Oct 19, 01:33
Size
26 KB
Mime Type
text/x-c
Expires
Mon, Oct 21, 01:33 (2 d)
Engine
blob
Format
Raw Data
Handle
21779410
Attached To
rLAMMPS lammps
meam_setup_done.cpp
View Options
#include "meam.h"
#include "math_special.h"
#include <algorithm>
using
namespace
LAMMPS_NS
;
void
MEAM
::
meam_setup_done
(
double
*
cutmax
)
{
int
nv2
,
nv3
,
m
,
n
,
p
;
// Force cutoff
this
->
cutforce
=
this
->
rc_meam
;
this
->
cutforcesq
=
this
->
cutforce
*
this
->
cutforce
;
// Pass cutoff back to calling program
*
cutmax
=
this
->
cutforce
;
// Augment t1 term
for
(
int
i
=
0
;
i
<
maxelt
;
i
++
)
this
->
t1_meam
[
i
]
=
this
->
t1_meam
[
i
]
+
this
->
augt1
*
3.0
/
5.0
*
this
->
t3_meam
[
i
];
// Compute off-diagonal alloy parameters
alloyparams
();
// indices and factors for Voight notation
nv2
=
0
;
nv3
=
0
;
for
(
m
=
0
;
m
<
3
;
m
++
)
{
for
(
n
=
m
;
n
<
3
;
n
++
)
{
this
->
vind2D
[
m
][
n
]
=
nv2
;
this
->
vind2D
[
n
][
m
]
=
nv2
;
nv2
=
nv2
+
1
;
for
(
p
=
n
;
p
<
3
;
p
++
)
{
this
->
vind3D
[
m
][
n
][
p
]
=
nv3
;
this
->
vind3D
[
m
][
p
][
n
]
=
nv3
;
this
->
vind3D
[
n
][
m
][
p
]
=
nv3
;
this
->
vind3D
[
n
][
p
][
m
]
=
nv3
;
this
->
vind3D
[
p
][
m
][
n
]
=
nv3
;
this
->
vind3D
[
p
][
n
][
m
]
=
nv3
;
nv3
=
nv3
+
1
;
}
}
}
this
->
v2D
[
0
]
=
1
;
this
->
v2D
[
1
]
=
2
;
this
->
v2D
[
2
]
=
2
;
this
->
v2D
[
3
]
=
1
;
this
->
v2D
[
4
]
=
2
;
this
->
v2D
[
5
]
=
1
;
this
->
v3D
[
0
]
=
1
;
this
->
v3D
[
1
]
=
3
;
this
->
v3D
[
2
]
=
3
;
this
->
v3D
[
3
]
=
3
;
this
->
v3D
[
4
]
=
6
;
this
->
v3D
[
5
]
=
3
;
this
->
v3D
[
6
]
=
1
;
this
->
v3D
[
7
]
=
3
;
this
->
v3D
[
8
]
=
3
;
this
->
v3D
[
9
]
=
1
;
nv2
=
0
;
for
(
m
=
0
;
m
<
this
->
neltypes
;
m
++
)
{
for
(
n
=
m
;
n
<
this
->
neltypes
;
n
++
)
{
this
->
eltind
[
m
][
n
]
=
nv2
;
this
->
eltind
[
n
][
m
]
=
nv2
;
nv2
=
nv2
+
1
;
}
}
// Compute background densities for reference structure
compute_reference_density
();
// Compute pair potentials and setup arrays for interpolation
this
->
nr
=
1000
;
this
->
dr
=
1.1
*
this
->
rc_meam
/
this
->
nr
;
compute_pair_meam
();
}
// ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
// Fill off-diagonal alloy parameters
void
MEAM
::
alloyparams
(
void
)
{
int
i
,
j
,
k
;
double
eb
;
// Loop over pairs
for
(
i
=
0
;
i
<
this
->
neltypes
;
i
++
)
{
for
(
j
=
0
;
j
<
this
->
neltypes
;
j
++
)
{
// Treat off-diagonal pairs
// If i>j, set all equal to i<j case (which has aready been set,
// here or in the input file)
if
(
i
>
j
)
{
this
->
re_meam
[
i
][
j
]
=
this
->
re_meam
[
j
][
i
];
this
->
Ec_meam
[
i
][
j
]
=
this
->
Ec_meam
[
j
][
i
];
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
];
// If i<j and term is unset, use default values (e.g. mean of i-i and
// j-j)
}
else
if
(
j
>
i
)
{
if
(
iszero
(
this
->
Ec_meam
[
i
][
j
]))
{
if
(
this
->
lattce_meam
[
i
][
j
]
==
L12
)
this
->
Ec_meam
[
i
][
j
]
=
(
3
*
this
->
Ec_meam
[
i
][
i
]
+
this
->
Ec_meam
[
j
][
j
])
/
4.0
-
this
->
delta_meam
[
i
][
j
];
else
if
(
this
->
lattce_meam
[
i
][
j
]
==
C11
)
{
if
(
this
->
lattce_meam
[
i
][
i
]
==
DIA
)
this
->
Ec_meam
[
i
][
j
]
=
(
2
*
this
->
Ec_meam
[
i
][
i
]
+
this
->
Ec_meam
[
j
][
j
])
/
3.0
-
this
->
delta_meam
[
i
][
j
];
else
this
->
Ec_meam
[
i
][
j
]
=
(
this
->
Ec_meam
[
i
][
i
]
+
2
*
this
->
Ec_meam
[
j
][
j
])
/
3.0
-
this
->
delta_meam
[
i
][
j
];
}
else
this
->
Ec_meam
[
i
][
j
]
=
(
this
->
Ec_meam
[
i
][
i
]
+
this
->
Ec_meam
[
j
][
j
])
/
2.0
-
this
->
delta_meam
[
i
][
j
];
}
if
(
iszero
(
this
->
alpha_meam
[
i
][
j
]))
this
->
alpha_meam
[
i
][
j
]
=
(
this
->
alpha_meam
[
i
][
i
]
+
this
->
alpha_meam
[
j
][
j
])
/
2.0
;
if
(
iszero
(
this
->
re_meam
[
i
][
j
]))
this
->
re_meam
[
i
][
j
]
=
(
this
->
re_meam
[
i
][
i
]
+
this
->
re_meam
[
j
][
j
])
/
2.0
;
}
}
}
// Cmin[i][k][j] is symmetric in i-j, but not k. For all triplets
// where i>j, set equal to the i<j element. Likewise for Cmax.
for
(
i
=
1
;
i
<
this
->
neltypes
;
i
++
)
{
for
(
j
=
0
;
j
<
i
;
j
++
)
{
for
(
k
=
0
;
k
<
this
->
neltypes
;
k
++
)
{
this
->
Cmin_meam
[
i
][
j
][
k
]
=
this
->
Cmin_meam
[
j
][
i
][
k
];
this
->
Cmax_meam
[
i
][
j
][
k
]
=
this
->
Cmax_meam
[
j
][
i
][
k
];
}
}
}
// ebound gives the squared distance such that, for rik2 or rjk2>ebound,
// atom k definitely lies outside the screening function ellipse (so
// there is no need to calculate its effects). Here, compute it for all
// triplets [i][j][k] so that ebound[i][j] is the maximized over k
for
(
i
=
0
;
i
<
this
->
neltypes
;
i
++
)
{
for
(
j
=
0
;
j
<
this
->
neltypes
;
j
++
)
{
for
(
k
=
0
;
k
<
this
->
neltypes
;
k
++
)
{
eb
=
(
this
->
Cmax_meam
[
i
][
j
][
k
]
*
this
->
Cmax_meam
[
i
][
j
][
k
])
/
(
4.0
*
(
this
->
Cmax_meam
[
i
][
j
][
k
]
-
1.0
));
this
->
ebound_meam
[
i
][
j
]
=
std
::
max
(
this
->
ebound_meam
[
i
][
j
],
eb
);
}
}
}
}
//-----------------------------------------------------------------------
// compute MEAM pair potential for each pair of element types
//
void
MEAM
::
compute_pair_meam
(
void
)
{
double
r
/*ununsed:, temp*/
;
int
j
,
a
,
b
,
nv2
;
double
astar
,
frac
,
phizbl
;
int
n
,
nmax
,
Z1
,
Z2
;
double
arat
,
rarat
,
scrn
,
scrn2
;
double
phiaa
,
phibb
/*unused:,phitmp*/
;
double
C
,
s111
,
s112
,
s221
,
S11
,
S22
;
// check for previously allocated arrays and free them
if
(
this
->
phir
!=
NULL
)
memory
->
destroy
(
this
->
phir
);
if
(
this
->
phirar
!=
NULL
)
memory
->
destroy
(
this
->
phirar
);
if
(
this
->
phirar1
!=
NULL
)
memory
->
destroy
(
this
->
phirar1
);
if
(
this
->
phirar2
!=
NULL
)
memory
->
destroy
(
this
->
phirar2
);
if
(
this
->
phirar3
!=
NULL
)
memory
->
destroy
(
this
->
phirar3
);
if
(
this
->
phirar4
!=
NULL
)
memory
->
destroy
(
this
->
phirar4
);
if
(
this
->
phirar5
!=
NULL
)
memory
->
destroy
(
this
->
phirar5
);
if
(
this
->
phirar6
!=
NULL
)
memory
->
destroy
(
this
->
phirar6
);
// allocate memory for array that defines the potential
memory
->
create
(
this
->
phir
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phir"
);
// allocate coeff memory
memory
->
create
(
this
->
phirar
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phirar"
);
memory
->
create
(
this
->
phirar1
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phirar1"
);
memory
->
create
(
this
->
phirar2
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phirar2"
);
memory
->
create
(
this
->
phirar3
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phirar3"
);
memory
->
create
(
this
->
phirar4
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phirar4"
);
memory
->
create
(
this
->
phirar5
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phirar5"
);
memory
->
create
(
this
->
phirar6
,
(
this
->
neltypes
*
(
this
->
neltypes
+
1
))
/
2
,
this
->
nr
,
"pair:phirar6"
);
// loop over pairs of element types
nv2
=
0
;
for
(
a
=
0
;
a
<
this
->
neltypes
;
a
++
)
{
for
(
b
=
a
;
b
<
this
->
neltypes
;
b
++
)
{
// loop over r values and compute
for
(
j
=
0
;
j
<
this
->
nr
;
j
++
)
{
r
=
j
*
this
->
dr
;
this
->
phir
[
nv2
][
j
]
=
phi_meam
(
r
,
a
,
b
);
// if using second-nearest neighbor, solve recursive problem
// (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
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
);
// 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
// pairs, but need to include NN2 contributions to those pairs as
// well.
if
(
this
->
lattce_meam
[
a
][
b
]
==
B1
||
this
->
lattce_meam
[
a
][
b
]
==
B2
||
this
->
lattce_meam
[
a
][
b
]
==
L12
)
{
rarat
=
r
*
arat
;
// phi_aa
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
);
nmax
=
10
;
if
(
scrn
>
0.0
)
{
for
(
n
=
1
;
n
<=
nmax
;
n
++
)
{
phiaa
=
phiaa
+
pow
((
-
Z2
*
scrn
/
Z1
),
n
)
*
phi_meam
(
rarat
*
pow
(
arat
,
n
),
a
,
a
);
}
}
// 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
);
nmax
=
10
;
if
(
scrn
>
0.0
)
{
for
(
n
=
1
;
n
<=
nmax
;
n
++
)
{
phibb
=
phibb
+
pow
((
-
Z2
*
scrn
/
Z1
),
n
)
*
phi_meam
(
rarat
*
pow
(
arat
,
n
),
b
,
b
);
}
}
if
(
this
->
lattce_meam
[
a
][
b
]
==
B1
||
this
->
lattce_meam
[
a
][
b
]
==
B2
)
{
// 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
->
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
->
phir
[
nv2
][
j
]
=
this
->
phir
[
nv2
][
j
]
-
Z2
*
scrn2
/
(
2
*
Z1
)
*
phibb
;
}
else
if
(
this
->
lattce_meam
[
a
][
b
]
==
L12
)
{
// The L12 case has one last trick; we have to be careful to
// compute
// the correct screening between 2nd-neighbor pairs. 1-1
// second-neighbor pairs are screened by 2 type 1 atoms and
// two type
// 2 atoms. 2-2 second-neighbor pairs are screened by 4 type
// 1
// atoms.
C
=
1.0
;
get_sijk
(
C
,
a
,
a
,
a
,
&
s111
);
get_sijk
(
C
,
a
,
a
,
b
,
&
s112
);
get_sijk
(
C
,
b
,
b
,
a
,
&
s221
);
S11
=
s111
*
s111
*
s112
*
s112
;
S22
=
pow
(
s221
,
4
);
this
->
phir
[
nv2
][
j
]
=
this
->
phir
[
nv2
][
j
]
-
0.75
*
S11
*
phiaa
-
0.25
*
S22
*
phibb
;
}
}
else
{
nmax
=
10
;
for
(
n
=
1
;
n
<=
nmax
;
n
++
)
{
this
->
phir
[
nv2
][
j
]
=
this
->
phir
[
nv2
][
j
]
+
pow
((
-
Z2
*
scrn
/
Z1
),
n
)
*
phi_meam
(
r
*
pow
(
arat
,
n
),
a
,
b
);
}
}
}
// For Zbl potential:
// if astar <= -3
// potential is zbl potential
// else if -3 < astar < -1
// potential is linear combination with zbl potential
// endif
if
(
this
->
zbl_meam
[
a
][
b
]
==
1
)
{
astar
=
this
->
alpha_meam
[
a
][
b
]
*
(
r
/
this
->
re_meam
[
a
][
b
]
-
1.0
);
if
(
astar
<=
-
3.0
)
this
->
phir
[
nv2
][
j
]
=
zbl
(
r
,
this
->
ielt_meam
[
a
],
this
->
ielt_meam
[
b
]);
else
if
(
astar
>
-
3.0
&&
astar
<
-
1.0
)
{
frac
=
fcut
(
1
-
(
astar
+
1.0
)
/
(
-
3.0
+
1.0
));
phizbl
=
zbl
(
r
,
this
->
ielt_meam
[
a
],
this
->
ielt_meam
[
b
]);
this
->
phir
[
nv2
][
j
]
=
frac
*
this
->
phir
[
nv2
][
j
]
+
(
1
-
frac
)
*
phizbl
;
}
}
}
// call interpolation
interpolate_meam
(
nv2
);
nv2
=
nv2
+
1
;
}
}
}
//----------------------------------------------------------------------c
// Compute MEAM pair potential for distance r, element types a and b
//
double
MEAM
::
phi_meam
(
double
r
,
int
a
,
int
b
)
{
/*unused:double a1,a2,a12;*/
double
t11av
,
t21av
,
t31av
,
t12av
,
t22av
,
t32av
;
double
G1
,
G2
,
s1
[
3
],
s2
[
3
],
rho0_1
,
rho0_2
;
double
Gam1
,
Gam2
,
Z1
,
Z2
;
double
rhobar1
,
rhobar2
,
F1
,
F2
;
double
rho01
,
rho11
,
rho21
,
rho31
;
double
rho02
,
rho12
,
rho22
,
rho32
;
double
scalfac
,
phiaa
,
phibb
;
double
Eu
;
double
arat
,
scrn
/*unused:,scrn2*/
;
int
Z12
,
errorflag
;
int
n
,
nmax
,
Z1nn
,
Z2nn
;
lattice_t
latta
/*unused:,lattb*/
;
double
rho_bkgd1
,
rho_bkgd2
;
double
phi_m
=
0.0
;
// Equation numbers below refer to:
// I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615
// get number of neighbors in the reference structure
// Nref[i][j] = # of i's neighbors of type j
Z12
=
get_Zij
(
this
->
lattce_meam
[
a
][
b
]);
get_densref
(
r
,
a
,
b
,
&
rho01
,
&
rho11
,
&
rho21
,
&
rho31
,
&
rho02
,
&
rho12
,
&
rho22
,
&
rho32
);
// if densities are too small, numerical problems may result; just return zero
if
(
rho01
<=
1e-14
&&
rho02
<=
1e-14
)
return
0.0
;
// calculate average weighting factors for the reference structure
if
(
this
->
lattce_meam
[
a
][
b
]
==
C11
)
{
if
(
this
->
ialloy
==
2
)
{
t11av
=
this
->
t1_meam
[
a
];
t12av
=
this
->
t1_meam
[
b
];
t21av
=
this
->
t2_meam
[
a
];
t22av
=
this
->
t2_meam
[
b
];
t31av
=
this
->
t3_meam
[
a
];
t32av
=
this
->
t3_meam
[
b
];
}
else
{
scalfac
=
1.0
/
(
rho01
+
rho02
);
t11av
=
scalfac
*
(
this
->
t1_meam
[
a
]
*
rho01
+
this
->
t1_meam
[
b
]
*
rho02
);
t12av
=
t11av
;
t21av
=
scalfac
*
(
this
->
t2_meam
[
a
]
*
rho01
+
this
->
t2_meam
[
b
]
*
rho02
);
t22av
=
t21av
;
t31av
=
scalfac
*
(
this
->
t3_meam
[
a
]
*
rho01
+
this
->
t3_meam
[
b
]
*
rho02
);
t32av
=
t31av
;
}
}
else
{
// average weighting factors for the reference structure, eqn. I.8
get_tavref
(
&
t11av
,
&
t21av
,
&
t31av
,
&
t12av
,
&
t22av
,
&
t32av
,
this
->
t1_meam
[
a
],
this
->
t2_meam
[
a
],
this
->
t3_meam
[
a
],
this
->
t1_meam
[
b
],
this
->
t2_meam
[
b
],
this
->
t3_meam
[
b
],
r
,
a
,
b
,
this
->
lattce_meam
[
a
][
b
]);
}
// for c11b structure, calculate background electron densities
if
(
this
->
lattce_meam
[
a
][
b
]
==
C11
)
{
latta
=
this
->
lattce_meam
[
a
][
a
];
if
(
latta
==
DIA
)
{
rhobar1
=
pow
(((
Z12
/
2
)
*
(
rho02
+
rho01
)),
2
)
+
t11av
*
pow
((
rho12
-
rho11
),
2
)
+
t21av
/
6.0
*
pow
(
rho22
+
rho21
,
2
)
+
121.0
/
40.0
*
t31av
*
pow
((
rho32
-
rho31
),
2
);
rhobar1
=
sqrt
(
rhobar1
);
rhobar2
=
pow
(
Z12
*
rho01
,
2
)
+
2.0
/
3.0
*
t21av
*
pow
(
rho21
,
2
);
rhobar2
=
sqrt
(
rhobar2
);
}
else
{
rhobar2
=
pow
(((
Z12
/
2
)
*
(
rho01
+
rho02
)),
2
)
+
t12av
*
pow
((
rho11
-
rho12
),
2
)
+
t22av
/
6.0
*
pow
(
rho21
+
rho22
,
2
)
+
121.0
/
40.0
*
t32av
*
pow
((
rho31
-
rho32
),
2
);
rhobar2
=
sqrt
(
rhobar2
);
rhobar1
=
pow
(
Z12
*
rho02
,
2
)
+
2.0
/
3.0
*
t22av
*
pow
(
rho22
,
2
);
rhobar1
=
sqrt
(
rhobar1
);
}
}
else
{
// for other structures, use formalism developed in Huang's paper
//
// composition-dependent scaling, equation I.7
// If using mixing rule for t, apply to reference structure; else
// use precomputed values
if
(
this
->
mix_ref_t
==
1
)
{
Z1
=
this
->
Z_meam
[
a
];
Z2
=
this
->
Z_meam
[
b
];
if
(
this
->
ibar_meam
[
a
]
<=
0
)
G1
=
1.0
;
else
{
get_shpfcn
(
this
->
lattce_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
);
Gam2
=
(
s2
[
0
]
*
t12av
+
s2
[
1
]
*
t22av
+
s2
[
2
]
*
t32av
)
/
(
Z2
*
Z2
);
G2
=
G_gam
(
Gam2
,
this
->
ibar_meam
[
b
],
errorflag
);
}
rho0_1
=
this
->
rho0_meam
[
a
]
*
Z1
*
G1
;
rho0_2
=
this
->
rho0_meam
[
b
]
*
Z2
*
G2
;
}
Gam1
=
(
t11av
*
rho11
+
t21av
*
rho21
+
t31av
*
rho31
);
if
(
rho01
<
1.0e-14
)
Gam1
=
0.0
;
else
Gam1
=
Gam1
/
(
rho01
*
rho01
);
Gam2
=
(
t12av
*
rho12
+
t22av
*
rho22
+
t32av
*
rho32
);
if
(
rho02
<
1.0e-14
)
Gam2
=
0.0
;
else
Gam2
=
Gam2
/
(
rho02
*
rho02
);
G1
=
G_gam
(
Gam1
,
this
->
ibar_meam
[
a
],
errorflag
);
G2
=
G_gam
(
Gam2
,
this
->
ibar_meam
[
b
],
errorflag
);
if
(
this
->
mix_ref_t
==
1
)
{
rho_bkgd1
=
rho0_1
;
rho_bkgd2
=
rho0_2
;
}
else
{
if
(
this
->
bkgd_dyn
==
1
)
{
rho_bkgd1
=
this
->
rho0_meam
[
a
]
*
this
->
Z_meam
[
a
];
rho_bkgd2
=
this
->
rho0_meam
[
b
]
*
this
->
Z_meam
[
b
];
}
else
{
rho_bkgd1
=
this
->
rho_ref_meam
[
a
];
rho_bkgd2
=
this
->
rho_ref_meam
[
b
];
}
}
rhobar1
=
rho01
/
rho_bkgd1
*
G1
;
rhobar2
=
rho02
/
rho_bkgd2
*
G2
;
}
// compute embedding functions, eqn I.5
if
(
iszero
(
rhobar1
))
F1
=
0.0
;
else
{
if
(
this
->
emb_lin_neg
==
1
&&
rhobar1
<=
0
)
F1
=
-
this
->
A_meam
[
a
]
*
this
->
Ec_meam
[
a
][
a
]
*
rhobar1
;
else
F1
=
this
->
A_meam
[
a
]
*
this
->
Ec_meam
[
a
][
a
]
*
rhobar1
*
log
(
rhobar1
);
}
if
(
iszero
(
rhobar2
))
F2
=
0.0
;
else
{
if
(
this
->
emb_lin_neg
==
1
&&
rhobar2
<=
0
)
F2
=
-
this
->
A_meam
[
b
]
*
this
->
Ec_meam
[
b
][
b
]
*
rhobar2
;
else
F2
=
this
->
A_meam
[
b
]
*
this
->
Ec_meam
[
b
][
b
]
*
rhobar2
*
log
(
rhobar2
);
}
// compute Rose function, I.16
Eu
=
erose
(
r
,
this
->
re_meam
[
a
][
b
],
this
->
alpha_meam
[
a
][
b
],
this
->
Ec_meam
[
a
][
b
],
this
->
repuls_meam
[
a
][
b
],
this
->
attrac_meam
[
a
][
b
],
this
->
erose_form
);
// calculate the pair energy
if
(
this
->
lattce_meam
[
a
][
b
]
==
C11
)
{
latta
=
this
->
lattce_meam
[
a
][
a
];
if
(
latta
==
DIA
)
{
phiaa
=
phi_meam
(
r
,
a
,
a
);
phi_m
=
(
3
*
Eu
-
F2
-
2
*
F1
-
5
*
phiaa
)
/
Z12
;
}
else
{
phibb
=
phi_meam
(
r
,
b
,
b
);
phi_m
=
(
3
*
Eu
-
F1
-
2
*
F2
-
5
*
phibb
)
/
Z12
;
}
}
else
if
(
this
->
lattce_meam
[
a
][
b
]
==
L12
)
{
phiaa
=
phi_meam
(
r
,
a
,
a
);
// account for second neighbor a-a potential here...
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
],
arat
,
scrn
);
nmax
=
10
;
if
(
scrn
>
0.0
)
{
for
(
n
=
1
;
n
<=
nmax
;
n
++
)
{
phiaa
=
phiaa
+
pow
((
-
Z2nn
*
scrn
/
Z1nn
),
n
)
*
phi_meam
(
r
*
pow
(
arat
,
n
),
a
,
a
);
}
}
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
if
(
iszero
(
r
))
{
phi_m
=
0.0
;
}
return
phi_m
;
}
//----------------------------------------------------------------------c
// Compute background density for reference structure of each element
void
MEAM
::
compute_reference_density
(
void
)
{
int
a
,
Z
,
Z2
,
errorflag
;
double
gam
,
Gbar
,
shp
[
3
];
double
rho0
,
rho0_2nn
,
arat
,
scrn
;
// loop over element types
for
(
a
=
0
;
a
<
this
->
neltypes
;
a
++
)
{
Z
=
(
int
)
this
->
Z_meam
[
a
];
if
(
this
->
ibar_meam
[
a
]
<=
0
)
Gbar
=
1.0
;
else
{
get_shpfcn
(
this
->
lattce_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
);
}
// The zeroth order density in the reference structure, with
// equilibrium spacing, is just the number of first neighbors times
// the rho0_meam coefficient...
rho0
=
this
->
rho0_meam
[
a
]
*
Z
;
// ...unless we have unscreened second neighbors, in which case we
// add on the contribution from those (accounting for partial
// 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
);
rho0_2nn
=
this
->
rho0_meam
[
a
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta0_meam
[
a
]
*
(
arat
-
1
));
rho0
=
rho0
+
Z2
*
rho0_2nn
*
scrn
;
}
this
->
rho_ref_meam
[
a
]
=
rho0
*
Gbar
;
}
}
//------------------------------------------------------------------------------c
// Average weighting factors for the reference structure
void
MEAM
::
get_tavref
(
double
*
t11av
,
double
*
t21av
,
double
*
t31av
,
double
*
t12av
,
double
*
t22av
,
double
*
t32av
,
double
t11
,
double
t21
,
double
t31
,
double
t12
,
double
t22
,
double
t32
,
double
r
,
int
a
,
int
b
,
lattice_t
latt
)
{
double
rhoa01
,
rhoa02
,
a1
,
a2
,
rho01
/*,rho02*/
;
// For ialloy = 2, no averaging is done
if
(
this
->
ialloy
==
2
)
{
*
t11av
=
t11
;
*
t21av
=
t21
;
*
t31av
=
t31
;
*
t12av
=
t12
;
*
t22av
=
t22
;
*
t32av
=
t32
;
}
else
{
if
(
latt
==
FCC
||
latt
==
BCC
||
latt
==
DIA
||
latt
==
HCP
||
latt
==
B1
||
latt
==
DIM
||
latt
==
B2
)
{
// all neighbors are of the opposite type
*
t11av
=
t12
;
*
t21av
=
t22
;
*
t31av
=
t32
;
*
t12av
=
t11
;
*
t22av
=
t21
;
*
t32av
=
t31
;
}
else
{
a1
=
r
/
this
->
re_meam
[
a
][
a
]
-
1.0
;
a2
=
r
/
this
->
re_meam
[
b
][
b
]
-
1.0
;
rhoa01
=
this
->
rho0_meam
[
a
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta0_meam
[
a
]
*
a1
);
rhoa02
=
this
->
rho0_meam
[
b
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta0_meam
[
b
]
*
a2
);
if
(
latt
==
L12
)
{
rho01
=
8
*
rhoa01
+
4
*
rhoa02
;
*
t11av
=
(
8
*
t11
*
rhoa01
+
4
*
t12
*
rhoa02
)
/
rho01
;
*
t12av
=
t11
;
*
t21av
=
(
8
*
t21
*
rhoa01
+
4
*
t22
*
rhoa02
)
/
rho01
;
*
t22av
=
t21
;
*
t31av
=
(
8
*
t31
*
rhoa01
+
4
*
t32
*
rhoa02
)
/
rho01
;
*
t32av
=
t31
;
}
else
{
// call error('Lattice not defined in get_tavref.')
}
}
}
}
//------------------------------------------------------------------------------c
void
MEAM
::
get_sijk
(
double
C
,
int
i
,
int
j
,
int
k
,
double
*
sijk
)
{
double
x
;
x
=
(
C
-
this
->
Cmin_meam
[
i
][
j
][
k
])
/
(
this
->
Cmax_meam
[
i
][
j
][
k
]
-
this
->
Cmin_meam
[
i
][
j
][
k
]);
*
sijk
=
fcut
(
x
);
}
//------------------------------------------------------------------------------c
// Calculate density functions, assuming reference configuration
void
MEAM
::
get_densref
(
double
r
,
int
a
,
int
b
,
double
*
rho01
,
double
*
rho11
,
double
*
rho21
,
double
*
rho31
,
double
*
rho02
,
double
*
rho12
,
double
*
rho22
,
double
*
rho32
)
{
double
a1
,
a2
;
double
s
[
3
];
lattice_t
lat
;
int
Zij2nn
;
double
rhoa01nn
,
rhoa02nn
;
double
rhoa01
,
rhoa11
,
rhoa21
,
rhoa31
;
double
rhoa02
,
rhoa12
,
rhoa22
,
rhoa32
;
double
arat
,
scrn
,
denom
;
double
C
,
s111
,
s112
,
s221
,
S11
,
S22
;
a1
=
r
/
this
->
re_meam
[
a
][
a
]
-
1.0
;
a2
=
r
/
this
->
re_meam
[
b
][
b
]
-
1.0
;
rhoa01
=
this
->
rho0_meam
[
a
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta0_meam
[
a
]
*
a1
);
rhoa11
=
this
->
rho0_meam
[
a
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta1_meam
[
a
]
*
a1
);
rhoa21
=
this
->
rho0_meam
[
a
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta2_meam
[
a
]
*
a1
);
rhoa31
=
this
->
rho0_meam
[
a
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta3_meam
[
a
]
*
a1
);
rhoa02
=
this
->
rho0_meam
[
b
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta0_meam
[
b
]
*
a2
);
rhoa12
=
this
->
rho0_meam
[
b
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta1_meam
[
b
]
*
a2
);
rhoa22
=
this
->
rho0_meam
[
b
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta2_meam
[
b
]
*
a2
);
rhoa32
=
this
->
rho0_meam
[
b
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta3_meam
[
b
]
*
a2
);
lat
=
this
->
lattce_meam
[
a
][
b
];
*
rho11
=
0.0
;
*
rho21
=
0.0
;
*
rho31
=
0.0
;
*
rho12
=
0.0
;
*
rho22
=
0.0
;
*
rho32
=
0.0
;
if
(
lat
==
FCC
)
{
*
rho01
=
12.0
*
rhoa02
;
*
rho02
=
12.0
*
rhoa01
;
}
else
if
(
lat
==
BCC
)
{
*
rho01
=
8.0
*
rhoa02
;
*
rho02
=
8.0
*
rhoa01
;
}
else
if
(
lat
==
B1
)
{
*
rho01
=
6.0
*
rhoa02
;
*
rho02
=
6.0
*
rhoa01
;
}
else
if
(
lat
==
DIA
)
{
*
rho01
=
4.0
*
rhoa02
;
*
rho02
=
4.0
*
rhoa01
;
*
rho31
=
32.0
/
9.0
*
rhoa32
*
rhoa32
;
*
rho32
=
32.0
/
9.0
*
rhoa31
*
rhoa31
;
}
else
if
(
lat
==
HCP
)
{
*
rho01
=
12
*
rhoa02
;
*
rho02
=
12
*
rhoa01
;
*
rho31
=
1.0
/
3.0
*
rhoa32
*
rhoa32
;
*
rho32
=
1.0
/
3.0
*
rhoa31
*
rhoa31
;
}
else
if
(
lat
==
DIM
)
{
get_shpfcn
(
DIM
,
s
);
*
rho01
=
rhoa02
;
*
rho02
=
rhoa01
;
*
rho11
=
s
[
0
]
*
rhoa12
*
rhoa12
;
*
rho12
=
s
[
0
]
*
rhoa11
*
rhoa11
;
*
rho21
=
s
[
1
]
*
rhoa22
*
rhoa22
;
*
rho22
=
s
[
1
]
*
rhoa21
*
rhoa21
;
*
rho31
=
s
[
2
]
*
rhoa32
*
rhoa32
;
*
rho32
=
s
[
2
]
*
rhoa31
*
rhoa31
;
}
else
if
(
lat
==
C11
)
{
*
rho01
=
rhoa01
;
*
rho02
=
rhoa02
;
*
rho11
=
rhoa11
;
*
rho12
=
rhoa12
;
*
rho21
=
rhoa21
;
*
rho22
=
rhoa22
;
*
rho31
=
rhoa31
;
*
rho32
=
rhoa32
;
}
else
if
(
lat
==
L12
)
{
*
rho01
=
8
*
rhoa01
+
4
*
rhoa02
;
*
rho02
=
12
*
rhoa01
;
if
(
this
->
ialloy
==
1
)
{
*
rho21
=
8.
/
3.
*
pow
(
rhoa21
*
this
->
t2_meam
[
a
]
-
rhoa22
*
this
->
t2_meam
[
b
],
2
);
denom
=
8
*
rhoa01
*
pow
(
this
->
t2_meam
[
a
],
2
)
+
4
*
rhoa02
*
pow
(
this
->
t2_meam
[
b
],
2
);
if
(
denom
>
0.
)
*
rho21
=
*
rho21
/
denom
*
*
rho01
;
}
else
*
rho21
=
8.
/
3.
*
(
rhoa21
-
rhoa22
)
*
(
rhoa21
-
rhoa22
);
}
else
if
(
lat
==
B2
)
{
*
rho01
=
8.0
*
rhoa02
;
*
rho02
=
8.0
*
rhoa01
;
}
else
{
// 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
);
a1
=
arat
*
r
/
this
->
re_meam
[
a
][
a
]
-
1.0
;
a2
=
arat
*
r
/
this
->
re_meam
[
b
][
b
]
-
1.0
;
rhoa01nn
=
this
->
rho0_meam
[
a
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta0_meam
[
a
]
*
a1
);
rhoa02nn
=
this
->
rho0_meam
[
b
]
*
MathSpecial
::
fm_exp
(
-
this
->
beta0_meam
[
b
]
*
a2
);
if
(
lat
==
L12
)
{
// As usual, L12 thinks it's special; we need to be careful computing
// the screening functions
C
=
1.0
;
get_sijk
(
C
,
a
,
a
,
a
,
&
s111
);
get_sijk
(
C
,
a
,
a
,
b
,
&
s112
);
get_sijk
(
C
,
b
,
b
,
a
,
&
s221
);
S11
=
s111
*
s111
*
s112
*
s112
;
S22
=
pow
(
s221
,
4
);
*
rho01
=
*
rho01
+
6
*
S11
*
rhoa01nn
;
*
rho02
=
*
rho02
+
6
*
S22
*
rhoa02nn
;
}
else
{
// For other cases, assume that second neighbor is of same type,
// first neighbor may be of different type
*
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
);
*
rho02
=
*
rho02
+
Zij2nn
*
scrn
*
rhoa02nn
;
}
}
}
void
MEAM
::
interpolate_meam
(
int
ind
)
{
int
j
;
double
drar
;
// map to coefficient space
this
->
nrar
=
this
->
nr
;
drar
=
this
->
dr
;
this
->
rdrar
=
1.0
/
drar
;
// phir interp
for
(
j
=
0
;
j
<
this
->
nrar
;
j
++
)
{
this
->
phirar
[
ind
][
j
]
=
this
->
phir
[
ind
][
j
];
}
this
->
phirar1
[
ind
][
0
]
=
this
->
phirar
[
ind
][
1
]
-
this
->
phirar
[
ind
][
0
];
this
->
phirar1
[
ind
][
1
]
=
0.5
*
(
this
->
phirar
[
ind
][
2
]
-
this
->
phirar
[
ind
][
0
]);
this
->
phirar1
[
ind
][
this
->
nrar
-
2
]
=
0.5
*
(
this
->
phirar
[
ind
][
this
->
nrar
-
1
]
-
this
->
phirar
[
ind
][
this
->
nrar
-
3
]);
this
->
phirar1
[
ind
][
this
->
nrar
-
1
]
=
0.0
;
for
(
j
=
2
;
j
<
this
->
nrar
-
2
;
j
++
)
{
this
->
phirar1
[
ind
][
j
]
=
((
this
->
phirar
[
ind
][
j
-
2
]
-
this
->
phirar
[
ind
][
j
+
2
])
+
8.0
*
(
this
->
phirar
[
ind
][
j
+
1
]
-
this
->
phirar
[
ind
][
j
-
1
]))
/
12.
;
}
for
(
j
=
0
;
j
<
this
->
nrar
-
1
;
j
++
)
{
this
->
phirar2
[
ind
][
j
]
=
3.0
*
(
this
->
phirar
[
ind
][
j
+
1
]
-
this
->
phirar
[
ind
][
j
])
-
2.0
*
this
->
phirar1
[
ind
][
j
]
-
this
->
phirar1
[
ind
][
j
+
1
];
this
->
phirar3
[
ind
][
j
]
=
this
->
phirar1
[
ind
][
j
]
+
this
->
phirar1
[
ind
][
j
+
1
]
-
2.0
*
(
this
->
phirar
[
ind
][
j
+
1
]
-
this
->
phirar
[
ind
][
j
]);
}
this
->
phirar2
[
ind
][
this
->
nrar
-
1
]
=
0.0
;
this
->
phirar3
[
ind
][
this
->
nrar
-
1
]
=
0.0
;
for
(
j
=
0
;
j
<
this
->
nrar
;
j
++
)
{
this
->
phirar4
[
ind
][
j
]
=
this
->
phirar1
[
ind
][
j
]
/
drar
;
this
->
phirar5
[
ind
][
j
]
=
2.0
*
this
->
phirar2
[
ind
][
j
]
/
drar
;
this
->
phirar6
[
ind
][
j
]
=
3.0
*
this
->
phirar3
[
ind
][
j
]
/
drar
;
}
}
//---------------------------------------------------------------------
// Compute Rose energy function, I.16
//
double
MEAM
::
compute_phi
(
double
rij
,
int
elti
,
int
eltj
)
{
double
pp
;
int
ind
,
kk
;
ind
=
this
->
eltind
[
elti
][
eltj
];
pp
=
rij
*
this
->
rdrar
;
kk
=
(
int
)
pp
;
kk
=
std
::
min
(
kk
,
this
->
nrar
-
2
);
pp
=
pp
-
kk
;
pp
=
std
::
min
(
pp
,
1.0
);
double
result
=
((
this
->
phirar3
[
ind
][
kk
]
*
pp
+
this
->
phirar2
[
ind
][
kk
])
*
pp
+
this
->
phirar1
[
ind
][
kk
])
*
pp
+
this
->
phirar
[
ind
][
kk
];
return
result
;
}
Event Timeline
Log In to Comment