Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82376014
pppm_intel.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
Wed, Sep 11, 01:47
Size
20 KB
Mime Type
text/x-c
Expires
Fri, Sep 13, 01:47 (2 d)
Engine
blob
Format
Raw Data
Handle
20691526
Attached To
rLAMMPS lammps
pppm_intel.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: Rodrigo Canales (RWTH Aachen University)
W. Michael Brown (Intel)
------------------------------------------------------------------------- */
#include <mpi.h>
#include <stdlib.h>
#include <math.h>
#include "pppm_intel.h"
#include "atom.h"
#include "error.h"
#include "gridcomm.h"
#include "math_const.h"
#include "math_special.h"
#include "memory.h"
#include "suffix.h"
using
namespace
LAMMPS_NS
;
using
namespace
MathConst
;
using
namespace
MathSpecial
;
#define MAXORDER 7
#define OFFSET 16384
#define LARGE 10000.0
#define SMALL 0.00001
#define EPS_HOC 1.0e-7
enum
{
REVERSE_RHO
};
enum
{
FORWARD_IK
,
FORWARD_AD
,
FORWARD_IK_PERATOM
,
FORWARD_AD_PERATOM
};
#ifdef FFT_SINGLE
#define ZEROF 0.0f
#define ONEF 1.0f
#else
#define ZEROF 0.0
#define ONEF 1.0
#endif
/* ---------------------------------------------------------------------- */
PPPMIntel
::
PPPMIntel
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
PPPM
(
lmp
,
narg
,
arg
)
{
suffix_flag
|=
Suffix
::
INTEL
;
}
PPPMIntel
::~
PPPMIntel
()
{
}
/* ----------------------------------------------------------------------
called once before run
------------------------------------------------------------------------- */
void
PPPMIntel
::
init
()
{
PPPM
::
init
();
int
ifix
=
modify
->
find_fix
(
"package_intel"
);
if
(
ifix
<
0
)
error
->
all
(
FLERR
,
"The 'package intel' command is required for /intel styles"
);
fix
=
static_cast
<
FixIntel
*>
(
modify
->
fix
[
ifix
]);
#ifdef _LMP_INTEL_OFFLOAD
_use_base
=
0
;
if
(
fix
->
offload_balance
()
!=
0.0
)
{
_use_base
=
1
;
return
;
}
#endif
fix
->
kspace_init_check
();
if
(
order
>
INTEL_P3M_MAXORDER
)
error
->
all
(
FLERR
,
"PPPM order greater than supported by USER-INTEL
\n
"
);
/*
if (fix->precision() == FixIntel::PREC_MODE_MIXED)
pack_force_const(force_const_single, fix->get_mixed_buffers());
else if (fix->precision() == FixIntel::PREC_MODE_DOUBLE)
pack_force_const(force_const_double, fix->get_double_buffers());
else
pack_force_const(force_const_single, fix->get_single_buffers());
*/
}
/* ----------------------------------------------------------------------
compute the PPPMIntel long-range force, energy, virial
------------------------------------------------------------------------- */
void
PPPMIntel
::
compute
(
int
eflag
,
int
vflag
)
{
#ifdef _LMP_INTEL_OFFLOAD
if
(
_use_base
)
{
PPPM
::
compute
(
eflag
,
vflag
);
return
;
}
#endif
compute_first
(
eflag
,
vflag
);
compute_second
(
eflag
,
vflag
);
}
/* ---------------------------------------------------------------------- */
void
PPPMIntel
::
compute_first
(
int
eflag
,
int
vflag
)
{
int
i
,
j
;
// set energy/virial flags
// invoke allocate_peratom() if needed for first time
if
(
eflag
||
vflag
)
ev_setup
(
eflag
,
vflag
);
else
evflag
=
evflag_atom
=
eflag_global
=
vflag_global
=
eflag_atom
=
vflag_atom
=
0
;
if
(
evflag_atom
&&
!
peratom_allocate_flag
)
{
allocate_peratom
();
cg_peratom
->
ghost_notify
();
cg_peratom
->
setup
();
}
// if atom count has changed, update qsum and qsqsum
if
(
atom
->
natoms
!=
natoms_original
)
{
qsum_qsq
();
natoms_original
=
atom
->
natoms
;
}
// return if there are no charges
if
(
qsqsum
==
0.0
)
return
;
// convert atoms from box to lamda coords
if
(
triclinic
==
0
)
boxlo
=
domain
->
boxlo
;
else
{
boxlo
=
domain
->
boxlo_lamda
;
domain
->
x2lamda
(
atom
->
nlocal
);
}
// extend size of per-atom arrays if necessary
if
(
atom
->
nmax
>
nmax
)
{
memory
->
destroy
(
part2grid
);
nmax
=
atom
->
nmax
;
memory
->
create
(
part2grid
,
nmax
,
3
,
"pppm:part2grid"
);
}
// find grid points for all my particles
// map my particle charge onto my local 3d density grid
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_MIXED
)
{
particle_map
<
float
,
double
>
(
fix
->
get_mixed_buffers
());
make_rho
<
float
,
double
>
(
fix
->
get_mixed_buffers
());
}
else
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_DOUBLE
)
{
particle_map
<
double
,
double
>
(
fix
->
get_double_buffers
());
make_rho
<
double
,
double
>
(
fix
->
get_double_buffers
());
}
else
{
particle_map
<
float
,
float
>
(
fix
->
get_single_buffers
());
make_rho
<
float
,
float
>
(
fix
->
get_single_buffers
());
}
// all procs communicate density values from their ghost cells
// to fully sum contribution in their 3d bricks
// remap from 3d decomposition to FFT decomposition
cg
->
reverse_comm
(
this
,
REVERSE_RHO
);
brick2fft
();
// compute potential gradient on my FFT grid and
// portion of e_long on this proc's FFT grid
// return gradients (electric fields) in 3d brick decomposition
// also performs per-atom calculations via poisson_peratom()
poisson
();
// all procs communicate E-field values
// to fill ghost cells surrounding their 3d bricks
if
(
differentiation_flag
==
1
)
cg
->
forward_comm
(
this
,
FORWARD_AD
);
else
cg
->
forward_comm
(
this
,
FORWARD_IK
);
// extra per-atom energy/virial communication
if
(
evflag_atom
)
{
if
(
differentiation_flag
==
1
&&
vflag_atom
)
cg_peratom
->
forward_comm
(
this
,
FORWARD_AD_PERATOM
);
else
if
(
differentiation_flag
==
0
)
cg_peratom
->
forward_comm
(
this
,
FORWARD_IK_PERATOM
);
}
}
/* ---------------------------------------------------------------------- */
void
PPPMIntel
::
compute_second
(
int
eflag
,
int
vflag
)
{
int
i
,
j
;
// calculate the force on my particles
if
(
differentiation_flag
==
1
)
{
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_MIXED
)
fieldforce_ad
<
float
,
double
>
(
fix
->
get_mixed_buffers
());
else
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_DOUBLE
)
fieldforce_ad
<
double
,
double
>
(
fix
->
get_double_buffers
());
else
fieldforce_ad
<
float
,
float
>
(
fix
->
get_single_buffers
());
}
else
{
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_MIXED
)
fieldforce_ik
<
float
,
double
>
(
fix
->
get_mixed_buffers
());
else
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_DOUBLE
)
fieldforce_ik
<
double
,
double
>
(
fix
->
get_double_buffers
());
else
fieldforce_ik
<
float
,
float
>
(
fix
->
get_single_buffers
());
}
// extra per-atom energy/virial communication
if
(
evflag_atom
)
fieldforce_peratom
();
// sum global energy across procs and add in volume-dependent term
const
double
qscale
=
qqrd2e
*
scale
;
if
(
eflag_global
)
{
double
energy_all
;
MPI_Allreduce
(
&
energy
,
&
energy_all
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
energy
=
energy_all
;
energy
*=
0.5
*
volume
;
energy
-=
g_ewald
*
qsqsum
/
MY_PIS
+
MY_PI2
*
qsum
*
qsum
/
(
g_ewald
*
g_ewald
*
volume
);
energy
*=
qscale
;
}
// sum global virial across procs
if
(
vflag_global
)
{
double
virial_all
[
6
];
MPI_Allreduce
(
virial
,
virial_all
,
6
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
for
(
i
=
0
;
i
<
6
;
i
++
)
virial
[
i
]
=
0.5
*
qscale
*
volume
*
virial_all
[
i
];
}
// per-atom energy/virial
// energy includes self-energy correction
// ntotal accounts for TIP4P tallying eatom/vatom for ghost atoms
if
(
evflag_atom
)
{
double
*
q
=
atom
->
q
;
int
nlocal
=
atom
->
nlocal
;
int
ntotal
=
nlocal
;
if
(
tip4pflag
)
ntotal
+=
atom
->
nghost
;
if
(
eflag_atom
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
eatom
[
i
]
*=
0.5
;
eatom
[
i
]
-=
g_ewald
*
q
[
i
]
*
q
[
i
]
/
MY_PIS
+
MY_PI2
*
q
[
i
]
*
qsum
/
(
g_ewald
*
g_ewald
*
volume
);
eatom
[
i
]
*=
qscale
;
}
for
(
i
=
nlocal
;
i
<
ntotal
;
i
++
)
eatom
[
i
]
*=
0.5
*
qscale
;
}
if
(
vflag_atom
)
{
for
(
i
=
0
;
i
<
ntotal
;
i
++
)
for
(
j
=
0
;
j
<
6
;
j
++
)
vatom
[
i
][
j
]
*=
0.5
*
qscale
;
}
}
// 2d slab correction
if
(
slabflag
==
1
)
slabcorr
();
// convert atoms back from lamda to box coords
if
(
triclinic
)
domain
->
lamda2x
(
atom
->
nlocal
);
}
/* ----------------------------------------------------------------------
find center grid pt for each of my particles
check that full stencil for the particle will fit in my 3d brick
store central grid pt indices in part2grid array
------------------------------------------------------------------------- */
template
<
class
flt_t
,
class
acc_t
>
void
PPPMIntel
::
particle_map
(
IntelBuffers
<
flt_t
,
acc_t
>
*
buffers
)
{
int
nx
,
ny
,
nz
;
ATOM_T
*
_noalias
const
x
=
buffers
->
get_x
(
0
);
int
nlocal
=
atom
->
nlocal
;
int
flag
=
0
;
if
(
!
ISFINITE
(
boxlo
[
0
])
||
!
ISFINITE
(
boxlo
[
1
])
||
!
ISFINITE
(
boxlo
[
2
]))
error
->
one
(
FLERR
,
"Non-numeric box dimensions - simulation unstable"
);
const
flt_t
lo0
=
boxlo
[
0
];
const
flt_t
lo1
=
boxlo
[
1
];
const
flt_t
lo2
=
boxlo
[
2
];
const
flt_t
xi
=
delxinv
;
const
flt_t
yi
=
delyinv
;
const
flt_t
zi
=
delzinv
;
const
flt_t
fshift
=
shift
;
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned
#pragma simd
#endif
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// current particle coord can be outside global and local box
// add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1
nx
=
static_cast
<
int
>
((
x
[
i
].
x
-
lo0
)
*
xi
+
fshift
)
-
OFFSET
;
ny
=
static_cast
<
int
>
((
x
[
i
].
y
-
lo1
)
*
yi
+
fshift
)
-
OFFSET
;
nz
=
static_cast
<
int
>
((
x
[
i
].
z
-
lo2
)
*
zi
+
fshift
)
-
OFFSET
;
part2grid
[
i
][
0
]
=
nx
;
part2grid
[
i
][
1
]
=
ny
;
part2grid
[
i
][
2
]
=
nz
;
// check that entire stencil around nx,ny,nz will fit in my 3d brick
if
(
nx
+
nlower
<
nxlo_out
||
nx
+
nupper
>
nxhi_out
||
ny
+
nlower
<
nylo_out
||
ny
+
nupper
>
nyhi_out
||
nz
+
nlower
<
nzlo_out
||
nz
+
nupper
>
nzhi_out
)
flag
=
1
;
}
if
(
flag
)
error
->
one
(
FLERR
,
"Out of range atoms - cannot compute PPPM"
);
}
/* ----------------------------------------------------------------------
create discretized "density" on section of global grid due to my particles
density(x,y,z) = charge "density" at grid points of my 3d brick
(nxlo:nxhi,nylo:nyhi,nzlo:nzhi) is extent of my brick (including ghosts)
in global grid
------------------------------------------------------------------------- */
template
<
class
flt_t
,
class
acc_t
>
void
PPPMIntel
::
make_rho
(
IntelBuffers
<
flt_t
,
acc_t
>
*
buffers
)
{
// clear 3d density array
memset
(
&
(
density_brick
[
nzlo_out
][
nylo_out
][
nxlo_out
]),
0
,
ngrid
*
sizeof
(
FFT_SCALAR
));
// loop over my charges, add their contribution to nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
ATOM_T
*
_noalias
const
x
=
buffers
->
get_x
(
0
);
flt_t
*
_noalias
const
q
=
buffers
->
get_q
(
0
);
int
nlocal
=
atom
->
nlocal
;
const
flt_t
lo0
=
boxlo
[
0
];
const
flt_t
lo1
=
boxlo
[
1
];
const
flt_t
lo2
=
boxlo
[
2
];
const
flt_t
xi
=
delxinv
;
const
flt_t
yi
=
delyinv
;
const
flt_t
zi
=
delzinv
;
const
flt_t
fshift
=
shift
;
const
flt_t
fshiftone
=
shiftone
;
const
flt_t
fdelvolinv
=
delvolinv
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
int
nx
=
part2grid
[
i
][
0
];
int
ny
=
part2grid
[
i
][
1
];
int
nz
=
part2grid
[
i
][
2
];
FFT_SCALAR
dx
=
nx
+
fshiftone
-
(
x
[
i
].
x
-
lo0
)
*
xi
;
FFT_SCALAR
dy
=
ny
+
fshiftone
-
(
x
[
i
].
y
-
lo1
)
*
yi
;
FFT_SCALAR
dz
=
nz
+
fshiftone
-
(
x
[
i
].
z
-
lo2
)
*
zi
;
flt_t
rho
[
3
][
INTEL_P3M_MAXORDER
];
for
(
int
k
=
nlower
;
k
<=
nupper
;
k
++
)
{
FFT_SCALAR
r1
,
r2
,
r3
;
r1
=
r2
=
r3
=
ZEROF
;
for
(
int
l
=
order
-
1
;
l
>=
0
;
l
--
)
{
r1
=
rho_coeff
[
l
][
k
]
+
r1
*
dx
;
r2
=
rho_coeff
[
l
][
k
]
+
r2
*
dy
;
r3
=
rho_coeff
[
l
][
k
]
+
r3
*
dz
;
}
rho
[
0
][
k
-
nlower
]
=
r1
;
rho
[
1
][
k
-
nlower
]
=
r2
;
rho
[
2
][
k
-
nlower
]
=
r3
;
}
FFT_SCALAR
z0
=
fdelvolinv
*
q
[
i
];
for
(
int
n
=
nlower
;
n
<=
nupper
;
n
++
)
{
int
mz
=
n
+
nz
;
FFT_SCALAR
y0
=
z0
*
rho
[
2
][
n
-
nlower
];
for
(
int
m
=
nlower
;
m
<=
nupper
;
m
++
)
{
int
my
=
m
+
ny
;
FFT_SCALAR
x0
=
y0
*
rho
[
1
][
m
-
nlower
];
for
(
int
l
=
nlower
;
l
<=
nupper
;
l
++
)
{
int
mx
=
l
+
nx
;
density_brick
[
mz
][
my
][
mx
]
+=
x0
*
rho
[
0
][
l
-
nlower
];
}
}
}
}
}
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles for ik
------------------------------------------------------------------------- */
template
<
class
flt_t
,
class
acc_t
>
void
PPPMIntel
::
fieldforce_ik
(
IntelBuffers
<
flt_t
,
acc_t
>
*
buffers
)
{
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
ATOM_T
*
_noalias
const
x
=
buffers
->
get_x
(
0
);
flt_t
*
_noalias
const
q
=
buffers
->
get_q
(
0
);
FORCE_T
*
_noalias
const
f
=
buffers
->
get_f
();
int
nlocal
=
atom
->
nlocal
;
const
flt_t
lo0
=
boxlo
[
0
];
const
flt_t
lo1
=
boxlo
[
1
];
const
flt_t
lo2
=
boxlo
[
2
];
const
flt_t
xi
=
delxinv
;
const
flt_t
yi
=
delyinv
;
const
flt_t
zi
=
delzinv
;
const
flt_t
fshiftone
=
shiftone
;
const
flt_t
fqqrd2es
=
qqrd2e
*
scale
;
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned nontemporal
#pragma simd
#endif
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
int
nx
=
part2grid
[
i
][
0
];
int
ny
=
part2grid
[
i
][
1
];
int
nz
=
part2grid
[
i
][
2
];
FFT_SCALAR
dx
=
nx
+
fshiftone
-
(
x
[
i
].
x
-
lo0
)
*
xi
;
FFT_SCALAR
dy
=
ny
+
fshiftone
-
(
x
[
i
].
y
-
lo1
)
*
yi
;
FFT_SCALAR
dz
=
nz
+
fshiftone
-
(
x
[
i
].
z
-
lo2
)
*
zi
;
flt_t
rho
[
3
][
INTEL_P3M_MAXORDER
];
for
(
int
k
=
nlower
;
k
<=
nupper
;
k
++
)
{
FFT_SCALAR
r1
=
rho_coeff
[
order
-
1
][
k
];
FFT_SCALAR
r2
=
rho_coeff
[
order
-
1
][
k
];
FFT_SCALAR
r3
=
rho_coeff
[
order
-
1
][
k
];
for
(
int
l
=
order
-
2
;
l
>=
0
;
l
--
)
{
r1
=
rho_coeff
[
l
][
k
]
+
r1
*
dx
;
r2
=
rho_coeff
[
l
][
k
]
+
r2
*
dy
;
r3
=
rho_coeff
[
l
][
k
]
+
r3
*
dz
;
}
rho
[
0
][
k
-
nlower
]
=
r1
;
rho
[
1
][
k
-
nlower
]
=
r2
;
rho
[
2
][
k
-
nlower
]
=
r3
;
}
FFT_SCALAR
ekx
,
eky
,
ekz
;
ekx
=
eky
=
ekz
=
ZEROF
;
for
(
int
n
=
nlower
;
n
<=
nupper
;
n
++
)
{
int
mz
=
n
+
nz
;
FFT_SCALAR
z0
=
rho
[
2
][
n
-
nlower
];
for
(
int
m
=
nlower
;
m
<=
nupper
;
m
++
)
{
int
my
=
m
+
ny
;
FFT_SCALAR
y0
=
z0
*
rho
[
1
][
m
-
nlower
];
for
(
int
l
=
nlower
;
l
<=
nupper
;
l
++
)
{
int
mx
=
l
+
nx
;
FFT_SCALAR
x0
=
y0
*
rho
[
0
][
l
-
nlower
];
ekx
-=
x0
*
vdx_brick
[
mz
][
my
][
mx
];
eky
-=
x0
*
vdy_brick
[
mz
][
my
][
mx
];
ekz
-=
x0
*
vdz_brick
[
mz
][
my
][
mx
];
}
}
}
// convert E-field to force
const
flt_t
qfactor
=
fqqrd2es
*
q
[
i
];
f
[
i
].
x
+=
qfactor
*
ekx
;
f
[
i
].
y
+=
qfactor
*
eky
;
if
(
slabflag
!=
2
)
f
[
i
].
z
+=
qfactor
*
ekz
;
}
}
/* ----------------------------------------------------------------------
interpolate from grid to get electric field & force on my particles for ad
------------------------------------------------------------------------- */
template
<
class
flt_t
,
class
acc_t
>
void
PPPMIntel
::
fieldforce_ad
(
IntelBuffers
<
flt_t
,
acc_t
>
*
buffers
)
{
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
// (dx,dy,dz) = distance to "lower left" grid pt
// (mx,my,mz) = global coords of moving stencil pt
// ek = 3 components of E-field on particle
ATOM_T
*
_noalias
const
x
=
buffers
->
get_x
(
0
);
const
flt_t
*
_noalias
const
q
=
buffers
->
get_q
(
0
);
FORCE_T
*
_noalias
const
f
=
buffers
->
get_f
();
int
nlocal
=
atom
->
nlocal
;
const
flt_t
ftwo_pi
=
MY_PI
*
2.0
;
const
flt_t
ffour_pi
=
MY_PI
*
4.0
;
const
flt_t
lo0
=
boxlo
[
0
];
const
flt_t
lo1
=
boxlo
[
1
];
const
flt_t
lo2
=
boxlo
[
2
];
const
flt_t
xi
=
delxinv
;
const
flt_t
yi
=
delyinv
;
const
flt_t
zi
=
delzinv
;
const
flt_t
fshiftone
=
shiftone
;
const
flt_t
fqqrd2es
=
qqrd2e
*
scale
;
const
double
*
prd
=
domain
->
prd
;
const
double
xprd
=
prd
[
0
];
const
double
yprd
=
prd
[
1
];
const
double
zprd
=
prd
[
2
];
const
flt_t
hx_inv
=
nx_pppm
/
xprd
;
const
flt_t
hy_inv
=
ny_pppm
/
yprd
;
const
flt_t
hz_inv
=
nz_pppm
/
zprd
;
const
flt_t
fsf_coeff0
=
sf_coeff
[
0
];
const
flt_t
fsf_coeff1
=
sf_coeff
[
1
];
const
flt_t
fsf_coeff2
=
sf_coeff
[
2
];
const
flt_t
fsf_coeff3
=
sf_coeff
[
3
];
const
flt_t
fsf_coeff4
=
sf_coeff
[
4
];
const
flt_t
fsf_coeff5
=
sf_coeff
[
5
];
#if defined(LMP_SIMD_COMPILER)
#pragma vector aligned nontemporal
#pragma simd
#endif
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
int
nx
=
part2grid
[
i
][
0
];
int
ny
=
part2grid
[
i
][
1
];
int
nz
=
part2grid
[
i
][
2
];
FFT_SCALAR
dx
=
nx
+
fshiftone
-
(
x
[
i
].
x
-
lo0
)
*
xi
;
FFT_SCALAR
dy
=
ny
+
fshiftone
-
(
x
[
i
].
y
-
lo1
)
*
yi
;
FFT_SCALAR
dz
=
nz
+
fshiftone
-
(
x
[
i
].
z
-
lo2
)
*
zi
;
flt_t
rho
[
3
][
INTEL_P3M_MAXORDER
];
flt_t
drho
[
3
][
INTEL_P3M_MAXORDER
];
for
(
int
k
=
nlower
;
k
<=
nupper
;
k
++
)
{
FFT_SCALAR
r1
,
r2
,
r3
,
dr1
,
dr2
,
dr3
;
dr1
=
dr2
=
dr3
=
ZEROF
;
r1
=
rho_coeff
[
order
-
1
][
k
];
r2
=
rho_coeff
[
order
-
1
][
k
];
r3
=
rho_coeff
[
order
-
1
][
k
];
for
(
int
l
=
order
-
2
;
l
>=
0
;
l
--
)
{
r1
=
rho_coeff
[
l
][
k
]
+
r1
*
dx
;
r2
=
rho_coeff
[
l
][
k
]
+
r2
*
dy
;
r3
=
rho_coeff
[
l
][
k
]
+
r3
*
dz
;
dr1
=
drho_coeff
[
l
][
k
]
+
dr1
*
dx
;
dr2
=
drho_coeff
[
l
][
k
]
+
dr2
*
dy
;
dr3
=
drho_coeff
[
l
][
k
]
+
dr3
*
dz
;
}
rho
[
0
][
k
-
nlower
]
=
r1
;
rho
[
1
][
k
-
nlower
]
=
r2
;
rho
[
2
][
k
-
nlower
]
=
r3
;
drho
[
0
][
k
-
nlower
]
=
dr1
;
drho
[
1
][
k
-
nlower
]
=
dr2
;
drho
[
2
][
k
-
nlower
]
=
dr3
;
}
FFT_SCALAR
ekx
,
eky
,
ekz
;
ekx
=
eky
=
ekz
=
ZEROF
;
for
(
int
n
=
nlower
;
n
<=
nupper
;
n
++
)
{
int
mz
=
n
+
nz
;
for
(
int
m
=
nlower
;
m
<=
nupper
;
m
++
)
{
int
my
=
m
+
ny
;
FFT_SCALAR
ekx_p
=
rho
[
1
][
m
-
nlower
]
*
rho
[
2
][
n
-
nlower
];
FFT_SCALAR
eky_p
=
drho
[
1
][
m
-
nlower
]
*
rho
[
2
][
n
-
nlower
];
FFT_SCALAR
ekz_p
=
rho
[
1
][
m
-
nlower
]
*
drho
[
2
][
n
-
nlower
];
for
(
int
l
=
nlower
;
l
<=
nupper
;
l
++
)
{
int
mx
=
l
+
nx
;
ekx
+=
drho
[
0
][
l
-
nlower
]
*
ekx_p
*
u_brick
[
mz
][
my
][
mx
];
eky
+=
rho
[
0
][
l
-
nlower
]
*
eky_p
*
u_brick
[
mz
][
my
][
mx
];
ekz
+=
rho
[
0
][
l
-
nlower
]
*
ekz_p
*
u_brick
[
mz
][
my
][
mx
];
}
}
}
ekx
*=
hx_inv
;
eky
*=
hy_inv
;
ekz
*=
hz_inv
;
// convert E-field to force
const
flt_t
qfactor
=
fqqrd2es
*
q
[
i
];
const
flt_t
twoqsq
=
(
flt_t
)
2.0
*
q
[
i
]
*
q
[
i
];
const
flt_t
s1
=
x
[
i
].
x
*
hx_inv
;
const
flt_t
s2
=
x
[
i
].
y
*
hy_inv
;
const
flt_t
s3
=
x
[
i
].
z
*
hz_inv
;
flt_t
sf
=
fsf_coeff0
*
sin
(
ftwo_pi
*
s1
);
sf
+=
fsf_coeff1
*
sin
(
ffour_pi
*
s1
);
sf
*=
twoqsq
;
f
[
i
].
x
+=
qfactor
*
ekx
-
fqqrd2es
*
sf
;
sf
=
fsf_coeff2
*
sin
(
ftwo_pi
*
s2
);
sf
+=
fsf_coeff3
*
sin
(
ffour_pi
*
s2
);
sf
*=
twoqsq
;
f
[
i
].
y
+=
qfactor
*
eky
-
fqqrd2es
*
sf
;
sf
=
fsf_coeff4
*
sin
(
ftwo_pi
*
s3
);
sf
+=
fsf_coeff5
*
sin
(
ffour_pi
*
s3
);
sf
*=
twoqsq
;
if
(
slabflag
!=
2
)
f
[
i
].
z
+=
qfactor
*
ekz
-
fqqrd2es
*
sf
;
}
}
/* ----------------------------------------------------------------------
Pack data into intel package buffers if using LRT mode
------------------------------------------------------------------------- */
void
PPPMIntel
::
pack_buffers
()
{
fix
->
start_watch
(
TIME_PACK
);
#if defined(_OPENMP)
#pragma omp parallel default(none)
#endif
{
int
ifrom
,
ito
,
tid
;
IP_PRE_omp_range_id_align
(
ifrom
,
ito
,
tid
,
atom
->
nlocal
+
atom
->
nghost
,
comm
->
nthreads
,
sizeof
(
IntelBuffers
<
float
,
double
>::
atom_t
));
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_MIXED
)
fix
->
get_mixed_buffers
()
->
thr_pack
(
ifrom
,
ito
,
1
);
else
if
(
fix
->
precision
()
==
FixIntel
::
PREC_MODE_DOUBLE
)
fix
->
get_double_buffers
()
->
thr_pack
(
ifrom
,
ito
,
1
);
else
fix
->
get_single_buffers
()
->
thr_pack
(
ifrom
,
ito
,
1
);
}
fix
->
stop_watch
(
TIME_PACK
);
}
/* ----------------------------------------------------------------------
Returns 0 if Intel optimizations for PPPM ignored due to offload
------------------------------------------------------------------------- */
#ifdef _LMP_INTEL_OFFLOAD
int
PPPMIntel
::
use_base
()
{
return
_use_base
;
}
#endif
Event Timeline
Log In to Comment