QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Enable separation of long- and short-range of Ewald sum in python

This commit is contained in:
Marek Nečada 2020-07-28 09:39:22 +03:00
parent 5fdbb362a3
commit 8abe7e3357
7 changed files with 68 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
qpms/cycommon.pyx
View File

@ -20,6 +20,12 @@ class BesselType(enum.IntEnum):
HANKEL_PLUS = QPMS_HANKEL_PLUS
HANKEL_MINUS = QPMS_HANKEL_MINUS
class EwaldPart(enum.IntEnum):
LONG_RANGE = QPMS_EWALD_LONG_RANGE
SHORT_RANGE = QPMS_EWALD_SHORT_RANGE
FULL = QPMS_EWALD_FULL
ZEROTERM = QPMS_EWALD_0TERM
class PointGroupClass(enum.IntEnum):
CN = QPMS_PGS_CN
S2N = QPMS_PGS_S2N

8
qpms/ewald.h
View File

@ -36,14 +36,6 @@
#include "lattices.h"
typedef enum {
QPMS_EWALD_LONG_RANGE = 1,
QPMS_EWALD_SHORT_RANGE = 2,
QPMS_EWALD_0TERM = 4,
QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
} qpms_ewald_part;
/// Use this handler to ignore underflows of incomplete gamma.
gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;

13
qpms/qpms_c.pyx
View File

@ -11,7 +11,7 @@ from .qpms_cdefs cimport *
from .cyquaternions cimport IRot3, CQuat
from .cybspec cimport BaseSpec
from .cycommon cimport make_c_string
from .cycommon import string_c2py, PointGroupClass, BesselType
from .cycommon import string_c2py, PointGroupClass, BesselType, EwaldPart
from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator
from .cymaterials cimport EpsMuGenerator, EpsMu
from libc.stdlib cimport malloc, free, calloc
@ -1026,7 +1026,8 @@ cdef class _ScatteringSystemAtOmegaK:
self.sswk.eta = eta
@boundscheck(False)
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS):
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS,
ewaldparts = EwaldPart.FULL):
"""Evaluate electric field for a given excitation coefficient vector (periodic system)
Parameters
@ -1045,6 +1046,7 @@ cdef class _ScatteringSystemAtOmegaK:
if(btyp != QPMS_HANKEL_PLUS):
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
cdef qpms_bessel_t btyp_c = BesselType(btyp)
cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3")
@ -1060,14 +1062,14 @@ cdef class _ScatteringSystemAtOmegaK:
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
res = qpms_scatsyswk_scattered_E(&self.sswk, btyp_c, &scv_view[0], pos)
res = qpms_scatsyswk_scattered_E_e(&self.sswk, btyp_c, &scv_view[0], pos, ewaldparts_c)
results[i,0] = res.x
results[i,1] = res.y
results[i,2] = res.z
return results.reshape(evalpos.shape)
@boundscheck(False)
def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS):
def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
# TODO examples
"""Evaluate scattered field "basis" (periodic system)
@ -1089,6 +1091,7 @@ cdef class _ScatteringSystemAtOmegaK:
if(btyp != QPMS_HANKEL_PLUS):
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
cdef qpms_bessel_t btyp_c = BesselType(btyp)
cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
cdef Py_ssize_t fecv_size = self.fecv_size
evalpos = np.array(evalpos, dtype=float, copy=False)
if evalpos.shape[-1] != 3:
@ -1104,7 +1107,7 @@ cdef class _ScatteringSystemAtOmegaK:
pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2]
qpms_scatsyswk_scattered_field_basis(res, &self.sswk, btyp_c, pos)
qpms_scatsyswk_scattered_field_basis_e(res, &self.sswk, btyp_c, pos, ewaldparts_c)
for j in range(fecv_size):
results[i,j,0] = res[j].x
results[i,j,1] = res[j].y

15
qpms/qpms_cdefs.pxd
View File

@ -150,6 +150,11 @@ cdef extern from "qpms_types.h":
cdouble mu
ctypedef enum qpms_coord_system_t:
pass
ctypedef enum qpms_ewald_part:
QPMS_EWALD_LONG_RANGE
QPMS_EWALD_SHORT_RANGE
QPMS_EWALD_FULL
QPMS_EWALD_0TERM
# maybe more if needed
cdef extern from "qpms_error.h":
@ -708,6 +713,8 @@ cdef extern from "scatsystem.h":
const cdouble *f_excitation_vector_full, cart3_t where) nogil
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
const cdouble *f_excitation_vector_full, cart3_t where) nogil
ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
const cdouble *f_excitation_vector_full, cart3_t where, qpms_ewald_part parts) nogil
qpms_errno_t qpms_scatsys_scattered_field_basis(ccart3_t *target, const qpms_scatsys_t *ss,
qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
qpms_errno_t qpms_scatsysw_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
@ -718,6 +725,8 @@ cdef extern from "scatsystem.h":
qpms_ss_pi_t pi, qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, cart3_t where) nogil
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, cart3_t where, qpms_ewald_part parts) nogil
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil
@ -726,12 +735,6 @@ cdef extern from "ewald.h":
cdouble val
double err
ctypedef enum qpms_ewald_part:
QPMS_EWALD_LONG_RANGE
QPMS_EWALD_SHORT_RANGE
QPMS_EWALD_FULL
QPMS_EWALD_0TERM
struct qpms_ewald3_constants_t:
qpms_l_t lMax
qpms_y_t nelem_sc

7
qpms/qpms_types.h
View File

@ -426,5 +426,12 @@ typedef struct qpms_epsmu_t {
struct qpms_tolerance_spec_t; // See tolerances.h
typedef enum {
QPMS_EWALD_LONG_RANGE = 1,
QPMS_EWALD_SHORT_RANGE = 2,
QPMS_EWALD_0TERM = 4,
QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
} qpms_ewald_part;
#define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
#endif // QPMS_TYPES

30
qpms/scatsystem.c
View File

@ -2179,10 +2179,11 @@ ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
// For periodic lattices, we use directly the "alternative" implementation,
// using translation operator and regular dipole waves at zero
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp,
const complex double *cvf,
const cart3_t where
const cart3_t where,
const qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
@ -2214,13 +2215,13 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
ss->c, s, NULL,
&dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK));
maxR, maxK, parts));
for(size_t i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < dipspec.n; ++j){
@ -2234,11 +2235,17 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
return res;
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis(
ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t btyp, const complex double *cvf, const cart3_t where) {
return qpms_scatsyswk_scattered_E_e(sswk, btyp, cvf, where, QPMS_EWALD_FULL);
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
ccart3_t *target,
const qpms_scatsys_at_omega_k_t *sswk,
const qpms_bessel_t btyp,
const cart3_t where
const cart3_t where,
const qpms_ewald_part parts
) {
QPMS_UNTESTED;
if (btyp != QPMS_HANKEL_PLUS)
@ -2272,13 +2279,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
const double maxR = sqrt(ss->per.unitcell_volume) * 64;
const double maxK = 2048 * 2 * M_PI / maxR;
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
ss->c, s, NULL,
&dipspec, 1, bspec, dipspec.n,
sswk->eta, sswk->ssw->wavenumber,
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/,
maxR, maxK));
maxR, maxK, parts));
for(size_t i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < dipspec.n; ++j){
@ -2294,6 +2301,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
return QPMS_SUCCESS;
}
qpms_errno_t qpms_scatsyswk_scattered_field_basis(
ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
const qpms_bessel_t btyp, const cart3_t where) {
return qpms_scatsyswk_scattered_field_basis_e(
target, sswk, btyp, where, QPMS_EWALD_FULL);
}
#if 0
ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
qpms_iri_t iri, const complex double *cvr, cart3_t where) {

16
qpms/scatsystem.h
View File

@ -892,6 +892,14 @@ ccart3_t qpms_scatsyswk_scattered_E(
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
ccart3_t qpms_scatsyswk_scattered_E_e(
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the field is evaluated.
qpms_ewald_part parts
);
/// Evaluates "scattered" field basis functions in a periodic system.
/**
*
@ -904,6 +912,14 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the basis is evaluated.
);
qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
ccart3_t *target, ///< Target array of size sswk->ssw->ss->fecv_size
const qpms_scatsys_at_omega_k_t *sswk,
qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
qpms_ewald_part parts
);
/// Adjusted Ewadl parameter to avoid high-frequency breakdown.
// TODO DOC