fennol.utils.kspace

View Source

 1import numpy as np
 2import math
 3
 4def get_reciprocal_space_parameters(reciprocal_cells, cutoff, kmax=30, kthr=1.0e-6):
 5    # find optimal ewald parameters (preprocessing)
 6    eps = 1.0e-8
 7    ratio = eps + 1
 8    x = 0.5
 9    i = 0
10    # approximate value
11    while ratio > eps:
12        x *= 2
13        ratio = math.erfc(x * cutoff) / cutoff
14    # refine with binary search
15    k = i + 60
16    xlo = 0.0
17    xhi = x
18    for i in range(1, k + 1):
19        x = (xlo + xhi) / 2.0
20        ratio = math.erfc(x * cutoff) / cutoff
21        if ratio > eps:
22            xlo = x
23        else:
24            xhi = x
25    bewald = x
26
27    # set k points
28    kxs = np.arange(kmax + 1)
29    kxs = np.concatenate((kxs, -kxs[1:]))
30    k = np.array(np.meshgrid(kxs, kxs, kxs)).reshape(3, -1).T[1:]
31    # set exp factor
32    m2s = []
33    expfacs = []
34    ks = []
35    nks = []
36    for i, A in enumerate(range(reciprocal_cells.shape[0])):
37        A = reciprocal_cells[i].T
38        m2 = np.sum(
39            (
40                k[:, 0, None] * A[None, 0, :]
41                + k[:, 1, None] * A[None, 1, :]
42                + k[:, 2, None] * A[None, 2, :]
43            )
44            ** 2,
45            axis=-1,
46        )
47        a2 = (np.pi / bewald) ** 2
48        expfac = np.exp(-a2 * m2) / m2
49        isort = np.argsort(expfac)[::-1]
50        expfac = expfac[isort]
51        m2 = m2[isort]
52        ki = k[isort]
53        sel = (expfac > kthr).nonzero()[0]
54        nks.append(len(sel))
55        m2s.append(m2)
56        expfacs.append(expfac)
57        ks.append(ki)
58
59    ks = np.stack(ks, axis=0)
60    m2s = np.stack(m2s, axis=0)
61    expfacs = np.stack(expfacs, axis=0)
62    nks = np.array(nks, dtype=np.int64)
63    nk = np.max(nks)
64    return ks[:, :nk, :], expfacs[:, :nk], m2s[:, :nk], bewald

def get_reciprocal_space_parameters(reciprocal_cells, cutoff, kmax=30, kthr=1e-06): View Source

 5def get_reciprocal_space_parameters(reciprocal_cells, cutoff, kmax=30, kthr=1.0e-6):
 6    # find optimal ewald parameters (preprocessing)
 7    eps = 1.0e-8
 8    ratio = eps + 1
 9    x = 0.5
10    i = 0
11    # approximate value
12    while ratio > eps:
13        x *= 2
14        ratio = math.erfc(x * cutoff) / cutoff
15    # refine with binary search
16    k = i + 60
17    xlo = 0.0
18    xhi = x
19    for i in range(1, k + 1):
20        x = (xlo + xhi) / 2.0
21        ratio = math.erfc(x * cutoff) / cutoff
22        if ratio > eps:
23            xlo = x
24        else:
25            xhi = x
26    bewald = x
27
28    # set k points
29    kxs = np.arange(kmax + 1)
30    kxs = np.concatenate((kxs, -kxs[1:]))
31    k = np.array(np.meshgrid(kxs, kxs, kxs)).reshape(3, -1).T[1:]
32    # set exp factor
33    m2s = []
34    expfacs = []
35    ks = []
36    nks = []
37    for i, A in enumerate(range(reciprocal_cells.shape[0])):
38        A = reciprocal_cells[i].T
39        m2 = np.sum(
40            (
41                k[:, 0, None] * A[None, 0, :]
42                + k[:, 1, None] * A[None, 1, :]
43                + k[:, 2, None] * A[None, 2, :]
44            )
45            ** 2,
46            axis=-1,
47        )
48        a2 = (np.pi / bewald) ** 2
49        expfac = np.exp(-a2 * m2) / m2
50        isort = np.argsort(expfac)[::-1]
51        expfac = expfac[isort]
52        m2 = m2[isort]
53        ki = k[isort]
54        sel = (expfac > kthr).nonzero()[0]
55        nks.append(len(sel))
56        m2s.append(m2)
57        expfacs.append(expfac)
58        ks.append(ki)
59
60    ks = np.stack(ks, axis=0)
61    m2s = np.stack(m2s, axis=0)
62    expfacs = np.stack(expfacs, axis=0)
63    nks = np.array(nks, dtype=np.int64)
64    nk = np.max(nks)
65    return ks[:, :nk, :], expfacs[:, :nk], m2s[:, :nk], bewald