Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81450278
pair_peri_lps_omp.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
Fri, Sep 6, 03:27
Size
12 KB
Mime Type
text/x-c++
Expires
Sun, Sep 8, 03:27 (2 d)
Engine
blob
Format
Raw Data
Handle
20568861
Attached To
rLAMMPS lammps
pair_peri_lps_omp.cpp
View Options
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
This software is distributed under the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "math.h"
#include "float.h"
#include "pair_peri_lps_omp.h"
#include "fix.h"
#include "fix_peri_neigh.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "force.h"
#include "memory.h"
#include "lattice.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_const.h"
#include "suffix.h"
using
namespace
LAMMPS_NS
;
using
namespace
MathConst
;
/* ---------------------------------------------------------------------- */
PairPeriLPSOMP
::
PairPeriLPSOMP
(
LAMMPS
*
lmp
)
:
PairPeriLPS
(
lmp
),
ThrOMP
(
lmp
,
THR_PAIR
)
{
suffix_flag
|=
Suffix
::
OMP
;
respa_enable
=
0
;
}
/* ---------------------------------------------------------------------- */
void
PairPeriLPSOMP
::
compute
(
int
eflag
,
int
vflag
)
{
if
(
eflag
||
vflag
)
{
ev_setup
(
eflag
,
vflag
);
}
else
evflag
=
vflag_fdotr
=
eflag_global
=
eflag_atom
=
0
;
const
int
nall
=
atom
->
nlocal
+
atom
->
nghost
;
const
int
nthreads
=
comm
->
nthreads
;
const
int
inum
=
list
->
inum
;
// grow bond forces array if necessary
if
(
atom
->
nmax
>
nmax
)
{
memory
->
destroy
(
s0_new
);
memory
->
destroy
(
theta
);
nmax
=
atom
->
nmax
;
memory
->
create
(
s0_new
,
nmax
,
"pair:s0_new"
);
memory
->
create
(
theta
,
nmax
,
"pair:theta"
);
}
#if defined(_OPENMP)
#pragma omp parallel default(none) shared(eflag,vflag)
#endif
{
int
ifrom
,
ito
,
tid
;
loop_setup_thr
(
ifrom
,
ito
,
tid
,
inum
,
nthreads
);
ThrData
*
thr
=
fix
->
get_thr
(
tid
);
ev_setup_thr
(
eflag
,
vflag
,
nall
,
eatom
,
vatom
,
thr
);
if
(
evflag
)
{
if
(
eflag
)
{
if
(
force
->
newton_pair
)
eval
<
1
,
1
,
1
>
(
ifrom
,
ito
,
thr
);
else
eval
<
1
,
1
,
0
>
(
ifrom
,
ito
,
thr
);
}
else
{
if
(
force
->
newton_pair
)
eval
<
1
,
0
,
1
>
(
ifrom
,
ito
,
thr
);
else
eval
<
1
,
0
,
0
>
(
ifrom
,
ito
,
thr
);
}
}
else
{
if
(
force
->
newton_pair
)
eval
<
0
,
0
,
1
>
(
ifrom
,
ito
,
thr
);
else
eval
<
0
,
0
,
0
>
(
ifrom
,
ito
,
thr
);
}
reduce_thr
(
this
,
eflag
,
vflag
,
thr
);
}
// end of omp parallel region
}
template
<
int
EVFLAG
,
int
EFLAG
,
int
NEWTON_PAIR
>
void
PairPeriLPSOMP
::
eval
(
int
iifrom
,
int
iito
,
ThrData
*
const
thr
)
{
int
i
,
j
,
ii
,
jj
,
jnum
,
itype
,
jtype
;
double
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
;
double
xtmp0
,
ytmp0
,
ztmp0
,
delx0
,
dely0
,
delz0
,
rsq0
;
double
rsq
,
r
,
dr
,
rk
,
evdwl
,
fpair
,
fbond
;
int
*
ilist
,
*
jlist
,
*
numneigh
,
**
firstneigh
;
double
d_ij
,
delta
,
stretch
;
evdwl
=
0.0
;
const
double
*
const
*
const
x
=
atom
->
x
;
double
*
const
*
const
f
=
thr
->
get_f
();
const
int
*
const
type
=
atom
->
type
;
const
int
nlocal
=
atom
->
nlocal
;
double
fxtmp
,
fytmp
,
fztmp
;
double
*
vfrac
=
atom
->
vfrac
;
double
*
s0
=
atom
->
s0
;
double
**
x0
=
atom
->
x0
;
double
**
r0
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
r0
;
tagint
**
partner
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
partner
;
int
*
npartner
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
npartner
;
double
*
wvolume
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
wvolume
;
// lc = lattice constant
// init_style guarantees it's the same in x, y, and z
double
lc
=
domain
->
lattice
->
xlattice
;
double
half_lc
=
0.5
*
lc
;
double
vfrac_scale
=
1.0
;
// short-range forces
int
periodic
=
(
domain
->
xperiodic
||
domain
->
yperiodic
||
domain
->
zperiodic
);
ilist
=
list
->
ilist
;
numneigh
=
list
->
numneigh
;
firstneigh
=
list
->
firstneigh
;
// loop over neighbors of my atoms
// need minimg() for x0 difference since not ghosted
for
(
ii
=
iifrom
;
ii
<
iito
;
++
ii
)
{
i
=
ilist
[
ii
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
xtmp0
=
x0
[
i
][
0
];
ytmp0
=
x0
[
i
][
1
];
ztmp0
=
x0
[
i
][
2
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
fxtmp
=
fytmp
=
fztmp
=
0.0
;
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
;
delx0
=
xtmp0
-
x0
[
j
][
0
];
dely0
=
ytmp0
-
x0
[
j
][
1
];
delz0
=
ztmp0
-
x0
[
j
][
2
];
if
(
periodic
)
domain
->
minimum_image
(
delx0
,
dely0
,
delz0
);
rsq0
=
delx0
*
delx0
+
dely0
*
dely0
+
delz0
*
delz0
;
jtype
=
type
[
j
];
r
=
sqrt
(
rsq
);
// short-range interaction distance based on initial particle position
// 0.9 and 1.35 are constants
d_ij
=
MIN
(
0.9
*
sqrt
(
rsq0
),
1.35
*
lc
);
// short-range contact forces
// 15 is constant taken from the EMU Theory Manual
// Silling, 12 May 2005, p 18
if
(
r
<
d_ij
)
{
dr
=
r
-
d_ij
;
// kshort based upon short-range force constant
// of the bond-based theory used in PMB model
double
kshort
=
(
15.0
*
18.0
*
bulkmodulus
[
itype
][
itype
])
/
(
MY_PI
*
cutsq
[
itype
][
jtype
]
*
cutsq
[
itype
][
jtype
]);
rk
=
(
kshort
*
vfrac
[
j
])
*
(
dr
/
cut
[
itype
][
jtype
]);
if
(
r
>
0.0
)
fpair
=
-
(
rk
/
r
);
else
fpair
=
0.0
;
fxtmp
+=
delx
*
fpair
;
fytmp
+=
dely
*
fpair
;
fztmp
+=
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
=
0.5
*
rk
*
dr
;
if
(
EVFLAG
)
ev_tally_thr
(
this
,
i
,
j
,
nlocal
,
NEWTON_PAIR
,
evdwl
,
0.0
,
fpair
*
vfrac
[
i
],
delx
,
dely
,
delz
,
thr
);
}
}
f
[
i
][
0
]
+=
fxtmp
;
f
[
i
][
1
]
+=
fytmp
;
f
[
i
][
2
]
+=
fztmp
;
}
// wait until all threads are done since we
// need to distribute the work differently.
sync_threads
();
#if defined(_OPENMP)
// each thread works on a fixed chunk of atoms.
const
int
idelta
=
1
+
nlocal
/
comm
->
nthreads
;
iifrom
=
thr
->
get_tid
()
*
idelta
;
iito
=
((
iifrom
+
idelta
)
>
nlocal
)
?
nlocal
:
(
iifrom
+
idelta
);
#else
iifrom
=
0
;
iito
=
nlocal
;
#endif
// Compute the dilatation on each particle
if
(
iifrom
<
nlocal
)
compute_dilatation_thr
(
iifrom
,
iito
);
// wait until all threads are done before communication
sync_threads
();
#if defined(_OPENMP)
#pragma omp master
#endif
{
// communicate dilatation (theta) of each particle
comm
->
forward_comm_pair
(
this
);
// communicate weighted volume (wvolume) upon every reneighbor
if
(
neighbor
->
ago
==
0
)
comm
->
forward_comm_fix
(
modify
->
fix
[
ifix_peri
]);
}
sync_threads
();
// Volume-dependent part of the energy
if
(
EFLAG
)
{
for
(
i
=
iifrom
;
i
<
iito
;
i
++
)
{
itype
=
type
[
i
];
e_tally_thr
(
this
,
i
,
i
,
nlocal
,
NEWTON_PAIR
,
0.5
*
bulkmodulus
[
itype
][
itype
]
*
(
theta
[
i
]
*
theta
[
i
]),
0.0
,
thr
);
}
}
// loop over my particles and their partners
// partner list contains all bond partners, so I-J appears twice
// if bond already broken, skip this partner
// first = true if this is first neighbor of particle i
bool
first
;
double
omega_minus
,
omega_plus
;
for
(
i
=
iifrom
;
i
<
iito
;
++
i
)
{
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
xtmp0
=
x0
[
i
][
0
];
ytmp0
=
x0
[
i
][
1
];
ztmp0
=
x0
[
i
][
2
];
itype
=
type
[
i
];
jnum
=
npartner
[
i
];
first
=
true
;
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
if
(
partner
[
i
][
jj
]
==
0
)
continue
;
j
=
atom
->
map
(
partner
[
i
][
jj
]);
// check if lost a partner without first breaking bond
if
(
j
<
0
)
{
partner
[
i
][
jj
]
=
0
;
continue
;
}
// compute force density, add to PD equation of motion
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
if
(
periodic
)
domain
->
minimum_image
(
delx
,
dely
,
delz
);
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
delx0
=
xtmp0
-
x0
[
j
][
0
];
dely0
=
ytmp0
-
x0
[
j
][
1
];
delz0
=
ztmp0
-
x0
[
j
][
2
];
if
(
periodic
)
domain
->
minimum_image
(
delx0
,
dely0
,
delz0
);
jtype
=
type
[
j
];
delta
=
cut
[
itype
][
jtype
];
r
=
sqrt
(
rsq
);
dr
=
r
-
r0
[
i
][
jj
];
// avoid roundoff errors
if
(
fabs
(
dr
)
<
2.2204e-016
)
dr
=
0.0
;
// scale vfrac[j] if particle j near the horizon
if
((
fabs
(
r0
[
i
][
jj
]
-
delta
))
<=
half_lc
)
vfrac_scale
=
(
-
1.0
/
(
2
*
half_lc
))
*
(
r0
[
i
][
jj
])
+
(
1.0
+
((
delta
-
half_lc
)
/
(
2
*
half_lc
)
)
);
else
vfrac_scale
=
1.0
;
omega_plus
=
influence_function
(
-
1.0
*
delx0
,
-
1.0
*
dely0
,
-
1.0
*
delz0
);
omega_minus
=
influence_function
(
delx0
,
dely0
,
delz0
);
rk
=
(
(
3.0
*
bulkmodulus
[
itype
][
itype
])
-
(
5.0
*
shearmodulus
[
itype
][
itype
])
)
*
vfrac
[
j
]
*
vfrac_scale
*
(
(
omega_plus
*
theta
[
i
]
/
wvolume
[
i
])
+
(
omega_minus
*
theta
[
j
]
/
wvolume
[
j
]
)
)
*
r0
[
i
][
jj
];
rk
+=
15.0
*
(
shearmodulus
[
itype
][
itype
]
*
vfrac
[
j
]
*
vfrac_scale
)
*
(
(
omega_plus
/
wvolume
[
i
])
+
(
omega_minus
/
wvolume
[
j
])
)
*
dr
;
if
(
r
>
0.0
)
fbond
=
-
(
rk
/
r
);
else
fbond
=
0.0
;
f
[
i
][
0
]
+=
delx
*
fbond
;
f
[
i
][
1
]
+=
dely
*
fbond
;
f
[
i
][
2
]
+=
delz
*
fbond
;
// since I-J is double counted, set newton off & use 1/2 factor and I,I
double
deviatoric_extension
=
dr
-
(
theta
[
i
]
*
r0
[
i
][
jj
]
/
3.0
);
if
(
EFLAG
)
evdwl
=
0.5
*
15
*
(
shearmodulus
[
itype
][
itype
]
/
wvolume
[
i
])
*
omega_plus
*
(
deviatoric_extension
*
deviatoric_extension
)
*
vfrac
[
j
]
*
vfrac_scale
;
if
(
EVFLAG
)
ev_tally_thr
(
this
,
i
,
i
,
nlocal
,
0
,
0.5
*
evdwl
,
0.0
,
0.5
*
fbond
*
vfrac
[
i
],
delx
,
dely
,
delz
,
thr
);
// find stretch in bond I-J and break if necessary
// use s0 from previous timestep
stretch
=
dr
/
r0
[
i
][
jj
];
if
(
stretch
>
MIN
(
s0
[
i
],
s0
[
j
]))
partner
[
i
][
jj
]
=
0
;
// update s0 for next timestep
if
(
first
)
s0_new
[
i
]
=
s00
[
itype
][
jtype
]
-
(
alpha
[
itype
][
jtype
]
*
stretch
);
else
s0_new
[
i
]
=
MAX
(
s0_new
[
i
],
s00
[
itype
][
jtype
]
-
(
alpha
[
itype
][
jtype
]
*
stretch
));
first
=
false
;
}
}
sync_threads
();
// store new s0 (in parallel)
if
(
iifrom
<
nlocal
)
for
(
i
=
iifrom
;
i
<
iito
;
i
++
)
s0
[
i
]
=
s0_new
[
i
];
}
/* ---------------------------------------------------------------------- */
void
PairPeriLPSOMP
::
compute_dilatation_thr
(
int
ifrom
,
int
ito
)
{
int
i
,
j
,
jj
,
jnum
,
itype
,
jtype
;
double
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
;
double
xtmp0
,
ytmp0
,
ztmp0
,
delx0
,
dely0
,
delz0
;
double
rsq
,
r
,
dr
;
double
delta
;
double
**
x
=
atom
->
x
;
int
*
type
=
atom
->
type
;
double
**
x0
=
atom
->
x0
;
double
*
vfrac
=
atom
->
vfrac
;
double
vfrac_scale
=
1.0
;
double
lc
=
domain
->
lattice
->
xlattice
;
double
half_lc
=
0.5
*
lc
;
double
**
r0
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
r0
;
tagint
**
partner
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
partner
;
int
*
npartner
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
npartner
;
double
*
wvolume
=
((
FixPeriNeigh
*
)
modify
->
fix
[
ifix_peri
])
->
wvolume
;
int
periodic
=
domain
->
xperiodic
||
domain
->
yperiodic
||
domain
->
zperiodic
;
// compute the dilatation theta
for
(
i
=
ifrom
;
i
<
ito
;
i
++
)
{
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
xtmp0
=
x0
[
i
][
0
];
ytmp0
=
x0
[
i
][
1
];
ztmp0
=
x0
[
i
][
2
];
jnum
=
npartner
[
i
];
theta
[
i
]
=
0.0
;
itype
=
type
[
i
];
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
// if bond already broken, skip this partner
if
(
partner
[
i
][
jj
]
==
0
)
continue
;
// Look up local index of this partner particle
j
=
atom
->
map
(
partner
[
i
][
jj
]);
// Skip if particle is "lost"
if
(
j
<
0
)
continue
;
// Compute force density and add to PD equation of motion
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
if
(
periodic
)
domain
->
minimum_image
(
delx
,
dely
,
delz
);
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
delx0
=
xtmp0
-
x0
[
j
][
0
];
dely0
=
ytmp0
-
x0
[
j
][
1
];
delz0
=
ztmp0
-
x0
[
j
][
2
];
if
(
periodic
)
domain
->
minimum_image
(
delx0
,
dely0
,
delz0
);
r
=
sqrt
(
rsq
);
dr
=
r
-
r0
[
i
][
jj
];
if
(
fabs
(
dr
)
<
2.2204e-016
)
dr
=
0.0
;
jtype
=
type
[
j
];
delta
=
cut
[
itype
][
jtype
];
// scale vfrac[j] if particle j near the horizon
if
((
fabs
(
r0
[
i
][
jj
]
-
delta
))
<=
half_lc
)
vfrac_scale
=
(
-
1.0
/
(
2
*
half_lc
))
*
(
r0
[
i
][
jj
])
+
(
1.0
+
((
delta
-
half_lc
)
/
(
2
*
half_lc
)
)
);
else
vfrac_scale
=
1.0
;
theta
[
i
]
+=
influence_function
(
delx0
,
dely0
,
delz0
)
*
r0
[
i
][
jj
]
*
dr
*
vfrac
[
j
]
*
vfrac_scale
;
}
// if wvolume[i] is zero, then particle i has no bonds
// therefore, the dilatation is set to
if
(
wvolume
[
i
]
!=
0.0
)
theta
[
i
]
=
(
3.0
/
wvolume
[
i
])
*
theta
[
i
];
else
theta
[
i
]
=
0
;
}
}
/* ---------------------------------------------------------------------- */
double
PairPeriLPSOMP
::
memory_usage
()
{
double
bytes
=
memory_usage_thr
();
bytes
+=
PairPeriLPS
::
memory_usage
();
return
bytes
;
}
Event Timeline
Log In to Comment