Evaluate scattered electric fields in 2d-periodic system.
Former-commit-id: 36386215f9d330a3047cb9a294ccc1de55f2121f
This commit is contained in:
Marek Nečada
parent
3046f03734
commit
6009de6fa2
11
TODO.md
11
TODO.md
|
|
@ -1,11 +1,12 @@
|
|||
TODO list before public release
|
||||
===============================
|
||||
TODO list before 1.0 release
|
||||
============================
|
||||
|
||||
- Tests!
|
||||
- Docs!
|
||||
- Cross section calculations.
|
||||
- Field calculations.
|
||||
- Complex frequencies, n's, k's.
|
||||
- Cross section calculations. (Done in some Python scripts.)
|
||||
- Field calculations. (Partly done, needs more testing.)
|
||||
* Also test periodic vs. nonperiodic consistence (big finite lattice + absorbing medium vs. infinite lattice + absorbing medium).
|
||||
- Complex frequencies, n's, k's. (Mostly done.)
|
||||
- Transforming point (meta)generators.
|
||||
- Check whether moble's quaternions and my
|
||||
quaternions give the same results in tmatrices.py
|
||||
|
|
|
|||
43
qpms/ewald.c
43
qpms/ewald.c
|
|
@ -158,7 +158,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
|
|||
// ----- New generation 2D-in-3D constants ------
|
||||
// TODO it is not necessary to treat +|m| and -|m| cases separately
|
||||
// N.B. currently, this is only valid for EWALD32_CONSTANTS_AGNOSTIC (NOT CHECKED!)
|
||||
c->S1_constfacs = malloc(c->nelem_sc * sizeof(complex double *));
|
||||
c->S1_constfacs = malloc((1+c->nelem_sc) * sizeof(complex double *));
|
||||
//determine sizes
|
||||
size_t S1_constfacs_sz = 0;
|
||||
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
|
||||
|
|
@ -173,13 +173,14 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
c->S1_constfacs_base = malloc(S1_constfacs_sz * sizeof(complex double));
|
||||
size_t S1_constfacs_sz_cumsum = 0; // second count
|
||||
for (qpms_y_t y = 0; y < c->nelem_sc; ++y) {
|
||||
qpms_l_t n; qpms_m_t m; qpms_y2mn_sc_p(y, &m, &n);
|
||||
const complex double yfactor = -2 * ipow(n+1) * M_SQRTPI
|
||||
* factorial((n-m)/2) * factorial((n+m)/2);
|
||||
c->S1_constfacs[y] = c->s1_constfacs_base + S1_constfacs_sz_cumsum;
|
||||
c->S1_constfacs[y] = c->S1_constfacs_base + S1_constfacs_sz_cumsum;
|
||||
size_t coeffs_per_y = 0;
|
||||
const qpms_l_t L_M = n - abs(m);
|
||||
for(qpms_l_t j = 0; j <= L_M; ++j) { // outer sum
|
||||
|
|
@ -197,6 +198,7 @@ qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, cons
|
|||
S1_constfacs_sz_cumsum += coeffs_per_y;
|
||||
}
|
||||
QPMS_ASSERT(S1_constfacs_sz_cumsum = S1_constfacs_sz);
|
||||
c->S1_constfacs[c->nelem_sc] = c->S1_constfacs_base + S1_constfacs_sz; // For easier limit checks
|
||||
|
||||
// ------ the "z-axis constants" -----
|
||||
// determine sizes
|
||||
|
|
@ -246,6 +248,8 @@ void qpms_ewald3_constants_free(qpms_ewald3_constants_t *c) {
|
|||
free(c->legendre_minus1);
|
||||
free(c->s1_constfacs);
|
||||
free(c->s1_constfacs_base);
|
||||
free(c->S1_constfacs);
|
||||
free(c->S1_constfacs_base);
|
||||
free(c->s1_constfacs_1Dz_base);
|
||||
free(c->s1_constfacs_1Dz);
|
||||
free(c->s1_jMaxes);
|
||||
|
|
@ -347,7 +351,7 @@ int ewald3_21_xy_sigma_long (
|
|||
double Gamma_pq_err[lMax/2+1];
|
||||
|
||||
// CHOOSE POINT BEGIN
|
||||
// TODO mayby PGen_next_sph is not the best coordinate system choice here
|
||||
// TODO maybe PGen_next_sph is not the best coordinate system choice here
|
||||
while ((pgen_retdata = PGen_next_sph(pgen_K)).flags & PGEN_NOTDONE) { // BEGIN POINT LOOP
|
||||
cart3_t K_pq_cart;
|
||||
sph_t beta_pq_sph;
|
||||
|
|
@ -407,12 +411,14 @@ int ewald3_21_xy_sigma_long (
|
|||
// and just fetched for each n, m pair
|
||||
for(qpms_l_t n = 0; n <= lMax; ++n)
|
||||
for(qpms_m_t m = -n; m <= n; ++m) {
|
||||
if((m+n) % 2 != 0) // odd coefficients are zero.
|
||||
if((particle_shift.z == 0) && ((m+n) % 2 != 0)) // odd coefficients are zero.
|
||||
continue;
|
||||
const qpms_y_t y = qpms_mn2y_sc(m, n);
|
||||
size_t constidx = 0; // constants offset
|
||||
const complex double e_imalpha_pq = cexp(I*m*arg_pq);
|
||||
complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c);
|
||||
double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors?
|
||||
if (particle_shift.z == 0) { // TODO remove when the general case is stable and tested
|
||||
assert((n-abs(m))/2 == c->s1_jMaxes[y]);
|
||||
for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
|
||||
complex double summand = cpow_0lim_zi(rbeta_pq/k, n-2*j)
|
||||
|
|
@ -432,6 +438,35 @@ int ewald3_21_xy_sigma_long (
|
|||
kahanadd(&lsum, &lsum_c, cabs(jsum));
|
||||
#endif
|
||||
if(err) kahanadd(err + y, err_c + y, jsum_err);
|
||||
} else { // particle_shift.z != 0
|
||||
const qpms_l_t L_M = n - abs(m);
|
||||
for(qpms_l_t j = 0; j <= L_M; ++j) { // outer sum
|
||||
complex double ssum, ssum_c; ckahaninit(&ssum, &ssum_c);
|
||||
// TODO errors of ssum
|
||||
// inner sum: j <= s <= min(2*j, n - |m|), s has the same parity as n - |m|
|
||||
for(qpms_l_t s = j + (L_M - j) % 2;
|
||||
(s <= 2 * j) && (s <= L_M);
|
||||
s += 2) {
|
||||
complex double ssummand = c->S1_constfacs[y][constidx]
|
||||
* cpow(-k * particle_shift.z, 2*j - s) * cpow_0lim_zi(rbeta_pq / k, n - s);
|
||||
ckahanadd(&ssum, &ssum_c, ssummand);
|
||||
++constidx;
|
||||
}
|
||||
const complex double jfactor = e_imalpha_pq * Gamma_pq[j] * cpow(gamma_pq, 2*j - 1);
|
||||
if (err) { // FIXME include also other sources of error than Gamma_pq's relative error
|
||||
double jfactor_err = Gamma_pq_err[j] * pow(cabs(gamma_pq), 2*j - 1);
|
||||
kahanadd(&jsum_err, &jsum_err_c, jfactor_err * ssum);
|
||||
}
|
||||
complex double jsummand = jfactor * ssum;
|
||||
ckahanadd(&jsum, &jsum_c, jsummand);
|
||||
}
|
||||
jsum *= phasefac; // factor1d not here, off-axis sums not implemented/allowed.
|
||||
ckahanadd(target + y, target_c + y, jsum);
|
||||
#ifdef EWALD_AUTO_CUTOFF
|
||||
kahanadd(&lsum, &lsum_c, cabs(jsum));
|
||||
#endif
|
||||
if(err) kahanadd(err + y, err_c + y, jsum_err);
|
||||
}
|
||||
}
|
||||
#ifndef NDEBUG
|
||||
rbeta_pq_prev = rbeta_pq;
|
||||
|
|
|
|||
|
|
@ -971,6 +971,30 @@ cdef class _ScatteringSystemAtOmegaK:
|
|||
def __set__(self, double eta):
|
||||
self.sswk.eta = eta
|
||||
|
||||
def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS): # TODO DOC!!!
|
||||
if(btyp != QPMS_HANKEL_PLUS):
|
||||
raise NotImplementedError("Only first kind Bessel function-based fields are supported")
|
||||
cdef qpms_bessel_t btyp_c = BesselType(btyp)
|
||||
evalpos = np.array(evalpos, dtype=float, copy=False)
|
||||
if evalpos.shape[-1] != 3:
|
||||
raise ValueError("Last dimension of evalpos has to be 3")
|
||||
cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
|
||||
cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False)
|
||||
cdef cdouble[::1] scv_view = scv
|
||||
cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
|
||||
cdef ccart3_t res
|
||||
cdef cart3_t pos
|
||||
cdef size_t i
|
||||
for i in range(evalpos_a.shape[0]):
|
||||
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)
|
||||
results[i,0] = res.x
|
||||
results[i,1] = res.y
|
||||
results[i,2] = res.z
|
||||
return results.reshape(evalpos.shape)
|
||||
|
||||
cdef class _ScatteringSystemAtOmega:
|
||||
'''
|
||||
Wrapper over the C qpms_scatsys_at_omega_t structure
|
||||
|
|
|
|||
|
|
@ -703,6 +703,8 @@ cdef extern from "scatsystem.h":
|
|||
const cdouble *f_excitation_vector_full, cart3_t where)
|
||||
ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t btyp,
|
||||
const cdouble *f_excitation_vector_full, cart3_t where)
|
||||
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)
|
||||
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector);
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1232,7 +1232,7 @@ static inline int qpms_ss_ppair_W32xy(const qpms_scatsys_t *ss,
|
|||
eta, wavenumber,
|
||||
cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
|
||||
kvector,
|
||||
cart2_substract(cart3xy2cart2(ss->p[pdest].pos), cart3xy2cart2(ss->p[psrc].pos)),
|
||||
cart3_substract(ss->p[pdest].pos, ss->p[psrc].pos),
|
||||
maxR, maxK, parts);
|
||||
}
|
||||
|
||||
|
|
@ -2039,6 +2039,7 @@ ccart3_t qpms_scatsys_scattered_E(const qpms_scatsys_t *ss,
|
|||
ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
|
||||
qpms_bessel_t btyp,
|
||||
const complex double *cvf, const cart3_t where) {
|
||||
qpms_ss_ensure_nonperiodic_a(ssw->ss, "qpms_scatsyswk_scattered_E()");
|
||||
return qpms_scatsys_scattered_E(ssw->ss, btyp, ssw->wavenumber,
|
||||
cvf, where);
|
||||
}
|
||||
|
|
@ -2051,6 +2052,7 @@ ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss,
|
|||
const cart3_t where
|
||||
) {
|
||||
QPMS_UNTESTED;
|
||||
qpms_ss_ensure_nonperiodic(ss);
|
||||
ccart3_t res = {0,0,0};
|
||||
ccart3_t res_kc = {0,0,0}; // kahan sum compensation
|
||||
|
||||
|
|
@ -2108,6 +2110,83 @@ ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
|
|||
cvf, where);
|
||||
}
|
||||
|
||||
// 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,
|
||||
qpms_bessel_t btyp,
|
||||
const complex double *cvf,
|
||||
const cart3_t where
|
||||
) {
|
||||
QPMS_UNTESTED;
|
||||
if (btyp != QPMS_HANKEL_PLUS)
|
||||
QPMS_NOT_IMPLEMENTED("Only scattered field with first kind Hankel functions currently implemented.");
|
||||
const qpms_scatsys_t *ss = sswk->ssw->ss;
|
||||
if (ss->lattice_dimension != 2)
|
||||
QPMS_NOT_IMPLEMENTED("Only 2D-periodic lattices implemented");
|
||||
ccart3_t res = {0,0,0};
|
||||
ccart3_t res_kc = {0,0,0}; // kahan sum compensation
|
||||
|
||||
static const int dipspecn = 3; // We have three basis vectors
|
||||
// bspec containing only electric dipoles
|
||||
const qpms_vswf_set_spec_t dipspec = {
|
||||
.n = dipspecn,
|
||||
.ilist = (qpms_uvswfi_t[]){
|
||||
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1),
|
||||
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1),
|
||||
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1),
|
||||
},
|
||||
.lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1,
|
||||
.capacity=0,
|
||||
.norm = ss->c->normalisation,
|
||||
};
|
||||
|
||||
ccart3_t regdipoles_0[dipspecn]; {
|
||||
const sph_t origin_sph = {.r = 0, .theta = M_PI_2, .phi=0}; // Should work with any theta/phi (TESTWORTHY)
|
||||
csphvec_t regdipoles_0_sph[dipspecn];
|
||||
QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(regdipoles_0_sph, &dipspec,
|
||||
sph2csph(origin_sph), QPMS_BESSEL_REGULAR));
|
||||
for(int i = 0; i < dipspecn; ++i)
|
||||
regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph);
|
||||
}
|
||||
|
||||
complex double *s; // Translation matrix
|
||||
QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n);
|
||||
|
||||
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
|
||||
const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
|
||||
const cart3_t particle_pos = ss->p[pi].pos;
|
||||
const complex double *particle_cv = cvf + ss->fecv_pstarts[pi];
|
||||
|
||||
const cart3_t origin_cart = {0, 0, 0};
|
||||
|
||||
QPMS_ASSERT(sswk->k[2] == 0); // At least not implemented now
|
||||
QPMS_ASSERT(ss->per.lattice_basis[0].z == 0);
|
||||
QPMS_ASSERT(ss->per.lattice_basis[1].z == 0);
|
||||
|
||||
// Same choices as in qpms_ss_ppair_W32xy; TODO make it more dynamic
|
||||
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(
|
||||
ss->c, s, NULL,
|
||||
&dipspec, 1, bspec, dipspecn,
|
||||
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));
|
||||
|
||||
for(size_t i = 0; i < bspec->n; ++i)
|
||||
for(size_t j = 0; j < dipspecn; ++j){
|
||||
ccart3_t summand = ccart3_scale(particle_cv[i] * s[dipspecn*i+j], regdipoles_0[j]);
|
||||
ckahanadd(&(res.x), &(res_kc.x), summand.x);
|
||||
ckahanadd(&(res.y), &(res_kc.y), summand.y);
|
||||
ckahanadd(&(res.z), &(res_kc.z), summand.z);
|
||||
}
|
||||
}
|
||||
free(s);
|
||||
return res;
|
||||
}
|
||||
|
||||
#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) {
|
||||
|
|
|
|||
|
|
@ -708,7 +708,9 @@ complex double *qpms_scatsys_incident_field_vector_irrep_packed(
|
|||
*
|
||||
* \see qpms_scatsysw_scattered_E()
|
||||
*
|
||||
* Not yet implemented for periodic systems.
|
||||
* \see qpms_scatsyswk_scattered_E() for periodic systems.
|
||||
*
|
||||
*
|
||||
*/
|
||||
ccart3_t qpms_scatsys_scattered_E(
|
||||
const qpms_scatsys_t *ss,
|
||||
|
|
@ -727,7 +729,7 @@ ccart3_t qpms_scatsys_scattered_E(
|
|||
*
|
||||
* \see qpms_scatsys_scattered_E()
|
||||
*
|
||||
* Not yet implemented for periodic systems.
|
||||
* \see qpms_scatsyswk_scattered_E() for periodic systems.
|
||||
*/
|
||||
ccart3_t qpms_scatsysw_scattered_E(
|
||||
const qpms_scatsys_at_omega_t *ssw,
|
||||
|
|
@ -769,6 +771,25 @@ ccart3_t qpms_scatsysw_scattered_E__alt(
|
|||
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
|
||||
);
|
||||
|
||||
|
||||
/// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients. (Periodic system.)
|
||||
/**
|
||||
* This function evaluates the field \f$ \vect E (\vect r ) \f$, with any
|
||||
* given vector of the excitation coefficients \f$ \wckcout \f$.
|
||||
*
|
||||
* \return Complex electric field at the point defined by \a evalpoint.
|
||||
*
|
||||
* \bug Currently implemented only for btyp == QPMS_HANKEL_PLUS.
|
||||
*
|
||||
* \see qpms_scatsys_scattered_E(), qpms_scatsysw_scattered_E() for finite systems.
|
||||
*/
|
||||
ccart3_t qpms_scatsyswk_scattered_E(
|
||||
const qpms_scatsys_at_omega_k_t *sswk,
|
||||
qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS).
|
||||
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.
|
||||
);
|
||||
|
||||
/// Adjusted Ewadl parameter to avoid high-frequency breakdown.
|
||||
// TODO DOC
|
||||
double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]);
|
||||
|
|
|
|||
|
|
@ -730,7 +730,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul
|
|||
const double eta, const complex double k,
|
||||
cart2_t b1, cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK,
|
||||
const qpms_ewald_part parts
|
||||
)
|
||||
|
|
@ -784,7 +784,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
|
|||
const double eta, const complex double k,
|
||||
cart2_t b1, cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK
|
||||
)
|
||||
{
|
||||
|
|
@ -931,7 +931,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
|
|||
const double eta, const complex double k,
|
||||
const cart2_t b1, const cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK,
|
||||
const qpms_ewald_part parts
|
||||
)
|
||||
|
|
@ -940,7 +940,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
|
|||
const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax);
|
||||
//const qpms_y_t nelem = qpms_lMax2nelem(c->lMax);
|
||||
const bool doerr = Aerr || Berr;
|
||||
const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
|
||||
const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector
|
||||
|
||||
complex double *sigmas_short = malloc(sizeof(complex double)*nelem2_sc);
|
||||
complex double *sigmas_long = malloc(sizeof(complex double)*nelem2_sc);
|
||||
|
|
@ -972,7 +972,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
|
|||
#else
|
||||
false,
|
||||
#endif
|
||||
cart22cart3xy(beta), cart22cart3xy(particle_shift)));
|
||||
cart22cart3xy(beta), particle_shift));
|
||||
if(Kgen.stateData) // PGen not consumed entirely (converged earlier)
|
||||
PGen_destroy(&Kgen);
|
||||
}
|
||||
|
|
@ -980,11 +980,18 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
|
|||
if (parts & QPMS_EWALD_SHORT_RANGE) {
|
||||
PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL,
|
||||
#ifdef GEN_RSHIFTEDPOINTS
|
||||
cart2_scale(-1 /*CHECKSIGN*/, particle_shift),
|
||||
cart2_scale(-1 /*CHECKSIGN*/, cart3xy2cart2(particle_shift)),
|
||||
#else
|
||||
CART2_ZERO,
|
||||
#endif
|
||||
0, !do_sigma0, maxR, false);
|
||||
#ifdef GEN_RSHIFTEDPOINTS // rather ugly hacks, LPTODO cleanup
|
||||
if (particle_shift.z != 0) {
|
||||
const cart3_t zshift = {0, 0, -particle_shift.z /*CHECKSIGN*/};
|
||||
Rgen = Pgen_shifted_new(Rgen, zshift);
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k,
|
||||
LAT_2D_IN_3D_XYONLY, &Rgen,
|
||||
|
|
@ -993,7 +1000,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
|
|||
#else
|
||||
false,
|
||||
#endif
|
||||
cart22cart3xy(beta), cart22cart3xy(particle_shift)));
|
||||
cart22cart3xy(beta), particle_shift));
|
||||
|
||||
if(Rgen.stateData) // PGen not consumed entirely (converged earlier)
|
||||
PGen_destroy(&Rgen);
|
||||
|
|
@ -1078,7 +1085,7 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
|
|||
const double eta, const complex double k,
|
||||
const cart2_t b1, const cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK)
|
||||
{
|
||||
return qpms_trans_calculator_get_AB_arrays_e32_e(
|
||||
|
|
|
|||
|
|
@ -159,10 +159,10 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
|
|||
complex double *Adest, double *Aerr,
|
||||
complex double *Bdest, double *Berr,
|
||||
const ptrdiff_t deststride, const ptrdiff_t srcstride,
|
||||
const double eta, const complex double k,
|
||||
const double eta, const complex double wavenumber,
|
||||
cart2_t b1, cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK
|
||||
);
|
||||
|
||||
|
|
@ -170,10 +170,10 @@ int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
|
|||
complex double *Adest, double *Aerr,
|
||||
complex double *Bdest, double *Berr,
|
||||
const ptrdiff_t deststride, const ptrdiff_t srcstride,
|
||||
const double eta, const complex double k,
|
||||
const double eta, const complex double wavenumber,
|
||||
cart2_t b1, cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK,
|
||||
qpms_ewald_part parts
|
||||
);
|
||||
|
|
@ -185,10 +185,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
|
|||
const qpms_vswf_set_spec_t *destspec, size_t deststride,
|
||||
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
|
||||
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
|
||||
const double eta, const complex double k,
|
||||
const double eta, const complex double wavenumber,
|
||||
cart2_t b1, cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK
|
||||
);
|
||||
|
||||
|
|
@ -198,10 +198,10 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calcul
|
|||
const qpms_vswf_set_spec_t *destspec, size_t deststride,
|
||||
/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
|
||||
const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
|
||||
const double eta, const complex double k,
|
||||
const double eta, const complex double wavenumber,
|
||||
cart2_t b1, cart2_t b2,
|
||||
const cart2_t beta,
|
||||
const cart2_t particle_shift,
|
||||
const cart3_t particle_shift,
|
||||
double maxR, double maxK,
|
||||
qpms_ewald_part parts
|
||||
);
|
||||
|
|
|
|||
Loading…
Reference in New Issue