Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80783434
pair_cdeam.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
Mon, Sep 2, 14:06
Size
26 KB
Mime Type
text/x-c
Expires
Wed, Sep 4, 14:06 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
18531325
Attached To
rLAMMPS lammps
pair_cdeam.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Alexander Stukowski
Technical University of Darmstadt,
Germany Department of Materials Science
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_cdeam.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
// This is for debugging purposes. The ASSERT() macro is used in the code to check
// if everything runs as expected. Change this to #if 0 if you don't need the checking.
#if 0
#define ASSERT(cond) ((!(cond)) ? my_failure(error,__FILE__,__LINE__) : my_noop())
inline void my_noop() {}
inline void my_failure(Error* error, const char* file, int line) {
char str[1024];
sprintf(str,"Assertion failure: File %s, line %i", file, line);
error->one(FLERR,str);
}
#else
#define ASSERT(cond)
#endif
#define MAXLINE 1024
// This sets the maximum line length in EAM input files.
PairCDEAM
::
PairCDEAM
(
LAMMPS
*
lmp
,
int
_cdeamVersion
)
:
PairEAM
(
lmp
),
PairEAMAlloy
(
lmp
),
cdeamVersion
(
_cdeamVersion
)
{
single_enable
=
0
;
restartinfo
=
0
;
rhoB
=
NULL
;
D_values
=
NULL
;
hcoeff
=
NULL
;
// Set communication buffer sizes needed by this pair style.
if
(
cdeamVersion
==
1
)
{
comm_forward
=
4
;
comm_reverse
=
3
;
}
else
if
(
cdeamVersion
==
2
)
{
comm_forward
=
3
;
comm_reverse
=
2
;
}
else
{
error
->
all
(
FLERR
,
"Invalid CD-EAM potential version."
);
}
}
PairCDEAM
::~
PairCDEAM
()
{
memory
->
destroy
(
rhoB
);
memory
->
destroy
(
D_values
);
if
(
hcoeff
)
delete
[]
hcoeff
;
}
void
PairCDEAM
::
compute
(
int
eflag
,
int
vflag
)
{
int
i
,
j
,
ii
,
jj
,
inum
,
jnum
,
itype
,
jtype
;
double
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
,
evdwl
,
fpair
;
double
rsq
,
rhoip
,
rhojp
,
recip
,
phi
;
int
*
ilist
,
*
jlist
,
*
numneigh
,
**
firstneigh
;
evdwl
=
0.0
;
if
(
eflag
||
vflag
)
ev_setup
(
eflag
,
vflag
);
else
evflag
=
vflag_fdotr
=
eflag_global
=
eflag_atom
=
0
;
// Grow per-atom arrays if necessary
if
(
atom
->
nmax
>
nmax
)
{
memory
->
destroy
(
rho
);
memory
->
destroy
(
fp
);
memory
->
destroy
(
rhoB
);
memory
->
destroy
(
D_values
);
nmax
=
atom
->
nmax
;
memory
->
create
(
rho
,
nmax
,
"pair:rho"
);
memory
->
create
(
rhoB
,
nmax
,
"pair:rhoB"
);
memory
->
create
(
fp
,
nmax
,
"pair:fp"
);
memory
->
create
(
D_values
,
nmax
,
"pair:D_values"
);
}
double
**
x
=
atom
->
x
;
double
**
f
=
atom
->
f
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
int
newton_pair
=
force
->
newton_pair
;
inum
=
list
->
inum
;
ilist
=
list
->
ilist
;
numneigh
=
list
->
numneigh
;
firstneigh
=
list
->
firstneigh
;
// Zero out per-atom arrays.
int
m
=
nlocal
+
atom
->
nghost
;
for
(
i
=
0
;
i
<
m
;
i
++
)
{
rho
[
i
]
=
0.0
;
rhoB
[
i
]
=
0.0
;
D_values
[
i
]
=
0.0
;
}
// Stage I
// Compute rho and rhoB at each local atom site.
// Additionally calculate the D_i values here if we are using the one-site formulation.
// For the two-site formulation we have to calculate the D values in an extra loop (Stage II).
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
if
(
rsq
<
cutforcesq
)
{
jtype
=
type
[
j
];
double
r
=
sqrt
(
rsq
);
const
EAMTableIndex
index
=
radiusToTableIndex
(
r
);
double
localrho
=
RhoOfR
(
index
,
jtype
,
itype
);
rho
[
i
]
+=
localrho
;
if
(
jtype
==
speciesB
)
rhoB
[
i
]
+=
localrho
;
if
(
newton_pair
||
j
<
nlocal
)
{
localrho
=
RhoOfR
(
index
,
itype
,
jtype
);
rho
[
j
]
+=
localrho
;
if
(
itype
==
speciesB
)
rhoB
[
j
]
+=
localrho
;
}
if
(
cdeamVersion
==
1
&&
itype
!=
jtype
)
{
// Note: if the i-j interaction is not concentration dependent (because either
// i or j are not species A or B) then its contribution to D_i and D_j should
// be ignored.
// This if-clause is only required for a ternary.
if
((
itype
==
speciesA
&&
jtype
==
speciesB
)
||
(
jtype
==
speciesA
&&
itype
==
speciesB
))
{
double
Phi_AB
=
PhiOfR
(
index
,
itype
,
jtype
,
1.0
/
r
);
D_values
[
i
]
+=
Phi_AB
;
if
(
newton_pair
||
j
<
nlocal
)
D_values
[
j
]
+=
Phi_AB
;
}
}
}
}
}
// Communicate and sum densities.
if
(
newton_pair
)
{
communicationStage
=
1
;
comm
->
reverse_comm_pair
(
this
);
}
// fp = derivative of embedding energy at each atom
// phi = embedding energy at each atom
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
EAMTableIndex
index
=
rhoToTableIndex
(
rho
[
i
]);
fp
[
i
]
=
FPrimeOfRho
(
index
,
type
[
i
]);
if
(
eflag
)
{
phi
=
FofRho
(
index
,
type
[
i
]);
if
(
eflag_global
)
eng_vdwl
+=
phi
;
if
(
eflag_atom
)
eatom
[
i
]
+=
phi
;
}
}
// Communicate derivative of embedding function and densities
// and D_values (this for one-site formulation only).
communicationStage
=
2
;
comm
->
forward_comm_pair
(
this
);
// The electron densities may not drop to zero because then the concentration would no longer be defined.
// But the concentration is not needed anyway if there is no interaction with another atom, which is the case
// if the electron density is exactly zero. That's why the following lines have been commented out.
//
//for(i = 0; i < nlocal + atom->nghost; i++) {
// if(rho[i] == 0 && (type[i] == speciesA || type[i] == speciesB))
// error->one(FLERR,"CD-EAM potential routine: Detected atom with zero electron density.");
//}
// Stage II
// This is only required for the original two-site formulation of the CD-EAM potential.
if
(
cdeamVersion
==
2
)
{
// Compute intermediate value D_i for each atom.
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
// This code line is required for ternary alloys.
if
(
itype
!=
speciesA
&&
itype
!=
speciesB
)
continue
;
double
x_i
=
rhoB
[
i
]
/
rho
[
i
];
// Concentration at atom i.
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
jtype
=
type
[
j
];
if
(
itype
==
jtype
)
continue
;
// This code line is required for ternary alloys.
if
(
jtype
!=
speciesA
&&
jtype
!=
speciesB
)
continue
;
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
if
(
rsq
<
cutforcesq
)
{
double
r
=
sqrt
(
rsq
);
const
EAMTableIndex
index
=
radiusToTableIndex
(
r
);
// The concentration independent part of the cross pair potential.
double
Phi_AB
=
PhiOfR
(
index
,
itype
,
jtype
,
1.0
/
r
);
// Average concentration of two sites
double
x_ij
=
0.5
*
(
x_i
+
rhoB
[
j
]
/
rho
[
j
]);
// Calculate derivative of h(x_ij) polynomial function.
double
h_prime
=
evalHprime
(
x_ij
);
D_values
[
i
]
+=
h_prime
*
Phi_AB
/
(
2.0
*
rho
[
i
]
*
rho
[
i
]);
if
(
newton_pair
||
j
<
nlocal
)
D_values
[
j
]
+=
h_prime
*
Phi_AB
/
(
2.0
*
rho
[
j
]
*
rho
[
j
]);
}
}
}
// Communicate and sum D values.
if
(
newton_pair
)
{
communicationStage
=
3
;
comm
->
reverse_comm_pair
(
this
);
}
communicationStage
=
4
;
comm
->
forward_comm_pair
(
this
);
}
// Stage III
// Compute force acting on each atom.
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
// Concentration at site i
double
x_i
=
-
1.0
;
// The value -1 indicates: no concentration dependence for all interactions of atom i.
// It will be replaced by the concentration at site i if atom i is either A or B.
double
D_i
,
h_prime_i
;
// This if-clause is only required for ternary alloys.
if
((
itype
==
speciesA
||
itype
==
speciesB
)
&&
rho
[
i
]
!=
0.0
)
{
// Compute local concentration at site i.
x_i
=
rhoB
[
i
]
/
rho
[
i
];
ASSERT
(
x_i
>=
0
&&
x_i
<=
1.0
);
if
(
cdeamVersion
==
1
)
{
// Calculate derivative of h(x_i) polynomial function.
h_prime_i
=
evalHprime
(
x_i
);
D_i
=
D_values
[
i
]
*
h_prime_i
/
(
2.0
*
rho
[
i
]
*
rho
[
i
]);
}
else
if
(
cdeamVersion
==
2
)
{
D_i
=
D_values
[
i
];
}
else
{
ASSERT
(
false
);
}
}
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
if
(
rsq
<
cutforcesq
)
{
jtype
=
type
[
j
];
double
r
=
sqrt
(
rsq
);
const
EAMTableIndex
index
=
radiusToTableIndex
(
r
);
// rhoip = derivative of (density at atom j due to atom i)
// rhojp = derivative of (density at atom i due to atom j)
// psip needs both fp[i] and fp[j] terms since r_ij appears in two
// terms of embed eng: Fi(sum rho_ij) and Fj(sum rho_ji)
// hence embed' = Fi(sum rho_ij) rhojp + Fj(sum rho_ji) rhoip
rhoip
=
RhoPrimeOfR
(
index
,
itype
,
jtype
);
rhojp
=
RhoPrimeOfR
(
index
,
jtype
,
itype
);
fpair
=
fp
[
i
]
*
rhojp
+
fp
[
j
]
*
rhoip
;
recip
=
1.0
/
r
;
double
x_j
=
-
1
;
// The value -1 indicates: no concentration dependence for this i-j pair
// because atom j is not of species A nor B.
// This code line is required for ternary alloy.
if
(
jtype
==
speciesA
||
jtype
==
speciesB
)
{
ASSERT
(
rho
[
i
]
!=
0.0
);
ASSERT
(
rho
[
j
]
!=
0.0
);
// Compute local concentration at site j.
x_j
=
rhoB
[
j
]
/
rho
[
j
];
ASSERT
(
x_j
>=
0
&&
x_j
<=
1.0
);
double
D_j
=
0.0
;
if
(
cdeamVersion
==
1
)
{
// Calculate derivative of h(x_j) polynomial function.
double
h_prime_j
=
evalHprime
(
x_j
);
D_j
=
D_values
[
j
]
*
h_prime_j
/
(
2.0
*
rho
[
j
]
*
rho
[
j
]);
}
else
if
(
cdeamVersion
==
2
)
{
D_j
=
D_values
[
j
];
}
else
{
ASSERT
(
false
);
}
double
t2
=
-
rhoB
[
j
];
if
(
itype
==
speciesB
)
t2
+=
rho
[
j
];
fpair
+=
D_j
*
rhoip
*
t2
;
}
// This if-clause is only required for a ternary alloy.
// Actually we don't need it at all because D_i should be zero anyway if
// atom i has no concentration dependent interactions (because it is not species A or B).
if
(
x_i
!=
-
1.0
)
{
double
t1
=
-
rhoB
[
i
];
if
(
jtype
==
speciesB
)
t1
+=
rho
[
i
];
fpair
+=
D_i
*
rhojp
*
t1
;
}
double
phip
;
double
phi
=
PhiOfR
(
index
,
itype
,
jtype
,
recip
,
phip
);
if
(
itype
==
jtype
||
x_i
==
-
1.0
||
x_j
==
-
1.0
)
{
// Case of no concentration dependence.
fpair
+=
phip
;
}
else
{
// We have a concentration dependence for the i-j interaction.
double
h
=
0.0
;
if
(
cdeamVersion
==
1
)
{
// Calculate h(x_i) polynomial function.
double
h_i
=
evalH
(
x_i
);
// Calculate h(x_j) polynomial function.
double
h_j
=
evalH
(
x_j
);
h
=
0.5
*
(
h_i
+
h_j
);
}
else
if
(
cdeamVersion
==
2
)
{
// Average concentration.
double
x_ij
=
0.5
*
(
x_i
+
x_j
);
// Calculate h(x_ij) polynomial function.
h
=
evalH
(
x_ij
);
}
else
{
ASSERT
(
false
);
}
fpair
+=
h
*
phip
;
phi
*=
h
;
}
// Divide by r_ij and negate to get forces from gradient.
fpair
/=
-
r
;
f
[
i
][
0
]
+=
delx
*
fpair
;
f
[
i
][
1
]
+=
dely
*
fpair
;
f
[
i
][
2
]
+=
delz
*
fpair
;
if
(
newton_pair
||
j
<
nlocal
)
{
f
[
j
][
0
]
-=
delx
*
fpair
;
f
[
j
][
1
]
-=
dely
*
fpair
;
f
[
j
][
2
]
-=
delz
*
fpair
;
}
if
(
eflag
)
evdwl
=
phi
;
if
(
evflag
)
ev_tally
(
i
,
j
,
nlocal
,
newton_pair
,
evdwl
,
0.0
,
fpair
,
delx
,
dely
,
delz
);
}
}
}
if
(
vflag_fdotr
)
virial_fdotr_compute
();
}
/* ---------------------------------------------------------------------- */
void
PairCDEAM
::
coeff
(
int
narg
,
char
**
arg
)
{
PairEAMAlloy
::
coeff
(
narg
,
arg
);
// Make sure the EAM file is a CD-EAM binary alloy.
if
(
setfl
->
nelements
<
2
)
error
->
all
(
FLERR
,
"The EAM file must contain at least 2 elements to be used with the eam/cd pair style."
);
// Read in the coefficients of the h polynomial from the end of the EAM file.
read_h_coeff
(
arg
[
2
]);
// Determine which atom type is the A species and which is the B species in the alloy.
// By default take the first element (index 0) in the EAM file as the A species
// and the second element (index 1) in the EAM file as the B species.
speciesA
=
-
1
;
speciesB
=
-
1
;
for
(
int
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
{
if
(
map
[
i
]
==
0
)
{
if
(
speciesA
>=
0
)
error
->
all
(
FLERR
,
"The first element from the EAM file may only be mapped to a single atom type."
);
speciesA
=
i
;
}
if
(
map
[
i
]
==
1
)
{
if
(
speciesB
>=
0
)
error
->
all
(
FLERR
,
"The second element from the EAM file may only be mapped to a single atom type."
);
speciesB
=
i
;
}
}
if
(
speciesA
<
0
)
error
->
all
(
FLERR
,
"The first element from the EAM file must be mapped to exactly one atom type."
);
if
(
speciesB
<
0
)
error
->
all
(
FLERR
,
"The second element from the EAM file must be mapped to exactly one atom type."
);
}
/* ----------------------------------------------------------------------
Reads in the h(x) polynomial coefficients
------------------------------------------------------------------------- */
void
PairCDEAM
::
read_h_coeff
(
char
*
filename
)
{
if
(
comm
->
me
==
0
)
{
// Open potential file
FILE
*
fp
;
char
line
[
MAXLINE
];
char
nextline
[
MAXLINE
];
fp
=
force
->
open_potential
(
filename
);
if
(
fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open EAM potential file %s"
,
filename
);
error
->
one
(
FLERR
,
str
);
}
// h coefficients are stored at the end of the file.
// Skip to last line of file.
while
(
fgets
(
nextline
,
MAXLINE
,
fp
)
!=
NULL
)
{
strcpy
(
line
,
nextline
);
}
char
*
ptr
=
strtok
(
line
,
"
\t\n\r\f
"
);
int
degree
=
atoi
(
ptr
);
nhcoeff
=
degree
+
1
;
hcoeff
=
new
double
[
nhcoeff
];
int
i
=
0
;
while
((
ptr
=
strtok
(
NULL
,
"
\t\n\r\f
"
))
!=
NULL
&&
i
<
nhcoeff
)
{
hcoeff
[
i
++
]
=
atof
(
ptr
);
}
if
(
i
!=
nhcoeff
||
nhcoeff
<
1
)
error
->
one
(
FLERR
,
"Failed to read h(x) function coefficients from EAM file."
);
// Close the potential file.
fclose
(
fp
);
}
MPI_Bcast
(
&
nhcoeff
,
1
,
MPI_INT
,
0
,
world
);
if
(
comm
->
me
!=
0
)
hcoeff
=
new
double
[
nhcoeff
];
MPI_Bcast
(
hcoeff
,
nhcoeff
,
MPI_DOUBLE
,
0
,
world
);
}
/* ---------------------------------------------------------------------- */
int
PairCDEAM
::
pack_forward_comm
(
int
n
,
int
*
list
,
double
*
buf
,
int
pbc_flag
,
int
*
pbc
)
{
int
i
,
j
,
m
;
m
=
0
;
if
(
communicationStage
==
2
)
{
if
(
cdeamVersion
==
1
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
fp
[
j
];
buf
[
m
++
]
=
rho
[
j
];
buf
[
m
++
]
=
rhoB
[
j
];
buf
[
m
++
]
=
D_values
[
j
];
}
return
m
;
}
else
if
(
cdeamVersion
==
2
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
fp
[
j
];
buf
[
m
++
]
=
rho
[
j
];
buf
[
m
++
]
=
rhoB
[
j
];
}
return
m
;
}
else
{
ASSERT
(
false
);
return
0
;
}
}
else
if
(
communicationStage
==
4
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
D_values
[
j
];
}
return
m
;
}
else
return
0
;
}
/* ---------------------------------------------------------------------- */
void
PairCDEAM
::
unpack_forward_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
m
=
0
;
last
=
first
+
n
;
if
(
communicationStage
==
2
)
{
if
(
cdeamVersion
==
1
)
{
for
(
i
=
first
;
i
<
last
;
i
++
)
{
fp
[
i
]
=
buf
[
m
++
];
rho
[
i
]
=
buf
[
m
++
];
rhoB
[
i
]
=
buf
[
m
++
];
D_values
[
i
]
=
buf
[
m
++
];
}
}
else
if
(
cdeamVersion
==
2
)
{
for
(
i
=
first
;
i
<
last
;
i
++
)
{
fp
[
i
]
=
buf
[
m
++
];
rho
[
i
]
=
buf
[
m
++
];
rhoB
[
i
]
=
buf
[
m
++
];
}
}
else
{
ASSERT
(
false
);
}
}
else
if
(
communicationStage
==
4
)
{
for
(
i
=
first
;
i
<
last
;
i
++
)
{
D_values
[
i
]
=
buf
[
m
++
];
}
}
}
/* ---------------------------------------------------------------------- */
int
PairCDEAM
::
pack_reverse_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
m
=
0
;
last
=
first
+
n
;
if
(
communicationStage
==
1
)
{
if
(
cdeamVersion
==
1
)
{
for
(
i
=
first
;
i
<
last
;
i
++
)
{
buf
[
m
++
]
=
rho
[
i
];
buf
[
m
++
]
=
rhoB
[
i
];
buf
[
m
++
]
=
D_values
[
i
];
}
return
m
;
}
else
if
(
cdeamVersion
==
2
)
{
for
(
i
=
first
;
i
<
last
;
i
++
)
{
buf
[
m
++
]
=
rho
[
i
];
buf
[
m
++
]
=
rhoB
[
i
];
}
return
m
;
}
else
{
ASSERT
(
false
);
return
0
;
}
}
else
if
(
communicationStage
==
3
)
{
for
(
i
=
first
;
i
<
last
;
i
++
)
{
buf
[
m
++
]
=
D_values
[
i
];
}
return
m
;
}
else
return
0
;
}
/* ---------------------------------------------------------------------- */
void
PairCDEAM
::
unpack_reverse_comm
(
int
n
,
int
*
list
,
double
*
buf
)
{
int
i
,
j
,
m
;
m
=
0
;
if
(
communicationStage
==
1
)
{
if
(
cdeamVersion
==
1
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
rho
[
j
]
+=
buf
[
m
++
];
rhoB
[
j
]
+=
buf
[
m
++
];
D_values
[
j
]
+=
buf
[
m
++
];
}
}
else
if
(
cdeamVersion
==
2
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
rho
[
j
]
+=
buf
[
m
++
];
rhoB
[
j
]
+=
buf
[
m
++
];
}
}
else
{
ASSERT
(
false
);
}
}
else
if
(
communicationStage
==
3
)
{
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
D_values
[
j
]
+=
buf
[
m
++
];
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double
PairCDEAM
::
memory_usage
()
{
double
bytes
=
2
*
nmax
*
sizeof
(
double
);
return
PairEAMAlloy
::
memory_usage
()
+
bytes
;
}
Event Timeline
Log In to Comment