Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88546572
ewald_n.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
Sat, Oct 19, 10:30
Size
22 KB
Mime Type
text/x-c
Expires
Mon, Oct 21, 10:30 (2 d)
Engine
blob
Format
Raw Data
Handle
21788826
Attached To
rLAMMPS lammps
ewald_n.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: Pieter J. in 't Veld (SNL)
------------------------------------------------------------------------- */
#include "mpi.h"
#include "string.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "ewald_n.h"
#include "math_vector.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
using
namespace
LAMMPS_NS
;
#define KSPACE_ILLEGAL "Illegal kspace_style ewald/n command"
#define KSPACE_ORDER "Unsupported order in kspace_style ewald/n for"
#define KSPACE_MIX "Unsupported mixing rule in kspace_style ewald/n for"
enum
{
GEOMETRIC
,
ARITHMETIC
,
SIXTHPOWER
};
// same as in pair.cpp
//#define DEBUG
/* ---------------------------------------------------------------------- */
EwaldN
::
EwaldN
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
KSpace
(
lmp
,
narg
,
arg
)
{
if
(
narg
!=
1
)
error
->
all
(
FLERR
,
KSPACE_ILLEGAL
);
precision
=
fabs
(
atof
(
arg
[
0
]));
memset
(
function
,
0
,
EWALD_NORDER
*
sizeof
(
int
));
kenergy
=
kvirial
=
NULL
;
cek_local
=
cek_global
=
NULL
;
ekr_local
=
NULL
;
hvec
=
NULL
;
kvec
=
NULL
;
B
=
NULL
;
first_output
=
0
;
}
EwaldN
::~
EwaldN
()
{
deallocate
();
delete
[]
ekr_local
;
delete
[]
B
;
}
/* --------------------------------------------------------------------- */
void
EwaldN
::
init
()
{
nkvec
=
nkvec_max
=
nevec
=
nevec_max
=
0
;
nfunctions
=
nsums
=
sums
=
0
;
nbox
=
-
1
;
bytes
=
0.0
;
if
(
!
comm
->
me
)
{
// output message
if
(
screen
)
fprintf
(
screen
,
"EwaldN initialization ...
\n
"
);
if
(
logfile
)
fprintf
(
logfile
,
"EwaldN initialization ...
\n
"
);
}
if
(
domain
->
dimension
==
2
)
// check for errors
error
->
all
(
FLERR
,
"Cannot use EwaldN with 2d simulation"
);
if
(
slabflag
==
0
&&
domain
->
nonperiodic
>
0
)
error
->
all
(
FLERR
,
"Cannot use nonperiodic boundaries with EwaldN"
);
if
(
slabflag
==
1
)
{
if
(
domain
->
xperiodic
!=
1
||
domain
->
yperiodic
!=
1
||
domain
->
boundary
[
2
][
0
]
!=
1
||
domain
->
boundary
[
2
][
1
]
!=
1
)
error
->
all
(
FLERR
,
"Incorrect boundaries with slab EwaldN"
);
}
scale
=
1.0
;
//mumurd2e = force->mumurd2e;
//dielectric = force->dielectric;
mumurd2e
=
dielectric
=
1.0
;
int
tmp
;
Pair
*
pair
=
force
->
pair
;
int
*
ptr
=
pair
?
(
int
*
)
pair
->
extract
(
"ewald_order"
,
tmp
)
:
NULL
;
double
*
cutoff
=
pair
?
(
double
*
)
pair
->
extract
(
"cut_coul"
,
tmp
)
:
NULL
;
if
(
!
(
ptr
||
cutoff
))
error
->
all
(
FLERR
,
"KSpace style is incompatible with Pair style"
);
int
ewald_order
=
ptr
?
*
((
int
*
)
ptr
)
:
1
<<
1
;
int
ewald_mix
=
ptr
?
*
((
int
*
)
pair
->
extract
(
"ewald_mix"
,
tmp
))
:
GEOMETRIC
;
memset
(
function
,
0
,
EWALD_NFUNCS
*
sizeof
(
int
));
for
(
int
i
=
0
;
i
<=
EWALD_NORDER
;
++
i
)
// transcribe order
if
(
ewald_order
&
(
1
<<
i
))
{
// from pair_style
int
n
[]
=
EWALD_NSUMS
,
k
;
char
str
[
128
];
switch
(
i
)
{
case
1
:
k
=
0
;
break
;
case
3
:
k
=
3
;
break
;
case
6
:
if
(
ewald_mix
==
GEOMETRIC
)
{
k
=
1
;
break
;
}
else
if
(
ewald_mix
==
ARITHMETIC
)
{
k
=
2
;
break
;
}
sprintf
(
str
,
"%s pair_style %s"
,
KSPACE_MIX
,
force
->
pair_style
);
error
->
all
(
FLERR
,
str
);
default
:
sprintf
(
str
,
"%s pair_style %s"
,
KSPACE_ORDER
,
force
->
pair_style
);
error
->
all
(
FLERR
,
str
);
}
nfunctions
+=
function
[
k
]
=
1
;
nsums
+=
n
[
k
];
}
g_ewald
=
(
1.35
-
0.15
*
log
(
precision
))
/
*
cutoff
;
// determine resolution
g2_max
=
-
4.0
*
g_ewald
*
g_ewald
*
log
(
precision
);
if
(
!
comm
->
me
)
{
// output results
if
(
screen
)
fprintf
(
screen
,
" G vector = %g
\n
"
,
g_ewald
);
if
(
logfile
)
fprintf
(
logfile
,
" G vector = %g
\n
"
,
g_ewald
);
}
}
/* ----------------------------------------------------------------------
adjust EwaldN coeffs, called initially and whenever volume has changed
------------------------------------------------------------------------- */
void
EwaldN
::
setup
()
{
volume
=
shape_det
(
domain
->
h
)
*
slab_volfactor
;
// cell volume
memcpy
(
unit
,
domain
->
h_inv
,
sizeof
(
shape
));
// wave vector units
shape_scalar_mult
(
unit
,
2.0
*
M_PI
);
unit
[
2
]
/=
slab_volfactor
;
//int nbox_old = nbox, nkvec_old = nkvec;
if
(
precision
>=
1
)
nbox
=
0
;
else
{
vector
n
=
{
1.0
,
1.0
,
1.0
};
// based on cutoff
vec_scalar_mult
(
n
,
g_ewald
*
sqrt
(
-
log
(
precision
))
/
M_PI
);
shape_vec_dot
(
n
,
n
,
domain
->
h
);
n
[
2
]
*=
slab_volfactor
;
nbox
=
(
int
)
n
[
0
];
if
(
nbox
<
(
int
)
n
[
1
])
nbox
=
(
int
)
n
[
1
];
if
(
nbox
<
(
int
)
n
[
2
])
nbox
=
(
int
)
n
[
2
];
}
reallocate
();
coefficients
();
// compute coeffs
init_coeffs
();
init_coeff_sums
();
init_self
();
if
(
!
(
first_output
||
comm
->
me
))
{
// output on first
first_output
=
1
;
if
(
screen
)
fprintf
(
screen
,
" vectors: nbox = %d, nkvec = %d
\n
"
,
nbox
,
nkvec
);
if
(
logfile
)
fprintf
(
logfile
,
" vectors: nbox = %d, nkvec = %d
\n
"
,
nbox
,
nkvec
);
}
}
void
EwaldN
::
reallocate
()
// allocate memory
{
int
ix
,
iy
,
iz
;
int
nkvec_max
=
nkvec
;
vector
h
;
nkvec
=
0
;
// determine size(kvec)
int
kflag
[(
nbox
+
1
)
*
(
2
*
nbox
+
1
)
*
(
2
*
nbox
+
1
)],
*
flag
=
kflag
;
for
(
ix
=
0
;
ix
<=
nbox
;
++
ix
)
for
(
iy
=-
nbox
;
iy
<=
nbox
;
++
iy
)
for
(
iz
=-
nbox
;
iz
<=
nbox
;
++
iz
)
if
(
!
(
ix
||
iy
||
iz
))
*
(
flag
++
)
=
0
;
else
if
((
!
ix
)
&&
(
iy
<
0
))
*
(
flag
++
)
=
0
;
else
if
((
!
(
ix
||
iy
))
&&
(
iz
<
0
))
*
(
flag
++
)
=
0
;
// use symmetry
else
{
h
[
0
]
=
unit
[
0
]
*
ix
;
h
[
1
]
=
unit
[
5
]
*
ix
+
unit
[
1
]
*
iy
;
h
[
2
]
=
unit
[
4
]
*
ix
+
unit
[
3
]
*
iy
+
unit
[
2
]
*
iz
;
if
((
*
(
flag
++
)
=
h
[
0
]
*
h
[
0
]
+
h
[
1
]
*
h
[
1
]
+
h
[
2
]
*
h
[
2
]
<=
g2_max
))
++
nkvec
;
}
if
(
nkvec
>
nkvec_max
)
{
deallocate
();
// free memory
hvec
=
new
hvector
[
nkvec
];
// hvec
bytes
+=
(
nkvec
-
nkvec_max
)
*
sizeof
(
hvector
);
kvec
=
new
kvector
[
nkvec
];
// kvec
bytes
+=
(
nkvec
-
nkvec_max
)
*
sizeof
(
kvector
);
kenergy
=
new
double
[
nkvec
*
nfunctions
];
// kenergy
bytes
+=
(
nkvec
-
nkvec_max
)
*
nfunctions
*
sizeof
(
double
);
kvirial
=
new
double
[
6
*
nkvec
*
nfunctions
];
// kvirial
bytes
+=
6
*
(
nkvec
-
nkvec_max
)
*
nfunctions
*
sizeof
(
double
);
cek_local
=
new
complex
[
nkvec
*
nsums
];
// cek_local
bytes
+=
(
nkvec
-
nkvec_max
)
*
nsums
*
sizeof
(
complex
);
cek_global
=
new
complex
[
nkvec
*
nsums
];
// cek_global
bytes
+=
(
nkvec
-
nkvec_max
)
*
nsums
*
sizeof
(
complex
);
nkvec_max
=
nkvec
;
}
flag
=
kflag
;
// create index and
kvector
*
k
=
kvec
;
// wave vectors
hvector
*
hi
=
hvec
;
for
(
ix
=
0
;
ix
<=
nbox
;
++
ix
)
for
(
iy
=-
nbox
;
iy
<=
nbox
;
++
iy
)
for
(
iz
=-
nbox
;
iz
<=
nbox
;
++
iz
)
if
(
*
(
flag
++
))
{
hi
->
x
=
unit
[
0
]
*
ix
;
hi
->
y
=
unit
[
5
]
*
ix
+
unit
[
1
]
*
iy
;
(
hi
++
)
->
z
=
unit
[
4
]
*
ix
+
unit
[
3
]
*
iy
+
unit
[
2
]
*
iz
;
k
->
x
=
ix
+
nbox
;
k
->
y
=
iy
+
nbox
;
(
k
++
)
->
z
=
iz
+
nbox
;
}
}
void
EwaldN
::
reallocate_atoms
()
{
if
((
nevec
=
atom
->
nmax
*
(
2
*
nbox
+
1
))
<=
nevec_max
)
return
;
delete
[]
ekr_local
;
ekr_local
=
new
cvector
[
nevec
];
bytes
+=
(
nevec
-
nevec_max
)
*
sizeof
(
cvector
);
nevec_max
=
nevec
;
}
void
EwaldN
::
deallocate
()
// free memory
{
delete
[]
hvec
;
hvec
=
NULL
;
delete
[]
kvec
;
kvec
=
NULL
;
delete
[]
kenergy
;
kenergy
=
NULL
;
delete
[]
kvirial
;
kvirial
=
NULL
;
delete
[]
cek_local
;
cek_local
=
NULL
;
delete
[]
cek_global
;
cek_global
=
NULL
;
}
void
EwaldN
::
coefficients
()
// set up pre-factors
{
vector
h
;
hvector
*
hi
=
hvec
,
*
nh
;
double
eta2
=
0.25
/
(
g_ewald
*
g_ewald
);
double
b1
,
b2
,
expb2
,
h1
,
h2
,
c1
,
c2
;
double
*
ke
=
kenergy
,
*
kv
=
kvirial
;
int
func0
=
function
[
0
],
func12
=
function
[
1
]
||
function
[
2
],
func3
=
function
[
3
];
for
(
nh
=
(
hi
=
hvec
)
+
nkvec
;
hi
<
nh
;
++
hi
)
{
// wave vectors
memcpy
(
h
,
hi
,
sizeof
(
vector
));
expb2
=
exp
(
-
(
b2
=
(
h2
=
vec_dot
(
h
,
h
))
*
eta2
));
if
(
func0
)
{
// qi*qj/r coeffs
*
(
ke
++
)
=
c1
=
expb2
/
h2
;
*
(
kv
++
)
=
c1
-
(
c2
=
2.0
*
c1
*
(
1.0
+
b2
)
/
h2
)
*
h
[
0
]
*
h
[
0
];
*
(
kv
++
)
=
c1
-
c2
*
h
[
1
]
*
h
[
1
];
// lammps convention
*
(
kv
++
)
=
c1
-
c2
*
h
[
2
]
*
h
[
2
];
// instead of voigt
*
(
kv
++
)
=
-
c2
*
h
[
1
]
*
h
[
0
];
*
(
kv
++
)
=
-
c2
*
h
[
2
]
*
h
[
0
];
*
(
kv
++
)
=
-
c2
*
h
[
2
]
*
h
[
1
];
}
if
(
func12
)
{
// -Bij/r^6 coeffs
b1
=
sqrt
(
b2
);
// minus sign folded
h1
=
sqrt
(
h2
);
// into constants
*
(
ke
++
)
=
c1
=
-
h1
*
h2
*
((
c2
=
sqrt
(
M_PI
)
*
erfc
(
b1
))
+
(
0.5
/
b2
-
1.0
)
*
expb2
/
b1
);
*
(
kv
++
)
=
c1
-
(
c2
=
3.0
*
h1
*
(
c2
-
expb2
/
b1
))
*
h
[
0
]
*
h
[
0
];
*
(
kv
++
)
=
c1
-
c2
*
h
[
1
]
*
h
[
1
];
// lammps convention
*
(
kv
++
)
=
c1
-
c2
*
h
[
2
]
*
h
[
2
];
// instead of voigt
*
(
kv
++
)
=
-
c2
*
h
[
1
]
*
h
[
0
];
*
(
kv
++
)
=
-
c2
*
h
[
2
]
*
h
[
0
];
*
(
kv
++
)
=
-
c2
*
h
[
2
]
*
h
[
1
];
}
if
(
func3
)
{
// dipole coeffs
*
(
ke
++
)
=
c1
=
expb2
/
h2
;
*
(
kv
++
)
=
c1
-
(
c2
=
2.0
*
c1
*
(
1.0
+
b2
)
/
h2
)
*
h
[
0
]
*
h
[
0
];
*
(
kv
++
)
=
c1
-
c2
*
h
[
1
]
*
h
[
1
];
// lammps convention
*
(
kv
++
)
=
c1
-
c2
*
h
[
2
]
*
h
[
2
];
// instead of voigt
*
(
kv
++
)
=
-
c2
*
h
[
1
]
*
h
[
0
];
*
(
kv
++
)
=
-
c2
*
h
[
2
]
*
h
[
0
];
*
(
kv
++
)
=
-
c2
*
h
[
2
]
*
h
[
1
];
}
}
}
void
EwaldN
::
init_coeffs
()
// local pair coeffs
{
int
tmp
;
int
n
=
atom
->
ntypes
;
if
(
function
[
1
])
{
// geometric 1/r^6
double
**
b
=
(
double
**
)
force
->
pair
->
extract
(
"B"
,
tmp
);
delete
[]
B
;
B
=
new
double
[
n
+
1
];
bytes
+=
(
n
+
1
)
*
sizeof
(
double
);
for
(
int
i
=
0
;
i
<=
n
;
++
i
)
B
[
i
]
=
sqrt
(
fabs
(
b
[
i
][
i
]));
}
if
(
function
[
2
])
{
// arithmetic 1/r^6
double
**
epsilon
=
(
double
**
)
force
->
pair
->
extract
(
"epsilon"
,
tmp
);
double
**
sigma
=
(
double
**
)
force
->
pair
->
extract
(
"sigma"
,
tmp
);
if
(
!
(
epsilon
&&
sigma
))
error
->
all
(
FLERR
,
"epsilon or sigma reference not set by pair style in ewald/n"
);
double
eps_i
,
sigma_i
,
sigma_n
,
*
bi
=
B
=
new
double
[
7
*
n
+
7
];
double
c
[
7
]
=
{
1.0
,
sqrt
(
6.0
),
sqrt
(
15.0
),
sqrt
(
20.0
),
sqrt
(
15.0
),
sqrt
(
6.0
),
1.0
};
for
(
int
i
=
0
;
i
<=
n
;
++
i
)
{
eps_i
=
sqrt
(
epsilon
[
i
][
i
]);
sigma_i
=
sigma
[
i
][
i
];
sigma_n
=
1.0
;
for
(
int
j
=
0
;
j
<
7
;
++
j
)
{
*
(
bi
++
)
=
sigma_n
*
eps_i
*
c
[
j
];
sigma_n
*=
sigma_i
;
}
}
}
}
void
EwaldN
::
init_coeff_sums
()
// sums based on atoms
{
if
(
sums
)
return
;
// calculated only once
sums
=
1
;
Sum
sum_local
[
EWALD_MAX_NSUMS
];
memset
(
sum_local
,
0
,
EWALD_MAX_NSUMS
*
sizeof
(
Sum
));
if
(
function
[
0
])
{
// 1/r
double
*
q
=
atom
->
q
,
*
qn
=
q
+
atom
->
nlocal
;
for
(
double
*
i
=
q
;
i
<
qn
;
++
i
)
{
sum_local
[
0
].
x
+=
i
[
0
];
sum_local
[
0
].
x2
+=
i
[
0
]
*
i
[
0
];
}
}
if
(
function
[
1
])
{
// geometric 1/r^6
int
*
type
=
atom
->
type
,
*
ntype
=
type
+
atom
->
nlocal
;
for
(
int
*
i
=
type
;
i
<
ntype
;
++
i
)
{
sum_local
[
1
].
x
+=
B
[
i
[
0
]];
sum_local
[
1
].
x2
+=
B
[
i
[
0
]]
*
B
[
i
[
0
]];
}
}
if
(
function
[
2
])
{
// aritmetic 1/r^6
double
*
bi
;
int
*
type
=
atom
->
type
,
*
ntype
=
type
+
atom
->
nlocal
;
for
(
int
*
i
=
type
;
i
<
ntype
;
++
i
)
{
bi
=
B
+
7
*
i
[
0
];
sum_local
[
2
].
x2
+=
bi
[
0
]
*
bi
[
6
];
for
(
int
k
=
2
;
k
<
9
;
++
k
)
sum_local
[
k
].
x
+=
*
(
bi
++
);
}
}
if
(
function
[
3
])
{
// dipole
double
**
mu
=
atom
->
mu
;
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
sum_local
[
9
].
x2
+=
mu
[
i
][
3
]
*
mu
[
i
][
3
];
}
MPI_Allreduce
(
sum_local
,
sum
,
2
*
EWALD_MAX_NSUMS
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
}
void
EwaldN
::
init_self
()
{
double
g1
=
g_ewald
,
g2
=
g1
*
g1
,
g3
=
g1
*
g2
;
memset
(
energy_self
,
0
,
EWALD_NFUNCS
*
sizeof
(
double
));
// self energy
memset
(
virial_self
,
0
,
EWALD_NFUNCS
*
sizeof
(
double
));
const
double
qscale
=
force
->
qqrd2e
*
scale
;
if
(
function
[
0
])
{
// 1/r
virial_self
[
0
]
=
-
0.5
*
M_PI
*
qscale
/
(
g2
*
volume
)
*
sum
[
0
].
x
*
sum
[
0
].
x
;
energy_self
[
0
]
=
sum
[
0
].
x2
*
qscale
*
g1
/
sqrt
(
M_PI
)
-
virial_self
[
0
];
}
if
(
function
[
1
])
{
// geometric 1/r^6
virial_self
[
1
]
=
M_PI
*
sqrt
(
M_PI
)
*
g3
/
(
6.0
*
volume
)
*
sum
[
1
].
x
*
sum
[
1
].
x
;
energy_self
[
1
]
=
-
sum
[
1
].
x2
*
g3
*
g3
/
12.0
+
virial_self
[
1
];
}
if
(
function
[
2
])
{
// arithmetic 1/r^6
virial_self
[
2
]
=
M_PI
*
sqrt
(
M_PI
)
*
g3
/
(
48.0
*
volume
)
*
(
sum
[
2
].
x
*
sum
[
8
].
x
+
sum
[
3
].
x
*
sum
[
7
].
x
+
sum
[
4
].
x
*
sum
[
6
].
x
+
0.5
*
sum
[
5
].
x
*
sum
[
5
].
x
);
energy_self
[
2
]
=
-
sum
[
2
].
x2
*
g3
*
g3
/
3.0
+
virial_self
[
2
];
}
if
(
function
[
3
])
{
// dipole
virial_self
[
3
]
=
0
;
// in surface
energy_self
[
3
]
=
sum
[
9
].
x2
*
mumurd2e
*
2.0
*
g3
/
3.0
/
sqrt
(
M_PI
)
-
virial_self
[
3
];
}
}
/* ----------------------------------------------------------------------
compute the EwaldN long-range force, energy, virial
------------------------------------------------------------------------- */
void
EwaldN
::
compute
(
int
eflag
,
int
vflag
)
{
if
(
!
nbox
)
return
;
reallocate_atoms
();
compute_ek
();
compute_force
();
compute_surface
();
compute_energy
(
eflag
);
compute_virial
(
vflag
);
}
void
EwaldN
::
compute_ek
()
{
cvector
*
ekr
=
ekr_local
;
int
lbytes
=
(
2
*
nbox
+
1
)
*
sizeof
(
cvector
);
hvector
*
h
;
kvector
*
k
,
*
nk
=
kvec
+
nkvec
;
cvector
z1
,
z
[
2
*
nbox
+
1
],
*
zx
,
*
zy
,
*
zz
,
*
zn
=
z
+
2
*
nbox
;
complex
*
cek
,
zxyz
,
zxy
,
cx
;
vector
mui
;
double
*
x
=
atom
->
x
[
0
],
*
xn
=
x
+
3
*
atom
->
nlocal
,
*
q
=
atom
->
q
,
qi
,
bi
,
ci
[
7
];
double
*
mu
=
atom
->
mu
?
atom
->
mu
[
0
]
:
NULL
;
int
i
,
kx
,
ky
,
n
=
nkvec
*
nsums
,
*
type
=
atom
->
type
,
tri
=
domain
->
triclinic
;
int
func
[
EWALD_NFUNCS
];
memcpy
(
func
,
function
,
EWALD_NFUNCS
*
sizeof
(
int
));
memset
(
cek_local
,
0
,
n
*
sizeof
(
complex
));
// reset sums
while
(
x
<
xn
)
{
zx
=
(
zy
=
(
zz
=
z
+
nbox
)
+
1
)
-
2
;
C_SET
(
zz
->
x
,
1
,
0
);
C_SET
(
zz
->
y
,
1
,
0
);
C_SET
(
zz
->
z
,
1
,
0
);
// z[0]
if
(
tri
)
{
// triclinic z[1]
C_ANGLE
(
z1
.
x
,
unit
[
0
]
*
x
[
0
]
+
unit
[
5
]
*
x
[
1
]
+
unit
[
4
]
*
x
[
2
]);
C_ANGLE
(
z1
.
y
,
unit
[
1
]
*
x
[
1
]
+
unit
[
3
]
*
x
[
2
]);
C_ANGLE
(
z1
.
z
,
x
[
2
]
*
unit
[
2
]);
x
+=
3
;
}
else
{
// orthogonal z[1]
C_ANGLE
(
z1
.
x
,
*
(
x
++
)
*
unit
[
0
]);
C_ANGLE
(
z1
.
y
,
*
(
x
++
)
*
unit
[
1
]);
C_ANGLE
(
z1
.
z
,
*
(
x
++
)
*
unit
[
2
]);
}
for
(;
zz
<
zn
;
--
zx
,
++
zy
,
++
zz
)
{
// set up z[k]=e^(ik.r)
C_RMULT
(
zy
->
x
,
zz
->
x
,
z1
.
x
);
// 3D k-vector
C_RMULT
(
zy
->
y
,
zz
->
y
,
z1
.
y
);
C_CONJ
(
zx
->
y
,
zy
->
y
);
C_RMULT
(
zy
->
z
,
zz
->
z
,
z1
.
z
);
C_CONJ
(
zx
->
z
,
zy
->
z
);
}
kx
=
ky
=
-
1
;
cek
=
cek_local
;
if
(
func
[
0
])
qi
=
*
(
q
++
);
if
(
func
[
1
])
bi
=
B
[
*
type
];
if
(
func
[
2
])
memcpy
(
ci
,
B
+
7
*
type
[
0
],
7
*
sizeof
(
double
));
if
(
func
[
3
])
{
memcpy
(
mui
,
mu
,
sizeof
(
vector
));
mu
+=
3
;
vec_scalar_mult
(
mui
,
mu
[
3
]);
mu
+=
4
;
h
=
hvec
;
}
for
(
k
=
kvec
;
k
<
nk
;
++
k
)
{
// compute rho(k)
if
(
ky
!=
k
->
y
)
{
// based on order in
if
(
kx
!=
k
->
x
)
cx
=
z
[
kx
=
k
->
x
].
x
;
// reallocate
C_RMULT
(
zxy
,
z
[
ky
=
k
->
y
].
y
,
cx
);
}
C_RMULT
(
zxyz
,
z
[
k
->
z
].
z
,
zxy
);
if
(
func
[
0
])
{
cek
->
re
+=
zxyz
.
re
*
qi
;
(
cek
++
)
->
im
+=
zxyz
.
im
*
qi
;
}
if
(
func
[
1
])
{
cek
->
re
+=
zxyz
.
re
*
bi
;
(
cek
++
)
->
im
+=
zxyz
.
im
*
bi
;
}
if
(
func
[
2
])
for
(
i
=
0
;
i
<
7
;
++
i
)
{
cek
->
re
+=
zxyz
.
re
*
ci
[
i
];
(
cek
++
)
->
im
+=
zxyz
.
im
*
ci
[
i
];
}
if
(
func
[
3
])
{
register
double
muk
=
mui
[
0
]
*
h
->
x
+
mui
[
1
]
*
h
->
y
+
mui
[
2
]
*
h
->
z
;
++
h
;
cek
->
re
+=
zxyz
.
re
*
muk
;
(
cek
++
)
->
im
+=
zxyz
.
im
*
muk
;
}
}
ekr
=
(
cvector
*
)
((
char
*
)
memcpy
(
ekr
,
z
,
lbytes
)
+
lbytes
);
++
type
;
}
MPI_Allreduce
(
cek_local
,
cek_global
,
2
*
n
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
}
void
EwaldN
::
compute_force
()
{
kvector
*
k
;
hvector
*
h
,
*
nh
;
cvector
*
z
=
ekr_local
;
vector
sum
[
EWALD_MAX_NSUMS
],
mui
;
complex
*
cek
,
zc
,
zx
,
zxy
;
double
*
f
=
atom
->
f
[
0
],
*
fn
=
f
+
3
*
atom
->
nlocal
,
*
q
=
atom
->
q
,
*
t
=
NULL
;
double
*
mu
=
atom
->
mu
?
atom
->
mu
[
0
]
:
NULL
;
const
double
qscale
=
force
->
qqrd2e
*
scale
;
double
*
ke
,
c
[
EWALD_NFUNCS
]
=
{
8.0
*
M_PI
*
qscale
/
volume
,
2.0
*
M_PI
*
sqrt
(
M_PI
)
/
(
12.0
*
volume
),
2.0
*
M_PI
*
sqrt
(
M_PI
)
/
(
192.0
*
volume
),
8.0
*
M_PI
*
mumurd2e
/
volume
};
double
kt
=
4.0
*
pow
(
g_ewald
,
3.0
)
/
3.0
/
sqrt
(
M_PI
)
/
c
[
3
];
int
i
,
kx
,
ky
,
lbytes
=
(
2
*
nbox
+
1
)
*
sizeof
(
cvector
),
*
type
=
atom
->
type
;
int
func
[
EWALD_NFUNCS
];
if
(
atom
->
torque
)
t
=
atom
->
torque
[
0
];
memcpy
(
func
,
function
,
EWALD_NFUNCS
*
sizeof
(
int
));
memset
(
sum
,
0
,
EWALD_MAX_NSUMS
*
sizeof
(
vector
));
// fj = -dE/dr =
for
(;
f
<
fn
;
f
+=
3
)
{
// -i*qj*fac*
k
=
kvec
;
// Sum[conj(d)-d]
kx
=
ky
=
-
1
;
// d = k*conj(ekj)*ek
ke
=
kenergy
;
cek
=
cek_global
;
memset
(
sum
,
0
,
EWALD_MAX_NSUMS
*
sizeof
(
vector
));
if
(
func
[
3
])
{
register
double
di
=
mu
[
3
]
*
c
[
3
];
mui
[
0
]
=
di
*
(
mu
++
)[
0
];
mui
[
1
]
=
di
*
(
mu
++
)[
1
];
mui
[
2
]
=
di
*
(
mu
++
)[
2
];
mu
++
;
}
for
(
nh
=
(
h
=
hvec
)
+
nkvec
;
h
<
nh
;
++
h
,
++
k
)
{
if
(
ky
!=
k
->
y
)
{
// based on order in
if
(
kx
!=
k
->
x
)
zx
=
z
[
kx
=
k
->
x
].
x
;
// reallocate
C_RMULT
(
zxy
,
z
[
ky
=
k
->
y
].
y
,
zx
);
}
C_CRMULT
(
zc
,
z
[
k
->
z
].
z
,
zxy
);
if
(
func
[
0
])
{
// 1/r
register
double
im
=
*
(
ke
++
)
*
(
zc
.
im
*
cek
->
re
+
cek
->
im
*
zc
.
re
);
++
cek
;
sum
[
0
][
0
]
+=
h
->
x
*
im
;
sum
[
0
][
1
]
+=
h
->
y
*
im
;
sum
[
0
][
2
]
+=
h
->
z
*
im
;
}
if
(
func
[
1
])
{
// geometric 1/r^6
register
double
im
=
*
(
ke
++
)
*
(
zc
.
im
*
cek
->
re
+
cek
->
im
*
zc
.
re
);
++
cek
;
sum
[
1
][
0
]
+=
h
->
x
*
im
;
sum
[
1
][
1
]
+=
h
->
y
*
im
;
sum
[
1
][
2
]
+=
h
->
z
*
im
;
}
if
(
func
[
2
])
{
// arithmetic 1/r^6
register
double
im
,
c
=
*
(
ke
++
);
for
(
i
=
2
;
i
<
9
;
++
i
)
{
im
=
c
*
(
zc
.
im
*
cek
->
re
+
cek
->
im
*
zc
.
re
);
++
cek
;
sum
[
i
][
0
]
+=
h
->
x
*
im
;
sum
[
i
][
1
]
+=
h
->
y
*
im
;
sum
[
i
][
2
]
+=
h
->
z
*
im
;
}
}
if
(
func
[
3
])
{
// dipole
register
double
im
=
*
(
ke
++
)
*
(
zc
.
im
*
cek
->
re
+
cek
->
im
*
zc
.
re
)
*
(
mui
[
0
]
*
h
->
x
+
mui
[
1
]
*
h
->
y
+
mui
[
2
]
*
h
->
z
);
++
cek
;
sum
[
9
][
0
]
+=
h
->
x
*
im
;
sum
[
9
][
1
]
+=
h
->
y
*
im
;
sum
[
9
][
2
]
+=
h
->
z
*
im
;
}
}
if
(
func
[
0
])
{
// 1/r
register
double
qi
=
*
(
q
++
)
*
c
[
0
];
f
[
0
]
-=
sum
[
0
][
0
]
*
qi
;
f
[
1
]
-=
sum
[
0
][
1
]
*
qi
;
f
[
2
]
-=
sum
[
0
][
2
]
*
qi
;
}
if
(
func
[
1
])
{
// geometric 1/r^6
register
double
bi
=
B
[
*
type
]
*
c
[
1
];
f
[
0
]
-=
sum
[
1
][
0
]
*
bi
;
f
[
1
]
-=
sum
[
1
][
1
]
*
bi
;
f
[
2
]
-=
sum
[
1
][
2
]
*
bi
;
}
if
(
func
[
2
])
{
// arithmetic 1/r^6
register
double
*
bi
=
B
+
7
*
type
[
0
]
+
7
;
for
(
i
=
2
;
i
<
9
;
++
i
)
{
register
double
c2
=
(
--
bi
)[
0
]
*
c
[
2
];
f
[
0
]
-=
sum
[
i
][
0
]
*
c2
;
f
[
1
]
-=
sum
[
i
][
1
]
*
c2
;
f
[
2
]
-=
sum
[
i
][
2
]
*
c2
;
}
}
if
(
func
[
3
])
{
// dipole
f
[
0
]
-=
sum
[
9
][
0
];
f
[
1
]
-=
sum
[
9
][
1
];
f
[
2
]
-=
sum
[
9
][
2
];
*
(
t
++
)
-=
mui
[
1
]
*
sum
[
0
][
2
]
+
mui
[
2
]
*
sum
[
0
][
1
]
-
mui
[
0
]
*
kt
;
// torque
*
(
t
++
)
-=
mui
[
2
]
*
sum
[
0
][
0
]
+
mui
[
0
]
*
sum
[
0
][
2
]
-
mui
[
1
]
*
kt
;
*
(
t
++
)
-=
mui
[
0
]
*
sum
[
0
][
1
]
+
mui
[
1
]
*
sum
[
0
][
0
]
-
mui
[
2
]
*
kt
;
}
z
=
(
cvector
*
)
((
char
*
)
z
+
lbytes
);
++
type
;
}
}
void
EwaldN
::
compute_surface
()
{
if
(
!
function
[
3
])
return
;
vector
sum_local
=
VECTOR_NULL
,
sum_total
;
double
**
mu
=
atom
->
mu
;
int
nlocal
=
atom
->
nlocal
;
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
register
double
di
=
mu
[
i
][
3
];
sum_local
[
0
]
+=
di
*
mu
[
i
][
0
];
sum_local
[
1
]
+=
di
*
mu
[
i
][
1
];
sum_local
[
2
]
+=
di
*
mu
[
i
][
2
];
}
MPI_Allreduce
(
sum_local
,
sum_total
,
3
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
energy_self
[
3
]
+=
virial_self
[
3
];
virial_self
[
3
]
=
mumurd2e
*
(
2.0
*
M_PI
*
vec_dot
(
sum_total
,
sum_total
)
/
(
2.0
*
dielectric
+
1
)
/
volume
);
energy_self
[
3
]
-=
virial_self
[
3
];
}
void
EwaldN
::
compute_energy
(
int
eflag
)
{
energy
=
0.0
;
if
(
!
eflag
)
return
;
complex
*
cek
=
cek_global
;
double
*
ke
=
kenergy
;
const
double
qscale
=
force
->
qqrd2e
*
scale
;
double
c
[
EWALD_NFUNCS
]
=
{
4.0
*
M_PI
*
qscale
/
volume
,
2.0
*
M_PI
*
sqrt
(
M_PI
)
/
(
24.0
*
volume
),
2.0
*
M_PI
*
sqrt
(
M_PI
)
/
(
192.0
*
volume
),
4.0
*
M_PI
*
mumurd2e
/
volume
};
double
sum
[
EWALD_NFUNCS
];
int
func
[
EWALD_NFUNCS
];
memcpy
(
func
,
function
,
EWALD_NFUNCS
*
sizeof
(
int
));
memset
(
sum
,
0
,
EWALD_NFUNCS
*
sizeof
(
double
));
// reset sums
for
(
int
k
=
0
;
k
<
nkvec
;
++
k
)
{
// sum over k vectors
if
(
func
[
0
])
{
// 1/r
sum
[
0
]
+=
*
(
ke
++
)
*
(
cek
->
re
*
cek
->
re
+
cek
->
im
*
cek
->
im
);
++
cek
;
}
if
(
func
[
1
])
{
// geometric 1/r^6
sum
[
1
]
+=
*
(
ke
++
)
*
(
cek
->
re
*
cek
->
re
+
cek
->
im
*
cek
->
im
);
++
cek
;
}
if
(
func
[
2
])
{
// arithmetic 1/r^6
register
double
r
=
(
cek
[
0
].
re
*
cek
[
6
].
re
+
cek
[
0
].
im
*
cek
[
6
].
im
)
+
(
cek
[
1
].
re
*
cek
[
5
].
re
+
cek
[
1
].
im
*
cek
[
5
].
im
)
+
(
cek
[
2
].
re
*
cek
[
4
].
re
+
cek
[
2
].
im
*
cek
[
4
].
im
)
+
0.5
*
(
cek
[
3
].
re
*
cek
[
3
].
re
+
cek
[
3
].
im
*
cek
[
3
].
im
);
cek
+=
7
;
sum
[
2
]
+=
*
(
ke
++
)
*
r
;
}
if
(
func
[
3
])
{
// dipole
sum
[
3
]
+=
*
(
ke
++
)
*
(
cek
->
re
*
cek
->
re
+
cek
->
im
*
cek
->
im
);
++
cek
;
}
}
for
(
int
k
=
0
;
k
<
EWALD_NFUNCS
;
++
k
)
energy
+=
c
[
k
]
*
sum
[
k
]
-
energy_self
[
k
];
if
(
slabflag
)
compute_slabcorr
(
eflag
);
}
#define swap(a, b) { register double t = a; a= b; b = t; }
void
EwaldN
::
compute_virial
(
int
vflag
)
{
memset
(
virial
,
0
,
sizeof
(
shape
));
if
(
!
vflag
)
return
;
complex
*
cek
=
cek_global
;
double
*
kv
=
kvirial
;
const
double
qscale
=
force
->
qqrd2e
*
scale
;
double
c
[
EWALD_NFUNCS
]
=
{
4.0
*
M_PI
*
qscale
/
volume
,
2.0
*
M_PI
*
sqrt
(
M_PI
)
/
(
24.0
*
volume
),
2.0
*
M_PI
*
sqrt
(
M_PI
)
/
(
192.0
*
volume
),
4.0
*
M_PI
*
mumurd2e
/
volume
};
shape
sum
[
EWALD_NFUNCS
];
int
func
[
EWALD_NFUNCS
];
memcpy
(
func
,
function
,
EWALD_NFUNCS
*
sizeof
(
int
));
memset
(
sum
,
0
,
EWALD_NFUNCS
*
sizeof
(
shape
));
for
(
int
k
=
0
;
k
<
nkvec
;
++
k
)
{
// sum over k vectors
if
(
func
[
0
])
{
// 1/r
register
double
r
=
cek
->
re
*
cek
->
re
+
cek
->
im
*
cek
->
im
;
++
cek
;
sum
[
0
][
0
]
+=
*
(
kv
++
)
*
r
;
sum
[
0
][
1
]
+=
*
(
kv
++
)
*
r
;
sum
[
0
][
2
]
+=
*
(
kv
++
)
*
r
;
sum
[
0
][
3
]
+=
*
(
kv
++
)
*
r
;
sum
[
0
][
4
]
+=
*
(
kv
++
)
*
r
;
sum
[
0
][
5
]
+=
*
(
kv
++
)
*
r
;
}
if
(
func
[
1
])
{
// geometric 1/r^6
register
double
r
=
cek
->
re
*
cek
->
re
+
cek
->
im
*
cek
->
im
;
++
cek
;
sum
[
1
][
0
]
+=
*
(
kv
++
)
*
r
;
sum
[
1
][
1
]
+=
*
(
kv
++
)
*
r
;
sum
[
1
][
2
]
+=
*
(
kv
++
)
*
r
;
sum
[
1
][
3
]
+=
*
(
kv
++
)
*
r
;
sum
[
1
][
4
]
+=
*
(
kv
++
)
*
r
;
sum
[
1
][
5
]
+=
*
(
kv
++
)
*
r
;
}
if
(
func
[
2
])
{
// arithmetic 1/r^6
register
double
r
=
(
cek
[
0
].
re
*
cek
[
6
].
re
+
cek
[
0
].
im
*
cek
[
6
].
im
)
+
(
cek
[
1
].
re
*
cek
[
5
].
re
+
cek
[
1
].
im
*
cek
[
5
].
im
)
+
(
cek
[
2
].
re
*
cek
[
4
].
re
+
cek
[
2
].
im
*
cek
[
4
].
im
)
+
0.5
*
(
cek
[
3
].
re
*
cek
[
3
].
re
+
cek
[
3
].
im
*
cek
[
3
].
im
);
cek
+=
7
;
sum
[
2
][
0
]
+=
*
(
kv
++
)
*
r
;
sum
[
2
][
1
]
+=
*
(
kv
++
)
*
r
;
sum
[
2
][
2
]
+=
*
(
kv
++
)
*
r
;
sum
[
2
][
3
]
+=
*
(
kv
++
)
*
r
;
sum
[
2
][
4
]
+=
*
(
kv
++
)
*
r
;
sum
[
2
][
5
]
+=
*
(
kv
++
)
*
r
;
}
if
(
func
[
3
])
{
register
double
r
=
cek
->
re
*
cek
->
re
+
cek
->
im
*
cek
->
im
;
++
cek
;
sum
[
3
][
0
]
+=
*
(
kv
++
)
*
r
;
sum
[
3
][
1
]
+=
*
(
kv
++
)
*
r
;
sum
[
3
][
2
]
+=
*
(
kv
++
)
*
r
;
sum
[
3
][
3
]
+=
*
(
kv
++
)
*
r
;
sum
[
3
][
4
]
+=
*
(
kv
++
)
*
r
;
sum
[
3
][
5
]
+=
*
(
kv
++
)
*
r
;
}
}
for
(
int
k
=
0
;
k
<
EWALD_NFUNCS
;
++
k
)
if
(
func
[
k
])
{
shape
self
=
{
virial_self
[
k
],
virial_self
[
k
],
virial_self
[
k
],
0
,
0
,
0
};
shape_scalar_mult
(
sum
[
k
],
c
[
k
]);
shape_add
(
virial
,
sum
[
k
]);
shape_subtr
(
virial
,
self
);
}
}
/* ----------------------------------------------------------------------
Slab-geometry correction term to dampen inter-slab interactions between
periodically repeating slabs. Yields good approximation to 2-D EwaldN if
adequate empty space is left between repeating slabs (J. Chem. Phys.
111, 3155). Slabs defined here to be parallel to the xy plane.
------------------------------------------------------------------------- */
void
EwaldN
::
compute_slabcorr
(
int
eflag
)
{
// compute local contribution to global dipole moment
double
*
q
=
atom
->
q
;
double
*
x
=
atom
->
x
[
0
]
-
1
,
*
xn
=
x
+
3
*
atom
->
nlocal
-
1
;
double
dipole
=
0.0
,
dipole_all
;
while
((
x
+=
3
)
<
xn
)
dipole
+=
*
x
*
*
(
q
++
);
MPI_Allreduce
(
&
dipole
,
&
dipole_all
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
const
double
qscale
=
force
->
qqrd2e
*
scale
;
double
ffact
=
-
4.0
*
M_PI
*
qscale
*
dipole_all
/
volume
;
double
*
f
=
atom
->
f
[
0
]
-
1
,
*
fn
=
f
+
3
*
atom
->
nlocal
-
3
;
q
=
atom
->
q
;
while
((
f
+=
3
)
<
fn
)
*
f
+=
ffact
*
*
(
q
++
);
if
(
eflag
)
// energy correction
energy
+=
qscale
*
(
2.0
*
M_PI
*
dipole_all
*
dipole_all
/
volume
);
}
Event Timeline
Log In to Comment