jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Style changes

This commit is contained in:
rohskopf 2023-01-25 14:52:02 -07:00
parent f58aeecec0
commit d88a4a768d
18 changed files with 135 additions and 240 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

0
examples/meam/msmeam/HGaMS.meam → examples/meam/msmeam/HGa.meam
View File

2
examples/meam/msmeam/library.MSmeam → examples/meam/msmeam/library.msmeam
View File

@ -1,4 +1,4 @@
# DATE: 2018-09-22 MB DATE: 2016-09-21 CONTRIBUTOR: Steve Valone, smv@lanl.gov CITATION: Baskes, PRB 1992; smv, sr, mib, JNM 2010 # DATE: 2018-09-22 UNITS: metal CONTRIBUTOR: Steve Valone, smv@lanl.gov CITATION: Baskes, PRB 1992; smv, sr, mib, JNM 2010
# ms-meam data format May 2010 # ms-meam data format May 2010
# elt lat z ielement atwt # elt lat z ielement atwt
# alpha b0 b1 b2 b3 b1m b2m b3m alat esub asub # alpha b0 b1 b2 b3 b1m b2m b3m alat esub asub

2
examples/meam/msmeam/potential.mod
View File

@ -7,7 +7,7 @@ print "potential chosen ${Pu}"
pair_style meam ms # the `ms` flag turns on MS-MEAM pair_style meam ms # the `ms` flag turns on MS-MEAM
print "we just executed" print "we just executed"
pair_coeff * * library.MSmeam ${Pu} Ga4 HGaMS.meam ${Pu} Ga4 pair_coeff * * library.msmeam ${Pu} Ga4 HGa.meam ${Pu} Ga4
# Setup neighbor style # Setup neighbor style
neighbor 1.0 bin neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes neigh_modify once no every 1 delay 0 check yes

2
examples/meam/msmeam/test.data
View File

@ -12,7 +12,7 @@ Masses
1 1.0079 1 1.0079
2 69.723 2 69.723
Atoms # atomic Atoms # atomic/kk
1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0
2 2 2.2 0 0 0 0 0 2 2 2.2 0 0 0 0 0

24
examples/meam/msmeam/test.dump
View File

@ -1,24 +0,0 @@
ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
3
ITEM: BOX BOUNDS pp pp pp
-4.0000000000000000e+00 4.0000000000000000e+00
-4.0000000000000000e+00 4.0000000000000000e+00
-4.0000000000000000e+00 4.0000000000000000e+00
ITEM: ATOMS id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
1 0 0 0 -131.925 -88.3005 0 22.9153 -2.147e+08 -1.62661e+08 -0 -2.05301e+07 -0 -0
2 2.2 0 0 120.809 -0.482171 0 14.7692 -2.12028e+08 -0 -0 403352 -0 -0
3 0.3 2.3 0 11.1159 88.7827 0 8.61478 -2.67145e+06 -1.62661e+08 -0 -2.09335e+07 -0 -0
ITEM: TIMESTEP
1
ITEM: NUMBER OF ATOMS
3
ITEM: BOX BOUNDS pp pp pp
-4.0000000000000000e+00 4.0000000000000000e+00
-4.0000000000000000e+00 4.0000000000000000e+00
-4.0000000000000000e+00 4.0000000000000000e+00
ITEM: ATOMS id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
1 0 0 0 -131.925 -88.3005 0 22.9153 -2.147e+08 -1.62661e+08 -0 -2.05301e+07 -0 -0
2 2.2 0 0 120.809 -0.482171 0 14.7692 -2.12028e+08 -0 -0 403352 -0 -0
3 0.3 2.3 0 11.1159 88.7827 0 8.61478 -2.67145e+06 -1.62661e+08 -0 -2.09335e+07 -0 -0

123
examples/meam/msmeam/test.log
View File

@ -1,123 +0,0 @@
# Test of MEAM potential for HGa
# ------------------------ INITIALIZATION ----------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
variable latparam equal 4.646
variable ncell equal 3
# ----------------------- ATOM DEFINITION ----------------------------
region box block -4 4 -4 4 -4 4
create_box 2 box
Created orthogonal box = (-4 -4 -4) to (4 4 4)
1 by 1 by 1 MPI processor grid
#
include potential.mod
# NOTE: This script can be modified for different pair styles
# See in.elastic for more info.
variable Pu string H
print "potential chosen ${Pu}"
potential chosen H
# Choose potential
pair_style meam ms # the `ms` flag turns on MS-MEAM
print "we just executed"
we just executed
pair_coeff * * library.MSmeam ${Pu} Ga4 HGaMS.meam ${Pu} Ga4
pair_coeff * * library.MSmeam H Ga4 HGaMS.meam ${Pu} Ga4
pair_coeff * * library.MSmeam H Ga4 HGaMS.meam H Ga4
Reading MEAM library file library.MSmeam with DATE: 2018-09-22
# Setup neighbor style
neighbor 1.0 nsq
neigh_modify once no every 1 delay 0 check yes
# Setup minimization style
variable dmax equal 1.0e-2
min_style cg
min_modify dmax ${dmax} line quadratic
min_modify dmax 0.01 line quadratic
compute eng all pe/atom
compute eatoms all reduce sum c_eng
# Setup output
thermo 100
thermo_style custom step temp etotal press pxx pyy pzz pxy pxz pyz lx ly lz vol c_eatoms
thermo_modify norm yes
create_atoms 1 single 0 0 0 units box
Created 1 atoms
using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
create_atoms CPU = 0.000 seconds
create_atoms 2 single 2.2 0 0 units box
Created 1 atoms
using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
create_atoms CPU = 0.000 seconds
create_atoms 2 single 0.3 2.3 0 units box
Created 1 atoms
using box units in orthogonal box = (-4 -4 -4) to (4 4 4)
create_atoms CPU = 0.000 seconds
# ---------- Define Settings ---------------------
variable teng equal "c_eatoms"
compute pot_energy all pe/atom
compute stress all stress/atom NULL
dump 1 all custom 1 test.dump id x y z fx fy fz c_pot_energy c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6]
run 1
WARNING: No fixes with time integration, atoms won't move (src/verlet.cpp:60)
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 6.9
ghost atom cutoff = 6.9
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair meam, perpetual
attributes: full, newton on
pair build: full/nsq
stencil: none
bin: none
(2) pair meam, perpetual, half/full from (1)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
Per MPI rank memory allocation (min/avg/max) = 10.69 | 10.69 | 10.69 Mbytes
Step Temp TotEng Press Pxx Pyy Pzz Pxy Pxz Pyz Lx Ly Lz Volume c_eatoms
0 0 15.433079 491354.7 838670.96 635393.15 0 80195.797 0 0 8 8 8 512 15.433079
1 0 15.433079 491354.7 838670.96 635393.15 0 80195.797 0 0 8 8 8 512 15.433079
Loop time of 0.000110387 on 1 procs for 1 steps with 3 atoms
Performance: 782.701 ns/day, 0.031 hours/ns, 9059.038 timesteps/s, 27.177 katom-step/s
56.2% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.6411e-05 | 2.6411e-05 | 2.6411e-05 | 0.0 | 23.93
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 4.455e-06 | 4.455e-06 | 4.455e-06 | 0.0 | 4.04
Output | 7.4411e-05 | 7.4411e-05 | 7.4411e-05 | 0.0 | 67.41
Modify | 7.54e-07 | 7.54e-07 | 7.54e-07 | 0.0 | 0.68
Other | | 4.356e-06 | | | 3.95
Nlocal: 3 ave 3 max 3 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 78 ave 78 max 78 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 7 ave 7 max 7 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 14 ave 14 max 14 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 14
Ave neighs/atom = 4.6666667
Neighbor list builds = 0
Dangerous builds = 0
write_data test.data
System init for write_data ...
print "All done!"
All done!
Total wall time: 0:00:00

30
potentials/HGa.msmeam Normal file
View File

@ -0,0 +1,30 @@
bkgd_dyn = 1
emb_lin_neg = 1
augt1=0
ialloy=1
rc = 5.9
#H
attrac(1,1)=0.460
repuls(1,1)=0.460
Cmin(1,1,1)=1.3 # PuMS
Cmax(1,1,1)= 2.80
nn2(1,1)=1
#Ga
rho0(2) = 0.6
attrac(2,2)=0.097
repuls(2,2)=0.097
nn2(2,2)=1
#HGa
attrac(1,2)=0.300
repuls(1,2)=0.300
lattce(1,2)=l12
re(1,2)=3.19
delta(1,2)=-0.48
alpha(1,2)=6.6
Cmin(1,1,2)=2.0
Cmin(2,1,2)= 2.0
Cmin(1,2,1)=2.0
Cmin(2,2,1) = 1.4
Cmin(1,2,2) = 1.4
Cmin(1,1,2) = 1.4
nn2(1,2)=1

14
potentials/library.msmeam Normal file
View File

@ -0,0 +1,14 @@
# DATE: 2018-09-22 UNITS: metal CONTRIBUTOR: Steve Valone, smv@lanl.gov CITATION: Baskes, PRB 1992; smv, sr, mib, JNM 2010
# ms-meam data format May 2010
# elt lat z ielement atwt
# alpha b0 b1 b2 b3 b1m b2m b3m alat esub asub
# - t0 t1 t2 t3 t1m t2m t3m rozero ibar
# NOTE: leading character cannot be a space
'H' 'dim' 1.0 1 1.0079
2.960 2.960 3.0 1.0 1.0 1.0 3.0 1.0 0.741 2.235 2.50
1.0 0.44721 0.0 0.00 0.0 0.31623 0 6.70 0
'Ga4' 'fcc' 12.0 31 69.723
4.42 4.80 3.10 6.00 0.00 0.0 0.0 0.5 4.247 2.897 0.97
1.0 1.649 1.435 0.00 0.0 0.0 2.0 0.70 0

10
src/KOKKOS/meam_dens_final_kokkos.h
View File

@ -61,7 +61,7 @@ void MEAMKokkos<DeviceType>::operator()(TagMEAMDensFinal, const int &i, EV_FLOAT
if (elti >= 0) { if (elti >= 0) {
scaleii = d_scale(type[i],type[i]); scaleii = d_scale(type[i],type[i]);
d_rho1[i] = 0.0; d_rho1[i] = 0.0;
if (msmeamflag){ if (msmeamflag) {
d_rho2[i] = -1.0 / 3.0 * (d_arho2b[i] * d_arho2b[i] d_rho2[i] = -1.0 / 3.0 * (d_arho2b[i] * d_arho2b[i]
- d_arho2mb[i] * d_arho2mb[i]); - d_arho2mb[i] * d_arho2mb[i]);
} else{ } else{
@ -69,7 +69,7 @@ void MEAMKokkos<DeviceType>::operator()(TagMEAMDensFinal, const int &i, EV_FLOAT
} }
d_rho3[i] = 0.0; d_rho3[i] = 0.0;
for (int m = 0; m < 3; m++) { for (int m = 0; m < 3; m++) {
if (msmeamflag){ if (msmeamflag) {
d_rho1[i] = d_rho1[i] + d_arho1(i, m) * d_arho1(i, m) d_rho1[i] = d_rho1[i] + d_arho1(i, m) * d_arho1(i, m)
- d_arho1m(i, m) * d_arho1m(i, m); - d_arho1m(i, m) * d_arho1m(i, m);
d_rho3[i] = d_rho3[i] - 3.0 / 5.0 * (d_arho3b(i, m) * d_arho3b(i, m) d_rho3[i] = d_rho3[i] - 3.0 / 5.0 * (d_arho3b(i, m) * d_arho3b(i, m)
@ -80,7 +80,7 @@ void MEAMKokkos<DeviceType>::operator()(TagMEAMDensFinal, const int &i, EV_FLOAT
} }
} }
for (int m = 0; m < 6; m++){ for (int m = 0; m < 6; m++){
if (msmeamflag){ if (msmeamflag) {
d_rho2[i] = d_rho2[i] + v2D[m] * (d_arho2(i, m) * d_arho2(i, m) d_rho2[i] = d_rho2[i] + v2D[m] * (d_arho2(i, m) * d_arho2(i, m)
- d_arho2m(i, m) * d_arho2m(i, m)); - d_arho2m(i, m) * d_arho2m(i, m));
} else{ } else{
@ -88,14 +88,14 @@ void MEAMKokkos<DeviceType>::operator()(TagMEAMDensFinal, const int &i, EV_FLOAT
} }
} }
for (int m = 0; m < 10; m++) for (int m = 0; m < 10; m++)
if (msmeamflag){ if (msmeamflag) {
d_rho3[i] = d_rho3[i] + v3D[m] * (d_arho3(i, m) * d_arho3(i, m) d_rho3[i] = d_rho3[i] + v3D[m] * (d_arho3(i, m) * d_arho3(i, m)
- d_arho3m(i, m) * d_arho3m(i, m)); - d_arho3m(i, m) * d_arho3m(i, m));
} else{ } else{
d_rho3[i] += v3D[m] * d_arho3(i,m) * d_arho3(i,m); d_rho3[i] += v3D[m] * d_arho3(i,m) * d_arho3(i,m);
} }
if (msmeamflag){ if (msmeamflag) {
// with msmeam all t weights are already accounted for in rho // with msmeam all t weights are already accounted for in rho
d_gamma[i] = d_rho1[i] + d_rho2[i] + d_rho3[i]; d_gamma[i] = d_rho1[i] + d_rho2[i] + d_rho3[i];
} else{ } else{

16
src/KOKKOS/meam_dens_init_kokkos.h
View File

@ -43,16 +43,16 @@ void MEAMKokkos<DeviceType>::operator()(TagMEAMZero, const int &i) const {
d_rho0[i] = 0.0; d_rho0[i] = 0.0;
d_arho2b[i] = 0.0; d_arho2b[i] = 0.0;
d_arho1(i,0) = d_arho1(i,1) = d_arho1(i,2) = 0.0; d_arho1(i,0) = d_arho1(i,1) = d_arho1(i,2) = 0.0;
if (this->msmeamflag){ if (this->msmeamflag) {
d_arho2mb[i] = 0.0; d_arho2mb[i] = 0.0;
d_arho1m(i,0) = d_arho1m(i,1) = d_arho1m(i,2) = 0.0; d_arho1m(i,0) = d_arho1m(i,1) = d_arho1m(i,2) = 0.0;
} }
for (int j = 0; j < 6; j++){ for (int j = 0; j < 6; j++) {
d_arho2(i,j) = 0.0; d_arho2(i,j) = 0.0;
if (msmeamflag) if (msmeamflag)
d_arho2m(i,j) = 0.0; d_arho2m(i,j) = 0.0;
} }
for (int j = 0; j < 10; j++){ for (int j = 0; j < 10; j++) {
d_arho3(i,j) = 0.0; d_arho3(i,j) = 0.0;
if (msmeamflag) if (msmeamflag)
d_arho3m(i,j) = 0.0; d_arho3m(i,j) = 0.0;
@ -540,7 +540,7 @@ MEAMKokkos<DeviceType>::calc_rho1(int i, int /*ntype*/, typename AT::t_int_1d ty
double rhoa3i = ro0i * MathSpecialKokkos::fm_exp(-beta3_meam[elti] * ai) * sij; double rhoa3i = ro0i * MathSpecialKokkos::fm_exp(-beta3_meam[elti] * ai) * sij;
// msmeam // msmeam
double rhoa1mj, rhoa2mj, rhoa3mj, rhoa1mi, rhoa2mi, rhoa3mi; double rhoa1mj, rhoa2mj, rhoa3mj, rhoa1mi, rhoa2mi, rhoa3mi;
if (msmeamflag){ if (msmeamflag) {
rhoa1mj = ro0j * t1m_meam[eltj] * MathSpecialKokkos::fm_exp(-beta1m_meam[eltj] * aj) * sij; rhoa1mj = ro0j * t1m_meam[eltj] * MathSpecialKokkos::fm_exp(-beta1m_meam[eltj] * aj) * sij;
rhoa2mj = ro0j * t2m_meam[eltj] * MathSpecialKokkos::fm_exp(-beta2m_meam[eltj] * aj) * sij; rhoa2mj = ro0j * t2m_meam[eltj] * MathSpecialKokkos::fm_exp(-beta2m_meam[eltj] * aj) * sij;
rhoa3mj = ro0j * t3m_meam[eltj] * MathSpecialKokkos::fm_exp(-beta3m_meam[eltj] * aj) * sij; rhoa3mj = ro0j * t3m_meam[eltj] * MathSpecialKokkos::fm_exp(-beta3m_meam[eltj] * aj) * sij;
@ -585,7 +585,7 @@ MEAMKokkos<DeviceType>::calc_rho1(int i, int /*ntype*/, typename AT::t_int_1d ty
const double A2i = rhoa2i / rij2; const double A2i = rhoa2i / rij2;
const double A3i = rhoa3i / (rij2 * rij); const double A3i = rhoa3i / (rij2 * rij);
double A1mj, A2mj, A3mj, A1mi, A2mi, A3mi; double A1mj, A2mj, A3mj, A1mi, A2mi, A3mi;
if (msmeamflag){ if (msmeamflag) {
a_arho2mb[i] += rhoa2mj; a_arho2mb[i] += rhoa2mj;
a_arho2mb[j] += rhoa2mi; a_arho2mb[j] += rhoa2mi;
A1mj = rhoa1mj / rij; A1mj = rhoa1mj / rij;
@ -602,7 +602,7 @@ MEAMKokkos<DeviceType>::calc_rho1(int i, int /*ntype*/, typename AT::t_int_1d ty
a_arho1(j,m) += -A1i * delij[m]; a_arho1(j,m) += -A1i * delij[m];
a_arho3b(i,m) += rhoa3j * delij[m] / rij; a_arho3b(i,m) += rhoa3j * delij[m] / rij;
a_arho3b(j,m) += -rhoa3i * delij[m] / rij; a_arho3b(j,m) += -rhoa3i * delij[m] / rij;
if (msmeamflag){ if (msmeamflag) {
a_arho1m(i,m) += A1mj * delij[m]; a_arho1m(i,m) += A1mj * delij[m];
a_arho1m(j,m) += -A1mi * delij[m]; a_arho1m(j,m) += -A1mi * delij[m];
a_arho3mb(i,m) += rhoa3mj * delij[m] / rij; a_arho3mb(i,m) += rhoa3mj * delij[m] / rij;
@ -611,7 +611,7 @@ MEAMKokkos<DeviceType>::calc_rho1(int i, int /*ntype*/, typename AT::t_int_1d ty
for (int n = m; n < 3; n++) { for (int n = m; n < 3; n++) {
a_arho2(i,nv2) += A2j * delij[m] * delij[n]; a_arho2(i,nv2) += A2j * delij[m] * delij[n];
a_arho2(j,nv2) += A2i * delij[m] * delij[n]; a_arho2(j,nv2) += A2i * delij[m] * delij[n];
if (msmeamflag){ if (msmeamflag) {
a_arho2m(i,nv2) += A2mj * delij[m] * delij[n]; a_arho2m(i,nv2) += A2mj * delij[m] * delij[n];
a_arho2m(j,nv2) += A2mi * delij[m] * delij[n]; a_arho2m(j,nv2) += A2mi * delij[m] * delij[n];
} }
@ -619,7 +619,7 @@ MEAMKokkos<DeviceType>::calc_rho1(int i, int /*ntype*/, typename AT::t_int_1d ty
for (int p = n; p < 3; p++) { for (int p = n; p < 3; p++) {
a_arho3(i,nv3) += A3j * delij[m] * delij[n] * delij[p]; a_arho3(i,nv3) += A3j * delij[m] * delij[n] * delij[p];
a_arho3(j,nv3) += -A3i * delij[m] * delij[n] * delij[p]; a_arho3(j,nv3) += -A3i * delij[m] * delij[n] * delij[p];
if (msmeamflag){ if (msmeamflag) {
a_arho3m(i,nv3) += A3mj * delij[m] * delij[n] * delij[p]; a_arho3m(i,nv3) += A3mj * delij[m] * delij[n] * delij[p];
a_arho3m(j,nv3) += -A3mi * delij[m] * delij[n] * delij[p]; a_arho3m(j,nv3) += -A3mi * delij[m] * delij[n] * delij[p];
} }

29
src/KOKKOS/meam_force_kokkos.h
View File

@ -208,7 +208,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
drhoa2i = -beta2_meam[elti] * invrei * rhoa2i; drhoa2i = -beta2_meam[elti] * invrei * rhoa2i;
rhoa3i = ro0i * MathSpecialKokkos::fm_exp(-beta3_meam[elti] * ai); rhoa3i = ro0i * MathSpecialKokkos::fm_exp(-beta3_meam[elti] * ai);
drhoa3i = -beta3_meam[elti] * invrei * rhoa3i; drhoa3i = -beta3_meam[elti] * invrei * rhoa3i;
if (msmeamflag){ if (msmeamflag) {
rhoa1mi = ro0i * MathSpecialKokkos::fm_exp(-beta1m_meam[elti] * ai) * t1m_meam[elti]; rhoa1mi = ro0i * MathSpecialKokkos::fm_exp(-beta1m_meam[elti] * ai) * t1m_meam[elti];
drhoa1mi = -beta1m_meam[elti] * invrei * rhoa1mi; drhoa1mi = -beta1m_meam[elti] * invrei * rhoa1mi;
rhoa2mi = ro0i * MathSpecialKokkos::fm_exp(-beta2m_meam[elti] * ai) * t2m_meam[elti]; rhoa2mi = ro0i * MathSpecialKokkos::fm_exp(-beta2m_meam[elti] * ai) * t2m_meam[elti];
@ -229,7 +229,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
drhoa2j = -beta2_meam[eltj] * invrej * rhoa2j; drhoa2j = -beta2_meam[eltj] * invrej * rhoa2j;
rhoa3j = ro0j * MathSpecialKokkos::fm_exp(-beta3_meam[eltj] * aj); rhoa3j = ro0j * MathSpecialKokkos::fm_exp(-beta3_meam[eltj] * aj);
drhoa3j = -beta3_meam[eltj] * invrej * rhoa3j; drhoa3j = -beta3_meam[eltj] * invrej * rhoa3j;
if (msmeamflag){ if (msmeamflag) {
rhoa1mj = ro0j * MathSpecialKokkos::fm_exp(-beta1m_meam[eltj] * aj) * t1m_meam[eltj]; rhoa1mj = ro0j * MathSpecialKokkos::fm_exp(-beta1m_meam[eltj] * aj) * t1m_meam[eltj];
drhoa1mj = -beta1m_meam[eltj] * invrej * rhoa1mj; drhoa1mj = -beta1m_meam[eltj] * invrej * rhoa1mj;
rhoa2mj = ro0j * MathSpecialKokkos::fm_exp(-beta2m_meam[eltj] * aj) * t2m_meam[eltj]; rhoa2mj = ro0j * MathSpecialKokkos::fm_exp(-beta2m_meam[eltj] * aj) * t2m_meam[eltj];
@ -246,7 +246,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
drhoa2j = drhoa2i; drhoa2j = drhoa2i;
rhoa3j = rhoa3i; rhoa3j = rhoa3i;
drhoa3j = drhoa3i; drhoa3j = drhoa3i;
if (msmeamflag){ if (msmeamflag) {
rhoa1mj = rhoa1mi; rhoa1mj = rhoa1mi;
drhoa1mj = drhoa1mi; drhoa1mj = drhoa1mi;
rhoa2mj = rhoa2mi; rhoa2mj = rhoa2mi;
@ -400,7 +400,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * d_arho3b(j, m)) * rhoa3i; drho3drm2[m] = (-a3 * drho3drm2[m] + a3a * d_arho3b(j, m)) * rhoa3i;
} }
if (msmeamflag){ if (msmeamflag) {
// rho1m terms // rho1m terms
a1 = 2 * sij / rij; a1 = 2 * sij / rij;
drho1mdr1 = a1 * (drhoa1mj - rhoa1mj / rij) * arg1i1m; drho1mdr1 = a1 * (drhoa1mj - rhoa1mj / rij) * arg1i1m;
@ -455,9 +455,8 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
drho3mdrm1[m] = (a3 * drho3mdrm1[m] - a3a * d_arho3mb(i, m)) * rhoa3mj; drho3mdrm1[m] = (a3 * drho3mdrm1[m] - a3a * d_arho3mb(i, m)) * rhoa3mj;
drho3mdrm2[m] = (-a3 * drho3mdrm2[m] + a3a * d_arho3mb(j, m)) * rhoa3mi; drho3mdrm2[m] = (-a3 * drho3mdrm2[m] + a3a * d_arho3mb(j, m)) * rhoa3mi;
} }
} else{ } else {
for (m = 0; m < 3; m++) {
for (m = 0; m < 3; m++){
drho1mdrm1[m] = 0.0; drho1mdrm1[m] = 0.0;
drho1mdrm2[m] = 0.0; drho1mdrm2[m] = 0.0;
drho2mdrm1[m] = 0.0; drho2mdrm1[m] = 0.0;
@ -470,7 +469,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
// Compute derivatives of weighting functions t wrt rij // Compute derivatives of weighting functions t wrt rij
// Weighting functions t set to unity for msmeam // Weighting functions t set to unity for msmeam
if (msmeamflag){ if (msmeamflag) {
t1i = 1.0; t1i = 1.0;
t2i = 1.0; t2i = 1.0;
t3i = 1.0; t3i = 1.0;
@ -483,7 +482,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
dt2dr2 = 0.0; dt2dr2 = 0.0;
dt3dr1 = 0.0; dt3dr1 = 0.0;
dt3dr2 = 0.0; dt3dr2 = 0.0;
} else{ } else {
t1i = d_t_ave(i, 0); t1i = d_t_ave(i, 0);
t2i = d_t_ave(i, 1); t2i = d_t_ave(i, 1);
@ -538,7 +537,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
get_shpfcn(lattce_meam[eltj][eltj], stheta_meam[elti][elti], ctheta_meam[elti][elti], shpj); get_shpfcn(lattce_meam[eltj][eltj], stheta_meam[elti][elti], ctheta_meam[elti][elti], shpj);
if (msmeamflag){ if (msmeamflag) {
drhodr1 = d_dgamma1[i] * drho0dr1 + drhodr1 = d_dgamma1[i] * drho0dr1 +
d_dgamma2[i] * (dt1dr1 * d_rho1[i] + t1i * (drho1dr1 - drho1mdr1) + d_dgamma2[i] * (dt1dr1 * d_rho1[i] + t1i * (drho1dr1 - drho1mdr1) +
dt2dr1 * d_rho2[i] + t2i * (drho2dr1 - drho2mdr1) + dt2dr1 * d_rho2[i] + t2i * (drho2dr1 - drho2mdr1) +
@ -559,7 +558,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
t2j * (drho2drm2[m] - drho2mdrm2[m]) + t2j * (drho2drm2[m] - drho2mdrm2[m]) +
t3j * (drho3drm2[m] - drho3mdrm2[m]) ); t3j * (drho3drm2[m] - drho3mdrm2[m]) );
} }
} else{ } else {
drhodr1 = d_dgamma1[i] * drho0dr1 + drhodr1 = d_dgamma1[i] * drho0dr1 +
d_dgamma2[i] * d_dgamma2[i] *
(dt1dr1 * d_rho1[i] + t1i * drho1dr1 + dt2dr1 * d_rho2[i] + t2i * drho2dr1 + (dt1dr1 * d_rho1[i] + t1i * drho1dr1 + dt2dr1 * d_rho2[i] + t2i * drho2dr1 +
@ -595,7 +594,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3; drho3ds1 = a3 * rhoa3j * arg1i3 - a3a * rhoa3j * arg3i3;
drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3; drho3ds2 = a3 * rhoa3i * arg1j3 - a3a * rhoa3i * arg3j3;
if (msmeamflag){ if (msmeamflag) {
drho1mds1 = a1 * rhoa1mj * arg1i1m; drho1mds1 = a1 * rhoa1mj * arg1i1m;
drho1mds2 = a1 * rhoa1mi * arg1j1m; drho1mds2 = a1 * rhoa1mi * arg1j1m;
drho2mds1 = a2 * rhoa2mj * arg1i2m - 2.0 / 3.0 * d_arho2mb[i] * rhoa2mj; drho2mds1 = a2 * rhoa2mj * arg1i2m - 2.0 / 3.0 * d_arho2mb[i] * rhoa2mj;
@ -604,7 +603,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
drho3mds2 = a3 * rhoa3mi * arg1j3m - a3a * rhoa3mi * arg3j3m; drho3mds2 = a3 * rhoa3mi * arg1j3m - a3a * rhoa3mi * arg3j3m;
drho3mds1 *= -1; drho3mds1 *= -1;
drho3mds2 *= -1; drho3mds2 *= -1;
} else{ } else {
drho1mds1 = 0.0; drho1mds1 = 0.0;
drho1mds2 = 0.0; drho1mds2 = 0.0;
drho2mds1 = 0.0; drho2mds1 = 0.0;
@ -652,7 +651,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
dt3ds2 = aj * (t3mi - t3j); dt3ds2 = aj * (t3mi - t3j);
} }
if (msmeamflag){ if (msmeamflag) {
drhods1 = d_dgamma1[i] * drho0ds1 + drhods1 = d_dgamma1[i] * drho0ds1 +
d_dgamma2[i] * (dt1ds1 * d_rho1[i] + t1i * (drho1ds1 - drho1mds1) + d_dgamma2[i] * (dt1ds1 * d_rho1[i] + t1i * (drho1ds1 - drho1mds1) +
dt2ds1 * d_rho2[i] + t2i * (drho2ds1 - drho2mds1) + dt2ds1 * d_rho2[i] + t2i * (drho2ds1 - drho2mds1) +
@ -663,7 +662,7 @@ KOKKOS_INLINE_FUNCTION void MEAMKokkos<DeviceType>::operator()(TagMEAMForce<NEIG
dt2ds2 * d_rho2[j] + t2j * (drho2ds2 - drho2mds2) + dt2ds2 * d_rho2[j] + t2j * (drho2ds2 - drho2mds2) +
dt3ds2 * d_rho3[j] + t3j * (drho3ds2 - drho3mds2)) - dt3ds2 * d_rho3[j] + t3j * (drho3ds2 - drho3mds2)) -
d_dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2); d_dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
} else{ } else {
drhods1 = d_dgamma1[i] * drho0ds1 + drhods1 = d_dgamma1[i] * drho0ds1 +
d_dgamma2[i] * d_dgamma2[i] *
(dt1ds1 * d_rho1[i] + t1i * drho1ds1 + dt2ds1 * d_rho2[i] + t2i * drho2ds1 + (dt1ds1 * d_rho1[i] + t1i * drho1ds1 + dt2ds1 * d_rho2[i] + t2i * drho2ds1 +

36
src/KOKKOS/pair_meam_kokkos.cpp
View File

@ -167,7 +167,7 @@ void PairMEAMKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
meam_inst_kk->k_arho3b.template modify<DeviceType>(); meam_inst_kk->k_arho3b.template modify<DeviceType>();
meam_inst_kk->k_t_ave.template modify<DeviceType>(); meam_inst_kk->k_t_ave.template modify<DeviceType>();
meam_inst_kk->k_tsq_ave.template modify<DeviceType>(); meam_inst_kk->k_tsq_ave.template modify<DeviceType>();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.template modify<DeviceType>(); meam_inst_kk->k_arho2mb.template modify<DeviceType>();
meam_inst_kk->k_arho1m.template modify<DeviceType>(); meam_inst_kk->k_arho1m.template modify<DeviceType>();
meam_inst_kk->k_arho2m.template modify<DeviceType>(); meam_inst_kk->k_arho2m.template modify<DeviceType>();
@ -185,7 +185,7 @@ void PairMEAMKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
meam_inst_kk->k_arho3b.template sync<DeviceType>(); meam_inst_kk->k_arho3b.template sync<DeviceType>();
meam_inst_kk->k_t_ave.template sync<DeviceType>(); meam_inst_kk->k_t_ave.template sync<DeviceType>();
meam_inst_kk->k_tsq_ave.template sync<DeviceType>(); meam_inst_kk->k_tsq_ave.template sync<DeviceType>();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.template sync<DeviceType>(); meam_inst_kk->k_arho2mb.template sync<DeviceType>();
meam_inst_kk->k_arho1m.template sync<DeviceType>(); meam_inst_kk->k_arho1m.template sync<DeviceType>();
meam_inst_kk->k_arho2m.template sync<DeviceType>(); meam_inst_kk->k_arho2m.template sync<DeviceType>();
@ -215,7 +215,7 @@ void PairMEAMKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
meam_inst_kk->k_arho3b.template modify<DeviceType>(); meam_inst_kk->k_arho3b.template modify<DeviceType>();
meam_inst_kk->k_t_ave.template modify<DeviceType>(); meam_inst_kk->k_t_ave.template modify<DeviceType>();
meam_inst_kk->k_tsq_ave.template modify<DeviceType>(); meam_inst_kk->k_tsq_ave.template modify<DeviceType>();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.template modify<DeviceType>(); meam_inst_kk->k_arho2mb.template modify<DeviceType>();
meam_inst_kk->k_arho1m.template modify<DeviceType>(); meam_inst_kk->k_arho1m.template modify<DeviceType>();
meam_inst_kk->k_arho2m.template modify<DeviceType>(); meam_inst_kk->k_arho2m.template modify<DeviceType>();
@ -241,7 +241,7 @@ void PairMEAMKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
meam_inst_kk->k_arho3b.template sync<DeviceType>(); meam_inst_kk->k_arho3b.template sync<DeviceType>();
meam_inst_kk->k_t_ave.template sync<DeviceType>(); meam_inst_kk->k_t_ave.template sync<DeviceType>();
meam_inst_kk->k_tsq_ave.template sync<DeviceType>(); meam_inst_kk->k_tsq_ave.template sync<DeviceType>();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.template sync<DeviceType>(); meam_inst_kk->k_arho2mb.template sync<DeviceType>();
meam_inst_kk->k_arho1m.template sync<DeviceType>(); meam_inst_kk->k_arho1m.template sync<DeviceType>();
meam_inst_kk->k_arho2m.template sync<DeviceType>(); meam_inst_kk->k_arho2m.template sync<DeviceType>();
@ -383,7 +383,7 @@ void PairMEAMKokkos<DeviceType>::operator()(TagPairMEAMPackForwardComm, const in
v_buf[m++] = d_tsq_ave(j,0); v_buf[m++] = d_tsq_ave(j,0);
v_buf[m++] = d_tsq_ave(j,1); v_buf[m++] = d_tsq_ave(j,1);
v_buf[m++] = d_tsq_ave(j,2); v_buf[m++] = d_tsq_ave(j,2);
if (msmeamflag){ if (msmeamflag) {
v_buf[m++] = d_arho2mb[j]; v_buf[m++] = d_arho2mb[j];
v_buf[m++] = d_arho1m(j,0); v_buf[m++] = d_arho1m(j,0);
v_buf[m++] = d_arho1m(j,1); v_buf[m++] = d_arho1m(j,1);
@ -448,7 +448,7 @@ void PairMEAMKokkos<DeviceType>::operator()(TagPairMEAMUnpackForwardComm, const
d_tsq_ave(i+first,0) = v_buf[m++]; d_tsq_ave(i+first,0) = v_buf[m++];
d_tsq_ave(i+first,1) = v_buf[m++]; d_tsq_ave(i+first,1) = v_buf[m++];
d_tsq_ave(i+first,2) = v_buf[m++]; d_tsq_ave(i+first,2) = v_buf[m++];
if (msmeamflag){ if (msmeamflag) {
d_arho2mb[i+first] = v_buf[m++]; d_arho2mb[i+first] = v_buf[m++];
d_arho1m(i+first,0) = v_buf[m++]; d_arho1m(i+first,0) = v_buf[m++];
d_arho1m(i+first,1) = v_buf[m++]; d_arho1m(i+first,1) = v_buf[m++];
@ -488,7 +488,7 @@ int PairMEAMKokkos<DeviceType>::pack_forward_comm(int n, int *list, double *buf,
meam_inst_kk->k_arho3b.sync_host(); meam_inst_kk->k_arho3b.sync_host();
meam_inst_kk->k_t_ave.sync_host(); meam_inst_kk->k_t_ave.sync_host();
meam_inst_kk->k_tsq_ave.sync_host(); meam_inst_kk->k_tsq_ave.sync_host();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.sync_host(); meam_inst_kk->k_arho2mb.sync_host();
meam_inst_kk->k_arho1m.sync_host(); meam_inst_kk->k_arho1m.sync_host();
meam_inst_kk->k_arho2m.sync_host(); meam_inst_kk->k_arho2m.sync_host();
@ -528,7 +528,7 @@ int PairMEAMKokkos<DeviceType>::pack_forward_comm(int n, int *list, double *buf,
buf[m++] = meam_inst_kk->h_tsq_ave(j,0); buf[m++] = meam_inst_kk->h_tsq_ave(j,0);
buf[m++] = meam_inst_kk->h_tsq_ave(j,1); buf[m++] = meam_inst_kk->h_tsq_ave(j,1);
buf[m++] = meam_inst_kk->h_tsq_ave(j,2); buf[m++] = meam_inst_kk->h_tsq_ave(j,2);
if (msmeamflag){ if (msmeamflag) {
buf[m++] = meam_inst_kk->h_arho2mb[j]; buf[m++] = meam_inst_kk->h_arho2mb[j];
buf[m++] = meam_inst_kk->h_arho1m(j,0); buf[m++] = meam_inst_kk->h_arho1m(j,0);
buf[m++] = meam_inst_kk->h_arho1m(j,1); buf[m++] = meam_inst_kk->h_arho1m(j,1);
@ -570,7 +570,7 @@ void PairMEAMKokkos<DeviceType>::unpack_forward_comm(int n, int first, double *b
meam_inst_kk->k_arho3b.sync_host(); meam_inst_kk->k_arho3b.sync_host();
meam_inst_kk->k_t_ave.sync_host(); meam_inst_kk->k_t_ave.sync_host();
meam_inst_kk->k_tsq_ave.sync_host(); meam_inst_kk->k_tsq_ave.sync_host();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.sync_host(); meam_inst_kk->k_arho2mb.sync_host();
meam_inst_kk->k_arho1m.sync_host(); meam_inst_kk->k_arho1m.sync_host();
meam_inst_kk->k_arho2m.sync_host(); meam_inst_kk->k_arho2m.sync_host();
@ -610,7 +610,7 @@ void PairMEAMKokkos<DeviceType>::unpack_forward_comm(int n, int first, double *b
meam_inst_kk->h_tsq_ave(i,0) = buf[m++]; meam_inst_kk->h_tsq_ave(i,0) = buf[m++];
meam_inst_kk->h_tsq_ave(i,1) = buf[m++]; meam_inst_kk->h_tsq_ave(i,1) = buf[m++];
meam_inst_kk->h_tsq_ave(i,2) = buf[m++]; meam_inst_kk->h_tsq_ave(i,2) = buf[m++];
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->h_arho2mb[i] = buf[m++]; meam_inst_kk->h_arho2mb[i] = buf[m++];
meam_inst_kk->h_arho1m(i,0) = buf[m++]; meam_inst_kk->h_arho1m(i,0) = buf[m++];
meam_inst_kk->h_arho1m(i,1) = buf[m++]; meam_inst_kk->h_arho1m(i,1) = buf[m++];
@ -644,7 +644,7 @@ void PairMEAMKokkos<DeviceType>::unpack_forward_comm(int n, int first, double *b
meam_inst_kk->k_arho3b.modify_host(); meam_inst_kk->k_arho3b.modify_host();
meam_inst_kk->k_t_ave.modify_host(); meam_inst_kk->k_t_ave.modify_host();
meam_inst_kk->k_tsq_ave.modify_host(); meam_inst_kk->k_tsq_ave.modify_host();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.modify_host(); meam_inst_kk->k_arho2mb.modify_host();
meam_inst_kk->k_arho1m.modify_host(); meam_inst_kk->k_arho1m.modify_host();
meam_inst_kk->k_arho2m.modify_host(); meam_inst_kk->k_arho2m.modify_host();
@ -694,7 +694,7 @@ void PairMEAMKokkos<DeviceType>::operator()(TagPairMEAMPackReverseComm, const in
v_buf[m++] = d_tsq_ave(i+first,0); v_buf[m++] = d_tsq_ave(i+first,0);
v_buf[m++] = d_tsq_ave(i+first,1); v_buf[m++] = d_tsq_ave(i+first,1);
v_buf[m++] = d_tsq_ave(i+first,2); v_buf[m++] = d_tsq_ave(i+first,2);
if (msmeamflag){ if (msmeamflag) {
v_buf[m++] = d_arho2mb[i+first]; v_buf[m++] = d_arho2mb[i+first];
v_buf[m++] = d_arho1m(i+first,0); v_buf[m++] = d_arho1m(i+first,0);
v_buf[m++] = d_arho1m(i+first,1); v_buf[m++] = d_arho1m(i+first,1);
@ -725,7 +725,7 @@ int PairMEAMKokkos<DeviceType>::pack_reverse_comm(int n, int first, double *buf)
meam_inst_kk->k_arho3b.sync_host(); meam_inst_kk->k_arho3b.sync_host();
meam_inst_kk->k_t_ave.sync_host(); meam_inst_kk->k_t_ave.sync_host();
meam_inst_kk->k_tsq_ave.sync_host(); meam_inst_kk->k_tsq_ave.sync_host();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.sync_host(); meam_inst_kk->k_arho2mb.sync_host();
meam_inst_kk->k_arho1m.sync_host(); meam_inst_kk->k_arho1m.sync_host();
meam_inst_kk->k_arho2m.sync_host(); meam_inst_kk->k_arho2m.sync_host();
@ -757,7 +757,7 @@ int PairMEAMKokkos<DeviceType>::pack_reverse_comm(int n, int first, double *buf)
buf[m++] = meam_inst_kk->h_tsq_ave(i,0); buf[m++] = meam_inst_kk->h_tsq_ave(i,0);
buf[m++] = meam_inst_kk->h_tsq_ave(i,1); buf[m++] = meam_inst_kk->h_tsq_ave(i,1);
buf[m++] = meam_inst_kk->h_tsq_ave(i,2); buf[m++] = meam_inst_kk->h_tsq_ave(i,2);
if (msmeamflag){ if (msmeamflag) {
buf[m++] = meam_inst_kk->h_arho2mb[i]; buf[m++] = meam_inst_kk->h_arho2mb[i];
buf[m++] = meam_inst_kk->h_arho1m(i,0); buf[m++] = meam_inst_kk->h_arho1m(i,0);
buf[m++] = meam_inst_kk->h_arho1m(i,1); buf[m++] = meam_inst_kk->h_arho1m(i,1);
@ -819,7 +819,7 @@ void PairMEAMKokkos<DeviceType>::operator()(TagPairMEAMUnpackReverseComm, const
d_tsq_ave(j,0) += v_buf[m++]; d_tsq_ave(j,0) += v_buf[m++];
d_tsq_ave(j,1) += v_buf[m++]; d_tsq_ave(j,1) += v_buf[m++];
d_tsq_ave(j,2) += v_buf[m++]; d_tsq_ave(j,2) += v_buf[m++];
if (msmeamflag){ if (msmeamflag) {
d_arho2mb[j] += v_buf[m++]; d_arho2mb[j] += v_buf[m++];
d_arho1m(j,0) += v_buf[m++]; d_arho1m(j,0) += v_buf[m++];
d_arho1m(j,1) += v_buf[m++]; d_arho1m(j,1) += v_buf[m++];
@ -850,7 +850,7 @@ void PairMEAMKokkos<DeviceType>::unpack_reverse_comm(int n, int *list, double *b
meam_inst_kk->k_arho3b.sync_host(); meam_inst_kk->k_arho3b.sync_host();
meam_inst_kk->k_t_ave.sync_host(); meam_inst_kk->k_t_ave.sync_host();
meam_inst_kk->k_tsq_ave.sync_host(); meam_inst_kk->k_tsq_ave.sync_host();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.sync_host(); meam_inst_kk->k_arho2mb.sync_host();
meam_inst_kk->k_arho1m.sync_host(); meam_inst_kk->k_arho1m.sync_host();
meam_inst_kk->k_arho2m.sync_host(); meam_inst_kk->k_arho2m.sync_host();
@ -882,7 +882,7 @@ void PairMEAMKokkos<DeviceType>::unpack_reverse_comm(int n, int *list, double *b
meam_inst_kk->h_tsq_ave(j,0) += buf[m++]; meam_inst_kk->h_tsq_ave(j,0) += buf[m++];
meam_inst_kk->h_tsq_ave(j,1) += buf[m++]; meam_inst_kk->h_tsq_ave(j,1) += buf[m++];
meam_inst_kk->h_tsq_ave(j,2) += buf[m++]; meam_inst_kk->h_tsq_ave(j,2) += buf[m++];
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->h_arho2mb[j] += buf[m++]; meam_inst_kk->h_arho2mb[j] += buf[m++];
meam_inst_kk->h_arho1m(j,0) += buf[m++]; meam_inst_kk->h_arho1m(j,0) += buf[m++];
meam_inst_kk->h_arho1m(j,1) += buf[m++]; meam_inst_kk->h_arho1m(j,1) += buf[m++];
@ -908,7 +908,7 @@ void PairMEAMKokkos<DeviceType>::unpack_reverse_comm(int n, int *list, double *b
meam_inst_kk->k_arho3b.modify_host(); meam_inst_kk->k_arho3b.modify_host();
meam_inst_kk->k_t_ave.modify_host(); meam_inst_kk->k_t_ave.modify_host();
meam_inst_kk->k_tsq_ave.modify_host(); meam_inst_kk->k_tsq_ave.modify_host();
if (msmeamflag){ if (msmeamflag) {
meam_inst_kk->k_arho2mb.modify_host(); meam_inst_kk->k_arho2mb.modify_host();
meam_inst_kk->k_arho1m.modify_host(); meam_inst_kk->k_arho1m.modify_host();
meam_inst_kk->k_arho2m.modify_host(); meam_inst_kk->k_arho2m.modify_host();

4
src/MEAM/meam_dens_final.cpp
View File

@ -27,7 +27,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
// Complete the calculation of density // Complete the calculation of density
if (this->msmeamflag){ if (this->msmeamflag) {
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
elti = fmap[type[i]]; elti = fmap[type[i]];
if (elti >= 0) { if (elti >= 0) {
@ -131,7 +131,7 @@ MEAM::meam_dens_final(int nlocal, int eflag_either, int eflag_global, int eflag_
} }
} }
} }
} else{ } else {
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
elti = fmap[type[i]]; elti = fmap[type[i]];
if (elti >= 0) { if (elti >= 0) {

24
src/MEAM/meam_dens_init.cpp
View File

@ -46,7 +46,7 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
memory->destroy(t_ave); memory->destroy(t_ave);
memory->destroy(tsq_ave); memory->destroy(tsq_ave);
// msmeam params // msmeam params
if (this->msmeamflag){ if (this->msmeamflag) {
memory->destroy(arho1m); memory->destroy(arho1m);
memory->destroy(arho2m); memory->destroy(arho2m);
memory->destroy(arho3m); memory->destroy(arho3m);
@ -74,7 +74,7 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
memory->create(t_ave, nmax, 3, "pair:t_ave"); memory->create(t_ave, nmax, 3, "pair:t_ave");
memory->create(tsq_ave, nmax, 3, "pair:tsq_ave"); memory->create(tsq_ave, nmax, 3, "pair:tsq_ave");
// msmeam params // msmeam params
if (this->msmeamflag){ if (this->msmeamflag) {
memory->create(arho1m, nmax, 3, "pair:arho1m"); memory->create(arho1m, nmax, 3, "pair:arho1m");
memory->create(arho2m, nmax, 6, "pair:arho2m"); memory->create(arho2m, nmax, 6, "pair:arho2m");
memory->create(arho3m, nmax, 10, "pair:arho3m"); memory->create(arho3m, nmax, 10, "pair:arho3m");
@ -99,24 +99,24 @@ MEAM::meam_dens_setup(int atom_nmax, int nall, int n_neigh)
rho0[i] = 0.0; rho0[i] = 0.0;
arho2b[i] = 0.0; arho2b[i] = 0.0;
arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0; arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
if (this->msmeamflag){ if (this->msmeamflag) {
arho2mb[i] = 0.0; arho2mb[i] = 0.0;
arho1m[i][0] = arho1m[i][1] = arho1m[i][2] = 0.0; arho1m[i][0] = arho1m[i][1] = arho1m[i][2] = 0.0;
} }
for (j = 0; j < 6; j++){ for (j = 0; j < 6; j++) {
arho2[i][j] = 0.0; arho2[i][j] = 0.0;
if (this->msmeamflag){ if (this->msmeamflag) {
arho2m[i][j] = 0.0; arho2m[i][j] = 0.0;
} }
} }
for (j = 0; j < 10; j++){ for (j = 0; j < 10; j++) {
arho3[i][j] = 0.0; arho3[i][j] = 0.0;
if (this->msmeamflag){ if (this->msmeamflag) {
arho3m[i][j] = 0.0; arho3m[i][j] = 0.0;
} }
} }
arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0; arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
if (this->msmeamflag){ if (this->msmeamflag) {
arho3mb[i][0] = arho3mb[i][1] = arho3mb[i][2] = 0.0; arho3mb[i][0] = arho3mb[i][1] = arho3mb[i][2] = 0.0;
} }
t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0; t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
@ -394,7 +394,7 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn
A3i = rhoa3i / (rij2 * rij); A3i = rhoa3i / (rij2 * rij);
nv2 = 0; nv2 = 0;
nv3 = 0; nv3 = 0;
if (this->msmeamflag){ if (this->msmeamflag) {
arho2mb[i] = arho2mb[i] + rhoa2mj; arho2mb[i] = arho2mb[i] + rhoa2mj;
arho2mb[j] = arho2mb[j] + rhoa2mi; arho2mb[j] = arho2mb[j] + rhoa2mi;
A1mj = rhoa1mj/rij; A1mj = rhoa1mj/rij;
@ -409,7 +409,7 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn
arho1[j][m] = arho1[j][m] - A1i * delij[m]; arho1[j][m] = arho1[j][m] - A1i * delij[m];
arho3b[i][m] = arho3b[i][m] + rhoa3j * delij[m] / rij; arho3b[i][m] = arho3b[i][m] + rhoa3j * delij[m] / rij;
arho3b[j][m] = arho3b[j][m] - rhoa3i * delij[m] / rij; arho3b[j][m] = arho3b[j][m] - rhoa3i * delij[m] / rij;
if (this->msmeamflag){ if (this->msmeamflag) {
arho1m[i][m] = arho1m[i][m] + A1mj*delij[m]; arho1m[i][m] = arho1m[i][m] + A1mj*delij[m];
arho1m[j][m] = arho1m[j][m] - A1mi*delij[m]; arho1m[j][m] = arho1m[j][m] - A1mi*delij[m];
arho3mb[i][m] = arho3mb[i][m] + rhoa3mj*delij[m] / rij; arho3mb[i][m] = arho3mb[i][m] + rhoa3mj*delij[m] / rij;
@ -418,7 +418,7 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn
for (n = m; n < 3; n++) { for (n = m; n < 3; n++) {
arho2[i][nv2] = arho2[i][nv2] + A2j * delij[m] * delij[n]; arho2[i][nv2] = arho2[i][nv2] + A2j * delij[m] * delij[n];
arho2[j][nv2] = arho2[j][nv2] + A2i * delij[m] * delij[n]; arho2[j][nv2] = arho2[j][nv2] + A2i * delij[m] * delij[n];
if (this->msmeamflag){ if (this->msmeamflag) {
arho2m[i][nv2] = arho2m[i][nv2] + A2mj*delij[m] * delij[n]; arho2m[i][nv2] = arho2m[i][nv2] + A2mj*delij[m] * delij[n];
arho2m[j][nv2] = arho2m[j][nv2] + A2mi*delij[m] * delij[n]; arho2m[j][nv2] = arho2m[j][nv2] + A2mi*delij[m] * delij[n];
} }
@ -426,7 +426,7 @@ MEAM::calc_rho1(int i, int /*ntype*/, int* type, int* fmap, double** x, int numn
for (p = n; p < 3; p++) { for (p = n; p < 3; p++) {
arho3[i][nv3] = arho3[i][nv3] + A3j * delij[m] * delij[n] * delij[p]; arho3[i][nv3] = arho3[i][nv3] + A3j * delij[m] * delij[n] * delij[p];
arho3[j][nv3] = arho3[j][nv3] - A3i * delij[m] * delij[n] * delij[p]; arho3[j][nv3] = arho3[j][nv3] - A3i * delij[m] * delij[n] * delij[p];
if (this->msmeamflag){ if (this->msmeamflag) {
arho3m[i][nv3] = arho3m[i][nv3] + A3mj*delij[m]*delij[n]*delij[p]; arho3m[i][nv3] = arho3m[i][nv3] + A3mj*delij[m]*delij[n]*delij[p];
arho3m[j][nv3] = arho3m[j][nv3] - A3mi*delij[m]*delij[n]*delij[p]; arho3m[j][nv3] = arho3m[j][nv3] - A3mi*delij[m]*delij[n]*delij[p];
} }

10
src/MEAM/meam_force.cpp
View File

@ -138,7 +138,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai); rhoa3i = ro0i * MathSpecial::fm_exp(-this->beta3_meam[elti] * ai);
drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i; drhoa3i = -this->beta3_meam[elti] * invrei * rhoa3i;
if (this->msmeamflag){ if (this->msmeamflag) {
rhoa1mi = ro0i * MathSpecial::fm_exp(-this->beta1m_meam[elti] * ai) * t1m_meam[elti]; rhoa1mi = ro0i * MathSpecial::fm_exp(-this->beta1m_meam[elti] * ai) * t1m_meam[elti];
drhoa1mi = -this->beta1m_meam[elti] * invrei * rhoa1mi; drhoa1mi = -this->beta1m_meam[elti] * invrei * rhoa1mi;
rhoa2mi = ro0i * MathSpecial::fm_exp(-this->beta2m_meam[elti] * ai) * t2m_meam[elti]; rhoa2mi = ro0i * MathSpecial::fm_exp(-this->beta2m_meam[elti] * ai) * t2m_meam[elti];
@ -471,7 +471,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpi); get_shpfcn(this->lattce_meam[elti][elti], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpi);
get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpj); get_shpfcn(this->lattce_meam[eltj][eltj], this->stheta_meam[elti][elti], this->ctheta_meam[elti][elti], shpj);
if (this->msmeamflag){ if (this->msmeamflag) {
drhodr1 = dgamma1[i] * drho0dr1 + drhodr1 = dgamma1[i] * drho0dr1 +
dgamma2[i] * (dt1dr1 * rho1[i] + t1i * (drho1dr1 - drho1mdr1) + dgamma2[i] * (dt1dr1 * rho1[i] + t1i * (drho1dr1 - drho1mdr1) +
dt2dr1 * rho2[i] + t2i * (drho2dr1 - drho2mdr1) + dt2dr1 * rho2[i] + t2i * (drho2dr1 - drho2mdr1) +
@ -492,7 +492,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
t2j * (drho2drm2[m] - drho2mdrm2[m]) + t2j * (drho2drm2[m] - drho2mdrm2[m]) +
t3j * (drho3drm2[m] - drho3mdrm2[m]) ); t3j * (drho3drm2[m] - drho3mdrm2[m]) );
} }
} else{ } else {
drhodr1 = dgamma1[i] * drho0dr1 + drhodr1 = dgamma1[i] * drho0dr1 +
dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 + dgamma2[i] * (dt1dr1 * rho1[i] + t1i * drho1dr1 + dt2dr1 * rho2[i] + t2i * drho2dr1 +
@ -586,7 +586,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
dt3ds2 = aj * (t3mi - t3j); dt3ds2 = aj * (t3mi - t3j);
} }
if (this->msmeamflag){ if (this->msmeamflag) {
drhods1 = dgamma1[i] * drho0ds1 + drhods1 = dgamma1[i] * drho0ds1 +
dgamma2[i] * (dt1ds1 * rho1[i] + t1i * (drho1ds1 - drho1mds1) + dgamma2[i] * (dt1ds1 * rho1[i] + t1i * (drho1ds1 - drho1mds1) +
dt2ds1 * rho2[i] + t2i * (drho2ds1 - drho2mds1) + dt2ds1 * rho2[i] + t2i * (drho2ds1 - drho2mds1) +
@ -598,7 +598,7 @@ MEAM::meam_force(int i, int eflag_global, int eflag_atom, int vflag_global, int
dt3ds2 * rho3[j] + t3j * (drho3ds2 - drho3mds2)) - dt3ds2 * rho3[j] + t3j * (drho3ds2 - drho3mds2)) -
dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2); dgamma3[j] * (shpj[0] * dt1ds2 + shpj[1] * dt2ds2 + shpj[2] * dt3ds2);
} }
else{ else {
drhods1 = dgamma1[i] * drho0ds1 + drhods1 = dgamma1[i] * drho0ds1 +
dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + dt2ds1 * rho2[i] + t2i * drho2ds1 + dgamma2[i] * (dt1ds1 * rho1[i] + t1i * drho1ds1 + dt2ds1 * rho2[i] + t2i * drho2ds1 +
dt3ds1 * rho3[i] + t3i * drho3ds1) - dt3ds1 * rho3[i] + t3i * drho3ds1) -

2
src/MEAM/meam_impl.cpp
View File

@ -96,7 +96,7 @@ MEAM::~MEAM()
memory->destroy(this->fcpair); memory->destroy(this->fcpair);
// msmeam // msmeam
if (this->msmeamflag){ if (this->msmeamflag){
memory->destroy(this->arho1m); memory->destroy(this->arho1m);
memory->destroy(this->arho2m); memory->destroy(this->arho2m);
memory->destroy(this->arho3m); memory->destroy(this->arho3m);

27
src/MEAM/meam_setup_done.cpp
View File

@ -348,11 +348,11 @@ double MEAM::phi_meam(double r, int a, int b)
Z12 = get_Zij(this->lattce_meam[a][b]); Z12 = get_Zij(this->lattce_meam[a][b]);
// this function has extra args for msmeam // this function has extra args for msmeam
if (this->msmeamflag){ if (this->msmeamflag) {
get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32, get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32,
&rho1m1, &rho2m1, &rho3m1, &rho1m1, &rho2m1, &rho3m1,
&rho1m2, &rho2m2, &rho3m2); &rho1m2, &rho2m2, &rho3m2);
} else{ } else {
get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32, get_densref(r, a, b, &rho01, &rho11, &rho21, &rho31, &rho02, &rho12, &rho22, &rho32,
nullptr, nullptr, nullptr, nullptr, nullptr, nullptr,
nullptr, nullptr, nullptr); nullptr, nullptr, nullptr);
@ -385,7 +385,7 @@ double MEAM::phi_meam(double r, int a, int b)
this->t3_meam[a], this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], r, a, b, this->t3_meam[a], this->t1_meam[b], this->t2_meam[b], this->t3_meam[b], r, a, b,
this->lattce_meam[a][b]); this->lattce_meam[a][b]);
// with msmeam call twice with different sets of variables // with msmeam call twice with different sets of variables
if (this->msmeamflag){ if (this->msmeamflag) {
get_tavref(&t1m1av, &t2m1av, &t3m1av, &t1m2av, &t2m2av, &t3m2av, this->t1m_meam[a], this->t2m_meam[a], get_tavref(&t1m1av, &t2m1av, &t3m1av, &t1m2av, &t2m2av, &t3m2av, this->t1m_meam[a], this->t2m_meam[a],
this->t3m_meam[a], this->t1m_meam[b], this->t2m_meam[b], this->t3m_meam[b], r, a, b, this->t3m_meam[a], this->t1m_meam[b], this->t2m_meam[b], this->t3m_meam[b], r, a, b,
this->lattce_meam[a][b]); this->lattce_meam[a][b]);
@ -437,7 +437,7 @@ double MEAM::phi_meam(double r, int a, int b)
rho0_2 = this->rho0_meam[b] * Z2 * G2; rho0_2 = this->rho0_meam[b] * Z2 * G2;
} }
if (this->msmeamflag){ if (this->msmeamflag) {
// no additional use of t's here; all included in definitions of rho's for msmeam // no additional use of t's here; all included in definitions of rho's for msmeam
Gam1 = rho11 + rho21 + rho31 - (rho1m1 + rho2m1 + rho3m1); Gam1 = rho11 + rho21 + rho31 - (rho1m1 + rho2m1 + rho3m1);
if (rho01 < 1.0e-14) if (rho01 < 1.0e-14)
@ -450,7 +450,7 @@ double MEAM::phi_meam(double r, int a, int b)
else else
Gam2 = Gam2 / (rho02 * rho02); Gam2 = Gam2 / (rho02 * rho02);
} else{ } else {
Gam1 = (t11av * rho11 + t21av * rho21 + t31av * rho31); Gam1 = (t11av * rho11 + t21av * rho21 + t31av * rho31);
if (rho01 < 1.0e-14) if (rho01 < 1.0e-14)
Gam1 = 0.0; Gam1 = 0.0;
@ -708,7 +708,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1); rhoa01 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta0_meam[a] * a1);
if (this->msmeamflag){ if (this->msmeamflag) {
// the rho variables are multiplied by t here since ialloy not needed in msmeam // the rho variables are multiplied by t here since ialloy not needed in msmeam
rhoa11 = this->rho0_meam[a] * this->t1_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1); rhoa11 = this->rho0_meam[a] * this->t1_meam[a] * MathSpecial::fm_exp(-this->beta1_meam[a] * a1);
rhoa21 = this->rho0_meam[a] * this->t2_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1); rhoa21 = this->rho0_meam[a] * this->t2_meam[a] * MathSpecial::fm_exp(-this->beta2_meam[a] * a1);
@ -724,8 +724,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
rhoa1m2 = this->rho0_meam[b] * this->t1m_meam[b] * MathSpecial::fm_exp(-this->beta1m_meam[b] * a2); rhoa1m2 = this->rho0_meam[b] * this->t1m_meam[b] * MathSpecial::fm_exp(-this->beta1m_meam[b] * a2);
rhoa2m2 = this->rho0_meam[b] * this->t2m_meam[b] * MathSpecial::fm_exp(-this->beta2m_meam[b] * a2); rhoa2m2 = this->rho0_meam[b] * this->t2m_meam[b] * MathSpecial::fm_exp(-this->beta2m_meam[b] * a2);
rhoa3m2 = this->rho0_meam[b] * this->t3m_meam[b] * MathSpecial::fm_exp(-this->beta3m_meam[b] * a2); rhoa3m2 = this->rho0_meam[b] * this->t3m_meam[b] * MathSpecial::fm_exp(-this->beta3m_meam[b] * a2);
} } else {
else{
rhoa11 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta1_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); 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); rhoa31 = this->rho0_meam[a] * MathSpecial::fm_exp(-this->beta3_meam[a] * a1);
@ -745,7 +744,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho12 = 0.0; *rho12 = 0.0;
*rho22 = 0.0; *rho22 = 0.0;
*rho32 = 0.0; *rho32 = 0.0;
if (this->msmeamflag){ if (this->msmeamflag) {
*rho1m1 = 0.0; *rho1m1 = 0.0;
*rho2m1 = 0.0; *rho2m1 = 0.0;
*rho3m1 = 0.0; *rho3m1 = 0.0;
@ -775,7 +774,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho02 = 4.0 * rhoa01; *rho02 = 4.0 * rhoa01;
*rho31 = 32.0 / 9.0 * rhoa32 * rhoa32; *rho31 = 32.0 / 9.0 * rhoa32 * rhoa32;
*rho32 = 32.0 / 9.0 * rhoa31 * rhoa31; *rho32 = 32.0 / 9.0 * rhoa31 * rhoa31;
if (this->msmeamflag){ if (this->msmeamflag) {
*rho3m1 = 32.0 / 9.0 * rhoa3m2 * rhoa3m2; *rho3m1 = 32.0 / 9.0 * rhoa3m2 * rhoa3m2;
*rho3m2 = 32.0 / 9.0 * rhoa3m1 * rhoa3m1; *rho3m2 = 32.0 / 9.0 * rhoa3m1 * rhoa3m1;
} }
@ -785,7 +784,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho02 = 12 * rhoa01; *rho02 = 12 * rhoa01;
*rho31 = 1.0 / 3.0 * rhoa32 * rhoa32; *rho31 = 1.0 / 3.0 * rhoa32 * rhoa32;
*rho32 = 1.0 / 3.0 * rhoa31 * rhoa31; *rho32 = 1.0 / 3.0 * rhoa31 * rhoa31;
if (this->msmeamflag){ if (this->msmeamflag) {
*rho3m1 = 1.0 / 3.0 * rhoa3m2 * rhoa3m2; *rho3m1 = 1.0 / 3.0 * rhoa3m2 * rhoa3m2;
*rho3m2 = 1.0 / 3.0 * rhoa3m1 * rhoa3m1; *rho3m2 = 1.0 / 3.0 * rhoa3m1 * rhoa3m1;
} }
@ -800,7 +799,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho22 = s[1] * rhoa21 * rhoa21; *rho22 = s[1] * rhoa21 * rhoa21;
*rho31 = s[2] * rhoa32 * rhoa32; *rho31 = s[2] * rhoa32 * rhoa32;
*rho32 = s[2] * rhoa31 * rhoa31; *rho32 = s[2] * rhoa31 * rhoa31;
if (this->msmeamflag){ if (this->msmeamflag) {
*rho1m1 = s[0] * rhoa1m2 * rhoa1m2; *rho1m1 = s[0] * rhoa1m2 * rhoa1m2;
*rho1m2 = s[0] * rhoa1m1 * rhoa1m1; *rho1m2 = s[0] * rhoa1m1 * rhoa1m1;
*rho2m1 = s[1] * rhoa2m2 * rhoa2m2; *rho2m1 = s[1] * rhoa2m2 * rhoa2m2;
@ -818,7 +817,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho22 = rhoa22; *rho22 = rhoa22;
*rho31 = rhoa31; *rho31 = rhoa31;
*rho32 = rhoa32; *rho32 = rhoa32;
if (this->msmeamflag){ if (this->msmeamflag) {
*rho1m1 = rhoa1m1; *rho1m1 = rhoa1m1;
*rho1m2 = rhoa1m2; *rho1m2 = rhoa1m2;
*rho2m1 = rhoa2m1; *rho2m1 = rhoa2m1;
@ -837,7 +836,7 @@ void MEAM::get_densref(double r, int a, int b, double* rho01, double* rho11, dou
*rho21 = *rho21 / denom * *rho01; *rho21 = *rho21 / denom * *rho01;
} else } else
*rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22); *rho21 = 8. / 3. * (rhoa21 - rhoa22) * (rhoa21 - rhoa22);
if (this->msmeamflag){ if (this->msmeamflag) {
*rho2m1 = 8. / 3. * (rhoa2m1 - rhoa2m2) * (rhoa2m1 - rhoa2m2); *rho2m1 = 8. / 3. * (rhoa2m1 - rhoa2m2) * (rhoa2m1 - rhoa2m2);
} }
break; break;

20
src/MEAM/pair_meam.cpp
View File

@ -196,10 +196,10 @@ void PairMEAM::settings(int narg, char **arg)
// set comm size needed by this Pair // set comm size needed by this Pair
if (msmeamflag){ if (msmeamflag) {
comm_forward = 38+23; // plus 23 for msmeam comm_forward = 38+23; // plus 23 for msmeam
comm_reverse = 30+23; // plus 23 for msmeam comm_reverse = 30+23; // plus 23 for msmeam
} else{ } else {
comm_forward = 38; comm_forward = 38;
comm_reverse = 30; comm_reverse = 30;
} }
@ -449,7 +449,7 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
b1[index] = values.next_double(); b1[index] = values.next_double();
b2[index] = values.next_double(); b2[index] = values.next_double();
b3[index] = values.next_double(); b3[index] = values.next_double();
if (msmeamflag){ if (msmeamflag) {
b1m[index] = values.next_double(); b1m[index] = values.next_double();
b2m[index] = values.next_double(); b2m[index] = values.next_double();
b3m[index] = values.next_double(); b3m[index] = values.next_double();
@ -461,7 +461,7 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
t1[index] = values.next_double(); t1[index] = values.next_double();
t2[index] = values.next_double(); t2[index] = values.next_double();
t3[index] = values.next_double(); t3[index] = values.next_double();
if (msmeamflag){ if (msmeamflag) {
t1m[index] = values.next_double(); t1m[index] = values.next_double();
t2m[index] = values.next_double(); t2m[index] = values.next_double();
t3m[index] = values.next_double(); t3m[index] = values.next_double();
@ -524,13 +524,13 @@ void PairMEAM::read_global_meam_file(const std::string &globalfile)
// pass element parameters to MEAM package // pass element parameters to MEAM package
if (msmeamflag){ if (msmeamflag) {
meam_inst->meam_setup_global(nlibelements, lat.data(), ielement.data(), atwt.data(), meam_inst->meam_setup_global(nlibelements, lat.data(), ielement.data(), atwt.data(),
alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(), alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
alat.data(), esub.data(), asub.data(), t0.data(), t1.data(), alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
t2.data(), t3.data(), rozero.data(), ibar.data(), b1m.data(), t2.data(), t3.data(), rozero.data(), ibar.data(), b1m.data(),
b2m.data(), b3m.data(), t1m.data(), t2m.data(), t3m.data()); b2m.data(), b3m.data(), t1m.data(), t2m.data(), t3m.data());
} else{ } else {
meam_inst->meam_setup_global(nlibelements, lat.data(), ielement.data(), atwt.data(), meam_inst->meam_setup_global(nlibelements, lat.data(), ielement.data(), atwt.data(),
alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(), alpha.data(), b0.data(), b1.data(), b2.data(), b3.data(),
alat.data(), esub.data(), asub.data(), t0.data(), t1.data(), alat.data(), esub.data(), asub.data(), t0.data(), t1.data(),
@ -659,7 +659,7 @@ int PairMEAM::pack_forward_comm(int n, int *list, double *buf,
buf[m++] = meam_inst->tsq_ave[j][0]; buf[m++] = meam_inst->tsq_ave[j][0];
buf[m++] = meam_inst->tsq_ave[j][1]; buf[m++] = meam_inst->tsq_ave[j][1];
buf[m++] = meam_inst->tsq_ave[j][2]; buf[m++] = meam_inst->tsq_ave[j][2];
if (msmeamflag){ if (msmeamflag) {
buf[m++] = meam_inst->arho2mb[j]; buf[m++] = meam_inst->arho2mb[j];
buf[m++] = meam_inst->arho1m[j][0]; buf[m++] = meam_inst->arho1m[j][0];
buf[m++] = meam_inst->arho1m[j][1]; buf[m++] = meam_inst->arho1m[j][1];
@ -719,7 +719,7 @@ void PairMEAM::unpack_forward_comm(int n, int first, double *buf)
meam_inst->tsq_ave[i][0] = buf[m++]; meam_inst->tsq_ave[i][0] = buf[m++];
meam_inst->tsq_ave[i][1] = buf[m++]; meam_inst->tsq_ave[i][1] = buf[m++];
meam_inst->tsq_ave[i][2] = buf[m++]; meam_inst->tsq_ave[i][2] = buf[m++];
if (msmeamflag){ if (msmeamflag) {
meam_inst->arho2mb[i] = buf[m++]; meam_inst->arho2mb[i] = buf[m++];
meam_inst->arho1m[i][0] = buf[m++]; meam_inst->arho1m[i][0] = buf[m++];
meam_inst->arho1m[i][1] = buf[m++]; meam_inst->arho1m[i][1] = buf[m++];
@ -768,7 +768,7 @@ int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
buf[m++] = meam_inst->tsq_ave[i][0]; buf[m++] = meam_inst->tsq_ave[i][0];
buf[m++] = meam_inst->tsq_ave[i][1]; buf[m++] = meam_inst->tsq_ave[i][1];
buf[m++] = meam_inst->tsq_ave[i][2]; buf[m++] = meam_inst->tsq_ave[i][2];
if (msmeamflag){ if (msmeamflag) {
buf[m++] = meam_inst->arho2mb[i]; buf[m++] = meam_inst->arho2mb[i];
buf[m++] = meam_inst->arho1m[i][0]; buf[m++] = meam_inst->arho1m[i][0];
buf[m++] = meam_inst->arho1m[i][1]; buf[m++] = meam_inst->arho1m[i][1];
@ -819,7 +819,7 @@ void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
meam_inst->tsq_ave[j][0] += buf[m++]; meam_inst->tsq_ave[j][0] += buf[m++];
meam_inst->tsq_ave[j][1] += buf[m++]; meam_inst->tsq_ave[j][1] += buf[m++];
meam_inst->tsq_ave[j][2] += buf[m++]; meam_inst->tsq_ave[j][2] += buf[m++];
if (msmeamflag){ if (msmeamflag) {
meam_inst->arho2mb[j] += buf[m++]; meam_inst->arho2mb[j] += buf[m++];
meam_inst->arho1m[j][0] += buf[m++]; meam_inst->arho1m[j][0] += buf[m++];
meam_inst->arho1m[j][1] += buf[m++]; meam_inst->arho1m[j][1] += buf[m++];