Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69986056
npair_half_bin_newton_ssa.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
Thu, Jul 4, 15:53
Size
10 KB
Mime Type
text/x-c
Expires
Sat, Jul 6, 15:53 (2 d)
Engine
blob
Format
Raw Data
Handle
18768707
Attached To
rLAMMPS lammps
npair_half_bin_newton_ssa.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 authors:
James Larentzos and Timothy I. Mattox (Engility Corporation)
------------------------------------------------------------------------- */
#include "npair_half_bin_newton_ssa.h"
#include "neighbor.h"
#include "nstencil_ssa.h"
#include "nbin_ssa.h"
#include "neigh_list.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "domain.h"
#include "group.h"
#include "memory.h"
#include "my_page.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
/* ---------------------------------------------------------------------- */
NPairHalfBinNewtonSSA
::
NPairHalfBinNewtonSSA
(
LAMMPS
*
lmp
)
:
NPair
(
lmp
)
{
ssa_maxPhaseCt
=
0
;
ssa_maxPhaseLen
=
0
;
ssa_phaseCt
=
0
;
ssa_phaseLen
=
NULL
;
ssa_itemLoc
=
NULL
;
ssa_itemLen
=
NULL
;
}
/* ---------------------------------------------------------------------- */
NPairHalfBinNewtonSSA
::~
NPairHalfBinNewtonSSA
()
{
ssa_maxPhaseCt
=
0
;
ssa_maxPhaseLen
=
0
;
ssa_phaseCt
=
0
;
memory
->
destroy
(
ssa_phaseLen
);
memory
->
destroy
(
ssa_itemLoc
);
memory
->
destroy
(
ssa_itemLen
);
}
/* ----------------------------------------------------------------------
binned neighbor list construction with full Newton's 3rd law
for use by Shardlow Spliting Algorithm
each owned atom i checks its own bin and other bins in Newton stencil
every pair stored exactly once by some processor
------------------------------------------------------------------------- */
void
NPairHalfBinNewtonSSA
::
build
(
NeighList
*
list
)
{
int
i
,
j
,
k
,
n
,
itype
,
jtype
,
ibin
,
which
,
imol
,
iatom
,
moltemplate
;
tagint
tagprev
;
double
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
,
rsq
;
int
*
neighptr
;
double
**
x
=
atom
->
x
;
int
*
type
=
atom
->
type
;
int
*
mask
=
atom
->
mask
;
tagint
*
tag
=
atom
->
tag
;
tagint
*
molecule
=
atom
->
molecule
;
tagint
**
special
=
atom
->
special
;
int
**
nspecial
=
atom
->
nspecial
;
int
nlocal
=
atom
->
nlocal
;
if
(
includegroup
)
nlocal
=
atom
->
nfirst
;
int
*
molindex
=
atom
->
molindex
;
int
*
molatom
=
atom
->
molatom
;
Molecule
**
onemols
=
atom
->
avec
->
onemols
;
int
molecular
=
atom
->
molecular
;
if
(
molecular
==
2
)
moltemplate
=
1
;
else
moltemplate
=
0
;
int
*
ilist
=
list
->
ilist
;
int
*
numneigh
=
list
->
numneigh
;
int
**
firstneigh
=
list
->
firstneigh
;
MyPage
<
int
>
*
ipage
=
list
->
ipage
;
NStencilSSA
*
ns_ssa
=
dynamic_cast
<
NStencilSSA
*>
(
ns
);
if
(
!
ns_ssa
)
error
->
one
(
FLERR
,
"NStencil wasn't a NStencilSSA object"
);
int
*
nstencil_ssa
=
&
(
ns_ssa
->
nstencil_ssa
[
0
]);
int
nstencil_full
=
ns_ssa
->
nstencil
;
NBinSSA
*
nb_ssa
=
dynamic_cast
<
NBinSSA
*>
(
nb
);
if
(
!
nb_ssa
)
error
->
one
(
FLERR
,
"NBin wasn't a NBinSSA object"
);
int
*
bins
=
nb_ssa
->
bins
;
int
*
binhead
=
nb_ssa
->
binhead
;
int
*
gairhead_ssa
=
&
(
nb_ssa
->
gairhead_ssa
[
0
]);
int
inum
=
0
;
int
gnum
=
0
;
int
xbin
,
ybin
,
zbin
,
xbin2
,
ybin2
,
zbin2
;
int
**
stencilxyz
=
ns_ssa
->
stencilxyz
;
int
lbinxlo
=
nb_ssa
->
lbinxlo
;
int
lbinxhi
=
nb_ssa
->
lbinxhi
;
int
lbinylo
=
nb_ssa
->
lbinylo
;
int
lbinyhi
=
nb_ssa
->
lbinyhi
;
int
lbinzlo
=
nb_ssa
->
lbinzlo
;
int
lbinzhi
=
nb_ssa
->
lbinzhi
;
int
sx1
=
ns_ssa
->
sx
+
1
;
int
sy1
=
ns_ssa
->
sy
+
1
;
int
sz1
=
ns_ssa
->
sz
+
1
;
ssa_phaseCt
=
sz1
*
sy1
*
sx1
;
xbin
=
(
lbinxhi
-
lbinxlo
+
sx1
-
1
)
/
sx1
+
1
;
ybin
=
(
lbinyhi
-
lbinylo
+
sy1
-
1
)
/
sy1
+
1
;
zbin
=
(
lbinzhi
-
lbinzlo
+
sz1
-
1
)
/
sz1
+
1
;
int
phaseLenEstimate
=
xbin
*
ybin
*
zbin
;
if
(
ssa_phaseCt
>
ssa_maxPhaseCt
)
{
ssa_maxPhaseCt
=
ssa_phaseCt
;
ssa_maxPhaseLen
=
0
;
memory
->
destroy
(
ssa_phaseLen
);
memory
->
destroy
(
ssa_itemLoc
);
memory
->
destroy
(
ssa_itemLen
);
memory
->
create
(
ssa_phaseLen
,
ssa_maxPhaseCt
,
"NPairHalfBinNewtonSSA:ssa_phaseLen"
);
}
if
(
phaseLenEstimate
>
ssa_maxPhaseLen
)
{
ssa_maxPhaseLen
=
phaseLenEstimate
;
memory
->
destroy
(
ssa_itemLoc
);
memory
->
destroy
(
ssa_itemLen
);
memory
->
create
(
ssa_itemLoc
,
ssa_maxPhaseCt
,
ssa_maxPhaseLen
,
"NPairHalfBinNewtonSSA:ssa_itemLoc"
);
memory
->
create
(
ssa_itemLen
,
ssa_maxPhaseCt
,
ssa_maxPhaseLen
,
"NPairHalfBinNewtonSSA:ssa_itemLen"
);
}
ipage
->
reset
();
int
workPhase
=
0
;
// loop over bins with local atoms, storing half of the neighbors
for
(
int
zoff
=
ns_ssa
->
sz
;
zoff
>=
0
;
--
zoff
)
{
for
(
int
yoff
=
ns_ssa
->
sy
;
yoff
>=
0
;
--
yoff
)
{
for
(
int
xoff
=
ns_ssa
->
sx
;
xoff
>=
0
;
--
xoff
)
{
int
workItem
=
0
;
for
(
zbin
=
lbinzlo
+
zoff
;
zbin
<
lbinzhi
;
zbin
+=
sz1
)
{
for
(
ybin
=
lbinylo
+
yoff
-
ns_ssa
->
sy
;
ybin
<
lbinyhi
;
ybin
+=
sy1
)
{
for
(
xbin
=
lbinxlo
+
xoff
-
ns_ssa
->
sx
;
xbin
<
lbinxhi
;
xbin
+=
sx1
)
{
if
(
workItem
>=
phaseLenEstimate
)
error
->
one
(
FLERR
,
"phaseLenEstimate was too small"
);
ssa_itemLoc
[
workPhase
][
workItem
]
=
inum
;
// record where workItem starts in ilist
for
(
int
subphase
=
0
;
subphase
<
4
;
subphase
++
)
{
int
s_ybin
=
ybin
+
((
subphase
&
0x2
)
?
ns_ssa
->
sy
:
0
);
int
s_xbin
=
xbin
+
((
subphase
&
0x1
)
?
ns_ssa
->
sx
:
0
);
int
ibin
,
ct
;
if
((
s_ybin
<
lbinylo
)
||
(
s_ybin
>=
lbinyhi
))
continue
;
if
((
s_xbin
<
lbinxlo
)
||
(
s_xbin
>=
lbinxhi
))
continue
;
ibin
=
zbin
*
nb_ssa
->
mbiny
*
nb_ssa
->
mbinx
+
s_ybin
*
nb_ssa
->
mbinx
+
s_xbin
;
for
(
i
=
binhead
[
ibin
];
i
>=
0
;
i
=
bins
[
i
])
{
n
=
0
;
neighptr
=
ipage
->
vget
();
itype
=
type
[
i
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
if
(
moltemplate
)
{
imol
=
molindex
[
i
];
iatom
=
molatom
[
i
];
tagprev
=
tag
[
i
]
-
iatom
-
1
;
}
// loop over all local atoms in the current stencil "subphase"
for
(
k
=
nstencil_ssa
[
subphase
];
k
<
nstencil_ssa
[
subphase
+
1
];
k
++
)
{
const
int
jbin
=
ibin
+
stencil
[
k
];
if
(
jbin
!=
ibin
)
j
=
binhead
[
jbin
];
else
j
=
bins
[
i
];
// same bin as i, so start just past i in the bin
for
(;
j
>=
0
;
j
=
bins
[
j
])
{
jtype
=
type
[
j
];
if
(
exclude
&&
exclusion
(
i
,
j
,
itype
,
jtype
,
mask
,
molecule
))
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
<=
cutneighsq
[
itype
][
jtype
])
{
if
(
molecular
)
{
if
(
!
moltemplate
)
which
=
find_special
(
special
[
i
],
nspecial
[
i
],
tag
[
j
]);
else
if
(
imol
>=
0
)
which
=
find_special
(
onemols
[
imol
]
->
special
[
iatom
],
onemols
[
imol
]
->
nspecial
[
iatom
],
tag
[
j
]
-
tagprev
);
else
which
=
0
;
if
(
which
==
0
)
neighptr
[
n
++
]
=
j
;
else
if
(
domain
->
minimum_image_check
(
delx
,
dely
,
delz
))
neighptr
[
n
++
]
=
j
;
else
if
(
which
>
0
)
neighptr
[
n
++
]
=
j
^
(
which
<<
SBBITS
);
}
else
neighptr
[
n
++
]
=
j
;
}
}
}
if
(
n
>
0
)
{
firstneigh
[
inum
]
=
neighptr
;
numneigh
[
inum
]
=
n
;
ilist
[
inum
++
]
=
i
;
}
ipage
->
vgot
(
n
);
if
(
ipage
->
status
())
error
->
one
(
FLERR
,
"Neighbor list overflow, boost neigh_modify one"
);
}
}
// record where workItem ends in ilist
ssa_itemLen
[
workPhase
][
workItem
]
=
inum
-
ssa_itemLoc
[
workPhase
][
workItem
];
if
(
ssa_itemLen
[
workPhase
][
workItem
]
>
0
)
workItem
++
;
}
}
}
// record where workPhase ends
ssa_phaseLen
[
workPhase
++
]
=
workItem
;
}
}
}
if
(
ssa_phaseCt
!=
workPhase
)
error
->
one
(
FLERR
,
"ssa_phaseCt was wrong"
);
list
->
AIRct_ssa
[
0
]
=
list
->
inum
=
inum
;
// loop over AIR ghost atoms, storing their local neighbors
// since these are ghosts, must check if stencil bin is out of bounds
for
(
int
airnum
=
1
;
airnum
<=
7
;
airnum
++
)
{
int
locAIRct
=
0
;
for
(
i
=
gairhead_ssa
[
airnum
];
i
>=
0
;
i
=
bins
[
i
])
{
n
=
0
;
neighptr
=
ipage
->
vget
();
itype
=
type
[
i
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
ibin
=
coord2bin
(
x
[
i
],
xbin
,
ybin
,
zbin
);
// loop over AIR ghost atoms in all bins in "full" stencil
// Note: the non-AIR ghost atoms have already been filtered out
for
(
k
=
0
;
k
<
nstencil_full
;
k
++
)
{
xbin2
=
xbin
+
stencilxyz
[
k
][
0
];
ybin2
=
ybin
+
stencilxyz
[
k
][
1
];
zbin2
=
zbin
+
stencilxyz
[
k
][
2
];
// Skip it if this bin is outside the extent of local bins
if
(
xbin2
<
lbinxlo
||
xbin2
>=
lbinxhi
||
ybin2
<
lbinylo
||
ybin2
>=
lbinyhi
||
zbin2
<
lbinzlo
||
zbin2
>=
lbinzhi
)
continue
;
for
(
j
=
binhead
[
ibin
+
stencil
[
k
]];
j
>=
0
;
j
=
bins
[
j
])
{
jtype
=
type
[
j
];
if
(
exclude
&&
exclusion
(
i
,
j
,
itype
,
jtype
,
mask
,
molecule
))
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
<=
cutneighsq
[
itype
][
jtype
])
{
if
(
molecular
)
{
if
(
!
moltemplate
)
which
=
find_special
(
special
[
j
],
nspecial
[
j
],
tag
[
i
]);
else
{
int
jmol
=
molindex
[
j
];
if
(
jmol
>=
0
)
{
int
jatom
=
molatom
[
j
];
which
=
find_special
(
onemols
[
jmol
]
->
special
[
jatom
],
onemols
[
jmol
]
->
nspecial
[
jatom
],
tag
[
i
]
-
(
tag
[
j
]
-
jatom
-
1
));
}
else
which
=
0
;
}
if
(
which
==
0
)
neighptr
[
n
++
]
=
j
;
else
if
(
domain
->
minimum_image_check
(
delx
,
dely
,
delz
))
neighptr
[
n
++
]
=
j
;
else
if
(
which
>
0
)
neighptr
[
n
++
]
=
j
^
(
which
<<
SBBITS
);
}
else
neighptr
[
n
++
]
=
j
;
}
}
}
if
(
n
>
0
)
{
firstneigh
[
inum
+
gnum
]
=
neighptr
;
numneigh
[
inum
+
gnum
]
=
n
;
ilist
[
inum
+
(
gnum
++
)]
=
i
;
++
locAIRct
;
}
ipage
->
vgot
(
n
);
if
(
ipage
->
status
())
error
->
one
(
FLERR
,
"Neighbor (ghost) list overflow, boost neigh_modify one"
);
}
list
->
AIRct_ssa
[
airnum
]
=
locAIRct
;
}
list
->
gnum
=
gnum
;
}
Event Timeline
Log In to Comment