This is the full computation of the residual bias correction for our (mean-field) GF solution for the percolation problem (immobile "solute" with vacancy diffusion). It work through all of the matrix averages "analytically" (storing them as polynomials in the concentration of the immobile solute, $c_\text{B}$), and then brings everything together to express the residual bias correction as an analytic function with numerical coefficients for the square lattice.
In [365]:
import sys
sys.path.extend(['.','./Vacancy'])
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
%matplotlib inline
import scipy.sparse
import itertools
from numba import jit, njit, prange, guvectorize # faster runtime with update routines
from scipy.misc import comb
# from sympy import *
import onsager.PowerExpansion as PE
import onsager.crystal as crystal
import onsager.crystalStars as stars
import onsager.GFcalc as GFcalc
from tqdm import tnrange, tqdm_notebook
In [ ]:
# Turn off or on to run optional testing code in notebook:
# Also turns on / off progress bars
__TESTING__ = False
Now, we need to expand out our probability factors. Let $x$ be the concentration of solute B; imagine we have $N$ sites possible. Then, if there are $n$ B atoms, the probability factor is
$$P(n;N) = x^n (1-x)^{N-n} = x^n \sum_{j=0}^{N-n} \frac{(N-n)!}{j!(N-n-j)!} (-x)^j
= \sum_{j=0}^{N-n} \frac{(N-n)!}{j!(N-n-j)!} (-1)^j x^{n+j}$$
The factorial term is $N-n$ choose $j$, which is scipy.misc.comb.
We want to construct a probability matrix P[n,c] such that $P(n;N)$ is written as a sum over $x^c$ terms; $c=0\ldots N$.
In [2]:
def calc_P(N):
"""
Returns the probability matrix P[n,c] where the probability of seeing `n` atoms
of type B in `N` sites is sum(c=0..N, x^c P[n,c])
:param N: total number of sites
:returns P[n,c]: matrix of probablities, n=0..N, c=0..N
"""
P = np.zeros((N+1, N+1), dtype=int)
for n in range(N+1):
Nn = N-n
P[n,n:] = comb([Nn]*(Nn+1), [j for j in range(Nn+1)])
for j in range(Nn+1):
P[n,j+n] *= (-1)**j
return P
In [3]:
if __TESTING__:
calc_P(4)
Out[3]:
Normalization check: construct the $2^N$ states, and see if it averages to 1. Each state is a vector of length $N$, with entries that are 0 (A) or 1 (B). Here, we explicitly build our state space, and also do a quick count to determine $n_\text{B}$ for each state. Note: we prepend a value of 0, since this corresponds to the initial location of the vacancy.
New version: we now generate group operations for the square lattice, and take advantage of those to reduce the computational time.
In [4]:
N = 24
prob = calc_P(N)
states = np.array([(0,) + st for st in itertools.product((0,1), repeat=N)])
nB = np.sum(states, axis=1)
In [5]:
if __TESTING__:
norm = np.zeros(N+1, dtype=int)
for n in tqdm_notebook(nB):
norm += prob[n]
print(norm)
In [6]:
states.shape
Out[6]:
In [7]:
Pstates = np.array([prob[n] for n in nB])
Now, we do some analysis by constructing up to 3 jumps (corresponding to third power of our transition rate matrix $W$). We do this analysis by setting up some bookkeeping:
[dx_0, dx_1, dx_2, dx_3]This is all sufficient to construct a sparse version of $W$ (and $\Gamma$) for a given state $\chi$.
In [8]:
dxlist = [np.array([1,0]), np.array([-1,0]), np.array([0,1]), np.array([0,-1])]
In [9]:
Njump = 3
sites = [np.array([0,0])]
sitedict = {(0,0): 0}
lastsites = sites.copy()
for nj in range(Njump):
newsites = []
for dx in dxlist:
for x in lastsites:
y = x+dx
yt = tuple(y)
if yt not in sitedict:
sitedict[yt] = len(sites)
sites.append(y)
newsites.append(y)
lastsites = newsites
Nsite = len(sites)
Nsite0 = len(sites) - len(lastsites)
sites0 = sites[:Nsite0]
In [10]:
jumplist = []
for x in sites0:
jumplist.append([sitedict[tuple(x+dx)] for dx in dxlist])
if __TESTING__:
print(jumplist)
Out[10]:
In [11]:
basisfunc, basisdict = [], {}
for x in sites:
for y in sites:
d = x-y
dt = tuple(d)
if dt not in basisdict:
basisdict[dt] = len(basisfunc)
basisfunc.append(d)
Nbasis = len(basisfunc)
Some matrices and lists to manage conversion between sites and basis functions.
We also include a matrix that corresponds to "matching" basis functions as a function of endstate $x$. This is used to correct the outer product for "missing" basis functions, for when the missing basis functions map onto identical sites.
In [12]:
chibasisfound, chibasismiss, chibasismissmatch = [], [], []
chibasisfar = []
for x in sites:
xbdict = {}
xbmiss = []
xbmissmatch = {}
xbfar = {}
for bindex, b in enumerate(basisfunc):
bt = tuple(b)
y = x+b
yt = tuple(y)
if yt in basisdict:
xbfar[bindex] = basisdict[yt]
if yt in sitedict:
xbdict[bindex] = sitedict[yt]
else:
xbmiss.append(bindex)
if bt not in sitedict and yt in basisdict:
xbmissmatch[bindex] = basisdict[yt]
chibasisfound.append(xbdict)
chibasismiss.append(xbmiss)
chibasismissmatch.append(xbmissmatch)
chibasisfar.append(xbfar)
# make a set of "outside" and "inside" basis functions:
basisout = set([tuple(basisfunc[bindex]) for bindex in chibasismiss[0]])
basisin = set([tuple(bv) for bv in basisfunc if tuple(bv) not in basisout])
In [13]:
# converting chibasisfound and chibasismiss into matrices:
chibasisfound_mat = np.zeros((N+1, Nbasis, N+1), dtype=int)
# chibasisfound_sparse = [scipy.sparse.csr_matrix((Nbasis, N+1), dtype=int)
# for n in range(N+1)]
chibasismiss_mat = np.zeros((N+1, Nbasis), dtype=int)
chibasismissmatch_mat = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
chibasisfar_mat = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
for n, cbf, cbm, cbmm, cbfar in zip(itertools.count(), chibasisfound,
chibasismiss, chibasismissmatch, chibasisfar):
for bindex in cbm:
chibasismiss_mat[n, bindex] = 1
for bindex, siteindex in cbf.items():
chibasisfound_mat[n, bindex, siteindex] = 1
# chibasisfound_sparse[n][bindex, siteindex] = 1
for bindex, siteindex in cbmm.items():
chibasismissmatch_mat[bindex, siteindex, n] = 1
for bindex, siteindex in cbfar.items():
chibasisfar_mat[bindex, siteindex, n] = 1
For our 8 group operations, corresponding to the point group operations on a square, we're going to make a reduced state list that only contains one symmetry-unique representative. This requires mapping the group operations on Cartesian coordinates into corresponding group operations on our sites, and our basis functions.
In [14]:
groupops = [np.array([[1,0],[0,1]]), np.array([[0,-1],[1,0]]),
np.array([[-1,0],[0,-1]]), np.array([[0,1],[-1,0]]),
np.array([[-1,0],[0,1]]), np.array([[1,0],[0,-1]]),
np.array([[0,-1],[-1,0]]), np.array([[0,1],[1,0]])]
In [15]:
sitegroupops, basisgroupops = [], []
for g in groupops:
sg = np.zeros([Nsite, Nsite], dtype=int)
bg = np.zeros([Nbasis, Nbasis], dtype=int)
for n, x in enumerate(sites):
yt = tuple(np.dot(g, x))
sg[sitedict[yt], n] = 1
for n, x in enumerate(basisfunc):
yt = tuple(np.dot(g, x))
bg[basisdict[yt], n] = 1
sitegroupops.append(sg)
basisgroupops.append(bg)
In [17]:
foundstates = set([])
binary = np.array([2**n for n in range(Nsite)])
symmstateslist, symmPlist = [], []
for st, P in tqdm_notebook(zip(states, Pstates), total=(2**N), disable=not __TESTING__):
bc = np.dot(st, binary)
if bc not in foundstates:
symmstateslist.append(st)
equivset = set([np.dot(np.dot(g, st), binary) for g in sitegroupops])
foundstates.update(equivset)
symmPlist.append(len(equivset)*P)
symmstates = np.array(symmstateslist)
symmPstates = np.array(symmPlist)
In [18]:
symmstates.shape
Out[18]:
In [19]:
if __TESTING__:
np.sum(symmPstates, axis=0)
Out[19]:
Now, we need symmetrized versions of a lot of our information from above, in order to properly account for all of the symmetrized versions of our basis functions. This includes
We can group these in terms of what factor of concentration goes in front.
In [20]:
biasvec_mat = np.zeros((2, N+1), dtype=int)
for j, dx in enumerate(dxlist):
biasvec_mat[:, j+1] -= dx
if __TESTING__:
print(np.dot(biasvec_mat, states[8388608]), states[8388608])
In [22]:
def symmetrize(mat, groupops0, groupops1):
"""
Designed to symmetrize the first two entries of a matrix with the
corresponding group operations
"""
symmmat = np.zeros(mat.shape)
for g0, g1 in zip(groupops0, groupops1):
for i in range(mat.shape[0]):
for j in range(mat.shape[1]):
symmmat[i, j] += np.tensordot(np.tensordot(
mat, g1[j], axes=(1,0)), g0[i], axes=(0,0))
symmmat /= len(groupops0)
return symmmat
In [23]:
@njit(nogil=True, parallel=True)
def tripleouterupdate(summand, A, B, C):
"""Update summand[i,j,k] += A[i]*B[j]*C[k]"""
I, = A.shape
J, = B.shape
K, = C.shape
for i in prange(I):
for j in prange(J):
for k in prange(K):
summand[i, j, k] += A[i]*B[j]*C[k]
In [24]:
@njit(nogil=True, parallel=True)
def matrixouterupdate(summand, A, B):
"""Update summand[i,j,k] += A[i, j]*B[k]"""
I,J = A.shape
K, = B.shape
for i in prange(I):
for j in prange(J):
for k in prange(K):
summand[i, j, k] += A[i,j]*B[k]
Let's try some averages; first, without basis functions: $\langle \tau_\chi \mathbf{b}_\chi\cdot\mathbf{b}_\chi\rangle_\chi$, the average residual bias.
In [25]:
resbiasave = np.zeros(N+1, dtype=int)
for st, P in tqdm_notebook(zip(symmstates, symmPstates), total=symmstates.shape[0],
disable=not __TESTING__):
# bv = np.sum(dx for j, dx in enumerate(dxlist) if st[j+1] == 0)
bv = np.dot(biasvec_mat, st)
W = 4-np.sum(st[1:5])
if W>0:
resbiasave += P*(bv[0]*bv[0]+bv[1]*bv[1])*(12//W)
print(resbiasave/12)
Now, an average involving a single basis function: $\langle \mathbf{b}_\chi \phi_{\chi,\mathbf{x}}\rangle_\chi$.
In [26]:
biasvecbar = np.zeros((2, Nbasis, N+1), dtype=int)
Pc = np.zeros(N+1, dtype=int)
for st, P in tqdm_notebook(zip(symmstates, symmPstates), total=symmstates.shape[0],
disable=not __TESTING__):
# bv = np.sum(dx for j, dx in enumerate(dxlist) if st[j+1] == 0)
W = 4-np.sum(st[1:5])
if W==0 or W==4: continue
bv = np.dot(biasvec_mat, st)
Pc[1:] = P[:-1]
tripleouterupdate(biasvecbar, bv, np.dot(chibasisfound_mat[0], st), P)
tripleouterupdate(biasvecbar, bv, chibasismiss_mat[0], Pc)
symmbiasvecbar = symmetrize(biasvecbar, groupops, basisgroupops)
Now, let's try a basis / basis vector average: $\langle \sum_{\chi'} \phi_{\chi,\mathbf{x}} W_{\chi\chi'} \phi_{\chi',\mathbf{y}}\rangle_\chi$.
This gets a bit complicated with the "missing" basis functions for $\chi$, and especially when we consider those that are missing in both $\chi$ and $\chi'$. We also need to treat the "far" case, where both $\mathbf{x}$ and $\mathbf{y}$ are far away from the origin.
We ignore terms higher than $c^N$ ($N$=25); no contributions are found higher than 10.
In [29]:
# @njit(nogil=True, parallel=True)
@jit
def matrixupdate(mat_bar, mat_vec, chibasis, chibasis_miss,
chibasismissmatch_mat, P, Pc, Pcc):
chibasis0, chibasis1 = chibasis[0], chibasis_miss[0]
chipbasis0, chipbasis1 = np.dot(mat_vec, chibasis), np.dot(mat_vec, chibasis_miss)
tripleouterupdate(mat_bar, chibasis0, chipbasis0, P)
tripleouterupdate(mat_bar, chibasis1, chipbasis0, Pc)
tripleouterupdate(mat_bar, chibasis0, chipbasis1, Pc)
# note: this is a little confusing; if the two ("missing") basis functions are
# referencing *different* sites, then we pick up a x^2 term; but if they
# reference the same site, it is a factor of x.
tripleouterupdate(mat_bar, chibasis1, chipbasis1, Pcc)
matchouter = np.dot(chibasismissmatch_mat, mat_vec)
matrixouterupdate(mat_bar, matchouter, Pc-Pcc)
In [30]:
# I'm not entirely sure how this is supposed to read; the matching seems to be the key?
# @njit(nogil=True, parallel=True)
@jit
def farmatrixupdate(mat_bar, mat_vec, chibasis_far, Pc, Pcc):
# note: this is a little confusing; if the two ("missing") basis functions are
# referencing *different* sites, then we pick up a x^2 term; but if they
# reference the same site, it is a factor of x.
# tripleouterupdate(mat_bar, chibasis1, chipbasis1, Pcc)
matchouter = np.dot(chibasis_far, mat_vec)
matrixouterupdate(mat_bar, matchouter, Pc-Pcc)
In [31]:
# @njit(nogil=True, parallel=True)
@jit
def vectorupdate(vec_bar, bv, vec, chibasis, chibasis_miss, P, Pc):
# chibasis0, chibasis1 = chibasis[0], chibasis_miss[0]
chipbasis0, chipbasis1 = np.dot(vec, chibasis), np.dot(vec, chibasis_miss)
tripleouterupdate(vec_bar, bv, chipbasis0, P)
tripleouterupdate(vec_bar, bv, chipbasis1, Pc)
In [32]:
tauscale = 12
eye = tauscale*np.pad(np.eye(Nsite0, dtype=int), ((0,0), (0,Nsite-Nsite0)), 'constant')
onevec = np.array([1,] + [0,]*(Nsite-1))
# We don't expect to need c^N+1 or c^N+2 so we ignore those...
# Matrices: <sum_c' W_cc' chi chi'> and higher order (GG, and WGG terms...)
Wbar = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
WGbar = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
WGGbar = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
# far-field versions of the same; the matched versions, followed by the "summed" (baseline) version:
Wbar_far = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
WGbar_far = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
WGGbar_far = np.zeros((Nbasis, Nbasis, N+1), dtype=int)
Wbar_far0 = np.zeros(N+1, dtype=int)
WGbar_far0 = np.zeros(N+1, dtype=int)
WGGbar_far0 = np.zeros(N+1, dtype=int)
# bias vector versions, including products with gamma:
biasvecbar = np.zeros((2, Nbasis, N+1), dtype=int)
biasGvecbar = np.zeros((2, Nbasis, N+1), dtype=int)
biasGGvecbar = np.zeros((2, Nbasis, N+1), dtype=int)
# residual bias vector versions:
resbiasave = np.zeros(N+1, dtype=int)
resbiasGave = np.zeros(N+1, dtype=int)
Pc, Pcc = np.zeros(N+1, dtype=int), np.zeros(N+1, dtype=int)
for st, P in tqdm_notebook(zip(symmstates, symmPstates), total=symmstates.shape[0],
disable=not __TESTING__):
Pc[1:] = P[:-1]
Pcc[2:] = P[:-2]
# basis0: those inside \chi, basis1: those outside \chi
chibasis = np.dot(chibasisfound_mat, st)
# chibasis0, chibasis1 = np.dot(chibasisfound_mat[0], st), chibasismiss_mat[0]
# construct our transition matrix:
W = np.zeros((Nsite0, Nsite), dtype=int)
for n, jumps in enumerate(jumplist):
if st[n] == 1: continue
for m in jumps:
if st[m] == 0:
W[n,n] -= 1
W[n,m] = 1
tau = -np.diag(W) # will be tau multiplied by tauscale = 12 (== -12//W[n,n])
Gam = W.copy() # Gamma matrix multiplied by tauscale = 12.
for n in range(Nsite0):
if tau[n] > 0:
tau[n] = tauscale//tau[n]
Gam[n,n] = 0
Gam[n] *= tau[n]
WG = -W[0,0]*np.dot(Gam[0,:Nsite0], Gam)+tauscale*tauscale*W[0,0]*onevec
WGG = np.dot(W[0,:Nsite0], np.dot(Gam[:,:Nsite0], Gam - 2*eye))
matrixupdate(Wbar, W[0], chibasis, chibasismiss_mat, chibasismissmatch_mat,
P, Pc, Pcc)
matrixupdate(WGbar, WG, chibasis, chibasismiss_mat, chibasismissmatch_mat,
P, Pc, Pcc)
matrixupdate(WGGbar, WGG, chibasis, chibasismiss_mat, chibasismissmatch_mat,
P, Pc, Pcc)
# far-field contributions of same:
farmatrixupdate(Wbar_far, W[0], chibasisfar_mat, Pc, Pcc)
farmatrixupdate(WGbar_far, WG, chibasisfar_mat, Pc, Pcc)
farmatrixupdate(WGGbar_far, WGG, chibasisfar_mat, Pc, Pcc)
Wbar_far0 += np.sum(W[0])*Pcc
WGbar_far0 += np.sum(WG)*Pcc
WGGbar_far0 += np.sum(WGG)*Pcc
# bias contributions (only bother if there's non-zero bias)
if tau[0]==0: continue
bv = np.sum(dx for j, dx in enumerate(dxlist) if st[j+1] == 0)
vectorupdate(biasvecbar, bv, onevec, chibasis, chibasismiss_mat, P, Pc)
vectorupdate(biasGvecbar, bv, Gam[0], chibasis, chibasismiss_mat, P, Pc)
vectorupdate(biasGGvecbar, bv, np.dot(Gam[0,:Nsite0],Gam-2*eye),
chibasis, chibasismiss_mat, P, Pc)
resbiasave += P*(bv[0]*bv[0]+bv[1]*bv[1])*tau[0]
bb = 0
for j, G in enumerate(Gam[0]):
if G>0:
bvp = np.array([0,0])
for k, dx in zip(jumplist[j], dxlist):
if st[k] == 0: bvp += dx
bb += G*np.dot(bv, bvp)*tau[j]
resbiasGave += P*bb
In [34]:
if __TESTING__:
print(Wbar_far0, WGbar_far0, WGGbar_far0)
Out[34]:
In [35]:
# scaling and symmetrization
symmWbar = symmetrize(Wbar, basisgroupops, basisgroupops)
symmWGbar = symmetrize(WGbar, basisgroupops, basisgroupops)/(tauscale*tauscale)
symmWGGbar = symmetrize(WGGbar, basisgroupops, basisgroupops)/(tauscale*tauscale)
symmWbar_far = symmetrize(Wbar_far, basisgroupops, basisgroupops)
symmWGbar_far = symmetrize(WGbar_far, basisgroupops, basisgroupops)/(tauscale*tauscale)
symmWGGbar_far = symmetrize(WGGbar_far, basisgroupops, basisgroupops)/(tauscale*tauscale)
symmresbiasave = resbiasave/tauscale
symmresbiasGave = resbiasGave/(tauscale*tauscale)
symmbiasvecbar = symmetrize(biasvecbar, groupops, basisgroupops)
symmbiasGvecbar = symmetrize(biasGvecbar, groupops, basisgroupops)/tauscale
symmbiasGGvecbar = symmetrize(biasGGvecbar, groupops, basisgroupops)/(tauscale*tauscale)
In [36]:
symmresbiasave
Out[36]:
In [37]:
symmresbiasGave
Out[37]:
In [39]:
def truncate_vec(v):
"""Return a vector that's shortened by truncating the high-order 0 components"""
return v[:(np.max(np.nonzero(v))+1)]
In [40]:
def printvecbasis(VB):
"""Print out the components of a vector-basis matrix"""
for d in range(2):
print("dim {}".format(d+1))
for bv, v in zip(basisfunc, VB[d]):
if np.any(v) != 0:
print(bv, truncate_vec(v))
In [41]:
def printbasisbasis(BB, comp=None):
"""Print out the components of a basis-basis matrix"""
for bv0, BB0 in zip(basisfunc, BB):
if comp is not None and tuple(bv0) not in comp: continue
for bv1, B in zip(basisfunc, BB0):
if np.any(B) != 0:
print(bv0, bv1, truncate_vec(B))
In [42]:
printbasisbasis(symmWbar_far, {(0,0)})
In [43]:
printbasisbasis(symmWbar-symmWbar_far)
In [44]:
printbasisbasis(symmWGbar_far, {(0,0)})
In [45]:
printbasisbasis(symmWGbar-symmWGbar_far)
In [46]:
printbasisbasis(symmWGGbar_far, {(0,0)})
In [47]:
printbasisbasis(symmWGGbar-symmWGGbar_far)
In [48]:
printvecbasis(symmbiasvecbar)
In [49]:
printvecbasis(symmbiasGvecbar)
In [50]:
printvecbasis(symmbiasGGvecbar)
In [62]:
import h5py
rewriteFile = False
printFile = False
In [63]:
if rewriteFile:
with h5py.File('Neighbor-averaging.hdf5', 'w') as f:
f['dxlist'] = np.array(dxlist)
f['sites'] = np.array(sites)
f['jumplist'] = np.array(jumplist)
f['basisfunc'] = np.array(basisfunc)
f['symmWbar'] = symmWbar
f['symmWGbar'] = symmWGbar
f['symmWGGbar'] = symmWGGbar
f['symmWbar_far'] = symmWbar_far
f['symmWGbar_far'] = symmWGbar_far
f['symmWGGbar_far'] = symmWGGbar_far
f['symmresbias'] = symmresbiasave
f['symmresbiasGave'] = symmresbiasGave
f['symmbiasvecbar'] = symmbiasvecbar
f['symmbiasGvecbar'] = symmbiasGvecbar
f['symmbiasGGvecbar'] = symmbiasGGvecbar
In [64]:
if printFile:
with h5py.File('Neighbor-averaging.hdf5', 'r') as f:
for k, c in f.items():
print(k)
print(c.value)
In [65]:
def mpmesh(Ndiv, pre=np.pi):
"""
Generates a MP mesh for a square lattice.
:param Ndiv: number of divisions
:param pre: prefactor for edge of Brilloiun zone (pi/a_0)
:returns k[Nk,2]: k-points
:returns w[Nk]: weight
"""
prescale = pre/Ndiv
wscale = 1./(Ndiv*Ndiv)
Nk = (Ndiv*(Ndiv+1))//2
kpt, w = np.zeros((Nk,2)), np.zeros(Nk)
i = 0
for n in range(Ndiv):
for m in range(n+1):
kpt[i,0] = prescale*(n+0.5)
kpt[i,1] = prescale*(m+0.5)
if n==m:
w[i] = wscale
else:
w[i] = 2*wscale
i += 1
return kpt, w
In [68]:
square = crystal.Crystal(np.eye(2), [np.zeros(2)])
chem = 0
sitelist = square.sitelist(chem)
jumpnetwork = square.jumpnetwork(chem, 1.01) # [[((0,0), dx) for dx in dxlist]]
In [69]:
starset = stars.StarSet(jumpnetwork, square, chem, 3)
vecstarset = stars.VectorStarSet(starset)
if __TESTING__:
print(starset)
In [70]:
if __TESTING__:
for vR, vV in zip(vecstarset.vecpos, vecstarset.vecvec):
print('')
for R, v in zip(vR, vV):
print(starset.states[R] , v)
In [71]:
GF = GFcalc.GFCrystalcalc(square, chem, sitelist, jumpnetwork, kptwt = mpmesh(32))
GF.SetRates(np.ones(1), np.zeros(1), np.ones(1), np.zeros(1))
In [72]:
if __TESTING__:
print(GF)
In [73]:
GFmat, GFstarset = vecstarset.GFexpansion()
GF0array = np.array([GF(0,0,GFstarset.states[s[0]].R) for s in GFstarset.stars])
g0 = np.dot(GFmat, GF0array)
print(g0)
In [76]:
basis2state = [starset.stateindex(stars.PairState(0, 0, bv, bv)) for bv in basisfunc]
basis2star = [starset.starindex(stars.PairState(0, 0, bv, bv)) for bv in basisfunc]
In [79]:
if __TESTING__:
for bv, stateind, starind in zip(basisfunc, basis2state, basis2star):
print(bv, stateind, starind)
In [80]:
state2basis = [basis2state.index(n) for n in range(starset.Nstates)]
In [81]:
if __TESTING__:
print(state2basis)
Now the real conversion begins! We start by mapping all of the bias vectors and local functions onto our vectorBasis.
In [107]:
NVS = vecstarset.Nvstars
symmbiasvecVS = np.zeros((N+1, NVS))
symmbiasGvecVS = np.zeros((N+1, NVS))
symmbiasGGvecVS = np.zeros((N+1, NVS))
for i in range(vecstarset.Nvstars):
for Ri, vi in zip(vecstarset.vecpos[i], vecstarset.vecvec[i]):
bi = state2basis[Ri]
symmbiasvecVS[:, i] += np.dot(vi, symmbiasvecbar[:,bi,:])
symmbiasGvecVS[:, i] += np.dot(vi, symmbiasGvecbar[:,bi,:])
symmbiasGGvecVS[:, i] += np.dot(vi, symmbiasGGvecbar[:,bi,:])
stars.zeroclean(symmbiasvecVS);
stars.zeroclean(symmbiasGvecVS);
stars.zeroclean(symmbiasGGvecVS);
In [100]:
for nv in range(NVS):
if not np.allclose(symmbiasvecVS[:,nv], 0):
print(nv, truncate_vec(symmbiasvecVS[:,nv]))
In [101]:
for nv in range(NVS):
if not np.allclose(symmbiasGvecVS[:,nv], 0):
print(nv, truncate_vec(symmbiasGvecVS[:,nv]))
In [102]:
for nv in range(NVS):
if not np.allclose(symmbiasGGvecVS[:,nv], 0):
print(nv, truncate_vec(symmbiasGGvecVS[:,nv]))
In [108]:
symmWbarVS = np.zeros((N+1, NVS, NVS))
symmWGbarVS = np.zeros((N+1, NVS, NVS))
symmWGGbarVS = np.zeros((N+1, NVS, NVS))
for i in range(vecstarset.Nvstars):
for Ri, vi in zip(vecstarset.vecpos[i], vecstarset.vecvec[i]):
bi = state2basis[Ri]
for j in range(vecstarset.Nvstars):
for Rj, vj in zip(vecstarset.vecpos[j], vecstarset.vecvec[j]):
bj = state2basis[Rj]
vivj = np.dot(vi,vj)
symmWbarVS[:, i, j] += vivj*(symmWbar[bi,bj,:]-symmWbar_far[bi,bj,:])
symmWGbarVS[:, i, j] += vivj*(symmWGbar[bi,bj,:]-symmWGbar_far[bi,bj,:])
symmWGGbarVS[:, i, j] += vivj*(symmWGGbar[bi,bj,:]-symmWGGbar_far[bi,bj,:])
stars.zeroclean(symmWbarVS);
stars.zeroclean(symmWGbarVS);
stars.zeroclean(symmWGGbarVS);
In [104]:
for nv,mv in itertools.product(range(NVS), repeat=2):
if not np.allclose(symmWbarVS[:,nv,mv], 0):
print(nv, mv, truncate_vec(symmWbarVS[:,nv,mv]))
In [105]:
for nv,mv in itertools.product(range(NVS), repeat=2):
if not np.allclose(symmWGbarVS[:,nv,mv], 0):
print(nv, mv, truncate_vec(symmWGbarVS[:,nv,mv]))
In [106]:
for nv,mv in itertools.product(range(NVS), repeat=2):
if not np.allclose(symmWGGbarVS[:,nv,mv], 0):
print(nv, mv, truncate_vec(symmWGGbarVS[:,nv,mv]))
Our "far" functions represent the translationally invariant contributions, and this requires Fourier transforms, and Taylor expansions to then be made into local contributions.
Mathematically, we're attempting to compute $\eta_i\cdot M_{ij}\cdot\eta_j$; the issue is that $\eta_i$ does not go to zero in the far-field (it's not local), and $M$ can be written as a local function plus a translationally invariant function $M^0$. Only the latter is problematic. However, as $\eta_i$ comes from a Green function solution (using the Dyson equation), if we multiply by the $w^0$, we produce a local function. Hence, we can rewrite that matrix equation as $(w^0\eta)_i\cdot (g^0M^0g^0)_{ij}\cdot (w^0\eta_j)$. Now, then we "simply" need to evaluate $g^0M^0g^0$, which can be done using Fourier transforms, as it is the product of three translationally invariant functions.
In [117]:
def FT(mat, kptwt):
"""
(real) Fourier transform of translationally invariant function.
:param mat[Nbasis, N+1]: far-field version of matrix;
each Nbasis is relative to 0
:param kptwt: tuple of (kpt[Nkpt, 2], wt[Nkpt])
:returns matFT[Nkpt, N+1]: FT of matrix
"""
kpt = kptwt[0]
matFT = np.zeros((kpt.shape[0], N+1))
for bv, matv in zip(basisfunc, mat):
matFT += np.outer(np.cos(np.dot(kpt, bv)), matv)
return matFT
In [146]:
PE.Taylor2D(Lmax=6); # initialize
def Taylor(mat):
"""
(real) Taylor expansion of Fourier transform of translationally invariant function.
:param mat[Nbasis, N+1]: far-field version of matrix;
each Nbasis is relative to 0
:returns matTaylor: T2D version of FT Taylor expansion matrix
"""
pre = np.array([1., 0., -1/2, 0., 1/24]) # Taylor coefficients for cos()
matTaylor = PE.Taylor2D()
for bv, matv in zip(basisfunc, mat):
for ve in PE.Taylor2D.constructexpansion([(matv, bv)], pre=pre):
matTaylor += ve
matTaylor.reduce()
return matTaylor
In [118]:
if __TESTING__:
print(FT(symmWbar_far[0], mpmesh(4)))
Out[118]:
In [149]:
g0Taylor = (Taylor(symmWbar_far[0])[1]).inv() # extract out the "constant" term
In [150]:
print(g0Taylor)
In [167]:
g0WGbarTaylor = ( (g0Taylor*g0Taylor)*Taylor(symmWGbar_far[0])).reduce().truncate(0)
g0WGGbarTaylor = ( (g0Taylor*g0Taylor)*Taylor(symmWGGbar_far[0])).reduce().truncate(0)
In [168]:
print(g0WGbarTaylor)
In [169]:
print(g0WGGbarTaylor)
In [170]:
kpt, wt = mpmesh(32)
In [189]:
g0FT = 1./FT(symmWbar_far[0], (kpt, wt))[:,1]
WGbarFT = FT(symmWGbar_far[0], (kpt, wt))
WGGbarFT = FT(symmWGGbar_far[0], (kpt, wt))
In [179]:
if __TESTING__:
print(g0FT)
Out[179]:
In [183]:
pmax = np.sqrt(min([np.dot(G, G) for G in square.BZG]) / -np.log(1e-11))
prefactor = square.volume
g0Taylor_fnlp = {(n, l): GFcalc.Fnl_p(n, pmax) for (n, l) in g0Taylor.nl()}
g0Taylor_fnlu = {(n, l): GFcalc.Fnl_u(n, l, pmax, prefactor, d=2)
for (n, l) in g0Taylor.nl()}
if __TESTING__:
print(pmax)
Out[183]:
In [187]:
if __TESTING__:
print(g0Taylor.nl(), g0WGbarTaylor.nl(), g0WGGbarTaylor.nl())
Out[187]:
In [191]:
g0WGbarsc = np.zeros_like(g0WGbarFT)
g0WGGbarsc = np.zeros_like(g0WGGbarFT)
for i, k in enumerate(kpt):
g0WGbarsc[i] = (g0FT[i]**2)*g0WGbarFT[i] - g0WGbarTaylor(k, g0Taylor_fnlp).real
g0WGGbarsc[i] = (g0FT[i]**2)*g0WGGbarFT[i] - g0WGGbarTaylor(k, g0Taylor_fnlp).real
In [197]:
if __TESTING__:
print(truncate_vec(np.dot(wt, g0WGGbarsc)))
Out[197]:
In [227]:
# this list is a bit of overkill, but...
veclist = [GFstarset.states[s[0]].dx for s in GFstarset.stars]
g0WGbar, g0WGGbar = [], []
for x in veclist:
coskx = np.sum(np.cos(np.tensordot(kpt, np.dot(g, x), axes=(1, 0)))
for g in groupops) / 8
g0WGbar.append(np.dot(wt*coskx,g0WGbarsc) + g0WGbarTaylor(x, g0Taylor_fnlu).real)
g0WGGbar.append(np.dot(wt*coskx,g0WGGbarsc) + g0WGGbarTaylor(x, g0Taylor_fnlu).real)
In [230]:
for v, g in zip(veclist, g0WGbar):
print(v, truncate_vec(g))
In [231]:
for v, g in zip(veclist, g0WGGbar):
print(v, truncate_vec(g))
The Green function, and the correction $\eta$, end up having particularly simple expressions, that we will compute directly (it requires some simplification of the polynomial expressions which are more difficult to directly express here. It, unfortunately, also introduces a denominator polynomial which makes some of our expressions more complicated.
We have $$\eta_i = -2\frac{g^0_{i0}}{1+g^0_{i0} - (1+3g^0_{i0})c_\text{B}}$$
In [232]:
@njit(nogil=True, parallel=True)
def polymult(p, q):
"""
Multiplication of two polynomial coefficients, where
p(x) = sum_n p[n] * x^n
:param p: polynomial coefficients for p
:param q: polynomial coefficients for q
:returns pq: polynomial coefficients for pq
"""
P = p.shape[0]-1
Q = q.shape[0]-1
pq = np.zeros(P+Q+1)
for n in range(P+Q+1):
for i in range(max(0,n-Q), min(n,P)+1):
pq[n] += p[i]*q[n-i]
return pq
In [252]:
@njit(nogil=True, parallel=True)
def polydiv(p, a):
"""
Division of polynomial p(x) by (x-a)
:param p: polynomial coefficients for p
:param a: term in nomial (x-a)
:returns d, r: divisor d(x), and remainder r
"""
P = p.shape[0]-1
d = np.zeros(P)
d[P-1] = p[P]
for n in range(P-2,-1,-1):
d[n] = p[n+1] + a*d[n+1]
return d, p[0] + a*d[0]
In [278]:
divpoly = np.zeros(N+1)
divpoly[0], divpoly[1] = 1+g0[0,0], -(1+3*g0[0,0])
etabar_div = -2*g0[0] # this is etabar*div, so that etabar = etabar_div/div
etaW0_div = np.zeros(N+1)
etaW0_div[0] = -2 # this is W0*etabar*div (for the translational invariant terms)
In [268]:
# unbiased:
L0 = np.zeros(N+1)
L0[0], L0[1] = 1., -1.
In [354]:
# Note: vecstarset.outer[i,j, v1, v2] = 1/2 delta_ij delta_v1v2,
# so we can use dot-products throughout
# SCGF:
L1 = 0.5*np.dot(symmbiasvecVS, etabar_div)
L_SCGF = polymult(L0, divpoly)[:N+1] + L1
polydiv(L_SCGF, 1)
Out[354]:
In [355]:
# print(np.dot(GFmat[0,0], g0WGGbar))
PsiB = polymult(polymult(divpoly, divpoly), symmresbiasave)[:N+1] + \
-2*polymult(divpoly, np.dot(symmbiasGvecVS, etabar_div))[:N+1] + \
np.dot(np.dot(symmWGbarVS, etabar_div), etabar_div) + \
4*np.dot(GFmat[0,0], g0WGbar) # far-field; note: etaW0_div == 2, so factor of 4
print(PsiB)
In [360]:
WR = polymult(polymult(divpoly, divpoly), symmresbiasGave)[:N+1] - \
polymult(polymult(divpoly, divpoly), symmresbiasave)[:N+1] + \
-2*polymult(divpoly, L1)[:N+1] + \
-2*polymult(divpoly, np.dot(symmbiasGGvecVS, etabar_div))[:N+1] + \
np.dot(np.dot(symmWGGbarVS, etabar_div), etabar_div) + \
4*np.dot(GFmat[0,0], g0WGGbar)
print(WR)
In [366]:
# Now, to put it together, and do the division...
cBv = np.linspace(0.01,1,num=99,endpoint=False)
D1, D2 = [], []
for cB in cBv:
# print(cB)
cA = 1-cB
cpow = np.array([cB**n for n in range(N+1)])
L0c, divc, L1c = np.dot(cpow, L0), np.dot(cpow, divpoly), np.dot(cpow, L_SCGF)
L1c /= divc
PsiBc, WRc = np.dot(cpow, PsiB)/(divc*divc), np.dot(cpow, WR)/(divc*divc)
L2c = L1c + 0.5*PsiBc*PsiBc/WRc
D0c = L0c/cA
D1c = L1c/cA
D2c = L2c/cA
D1.append(D1c)
D2.append(D2c)
print(cB, D1c, D2c, D2c/D1c) #, PsiBc)
D1v, D2v = np.array(D1), np.array(D2)
In [376]:
plt.rcParams['figure.figsize'] = (8,8)
fig, ax = plt.subplots()
ax.plot(cBv, D1, 'b', label='GF')
ax.plot(cBv, D2, 'r', label='GF+resbias')
ax.set_ylabel('$D^{\\rm A}$', fontsize='x-large')
ax.set_xlabel('$c_{\\rm B}$', fontsize='x-large')
ax.legend(bbox_to_anchor=(0.5,0.5,0.5,0.3), ncol=1, shadow=True,
frameon=True, fontsize='x-large', framealpha=1.)
plt.tight_layout()
plt.show()
In [388]:
num_SCGF, denom_SCGF = truncate_vec(-polydiv(L_SCGF,1)[0]), truncate_vec(divpoly)
In [389]:
num_SCGFbc, denom_SCGFbc = \
truncate_vec(-polydiv(0.5*polymult(PsiB,PsiB),1)[0]), \
truncate_vec(polymult(polymult(divpoly, divpoly), WR))
In [391]:
# check remainders (should be 0 for both)
if __TESTING__:
print(polydiv(L_SCGF,1)[1], polydiv(0.5*polymult(PsiB,PsiB),1)[1])
Out[391]:
In [411]:
def print_fraction(numer, denom, powstring='**'):
"""
Returns a string representation of our polynomial ratio
"""
def format_pow(n):
if n==0:
return ''
if n==1:
return '*c'
return '*c' + powstring +'{}'.format(n)
# first, "divide" through until lowest order is constant on both:
while np.isclose(numer[0], 0) and np.isclose(denom[0], 0):
numer, denom = numer[1:], denom[1:]
# second, scale everything by lowest order term in denominator
scale = denom[np.min(np.nonzero(denom))]
numer /= scale
denom /= scale
s = '('
for n, coeff in enumerate(numer):
if not np.isclose(coeff, 0):
s += '{:+.10g}'.format(coeff) + format_pow(n)
s += ')/('
for n, coeff in enumerate(denom):
if not np.isclose(coeff, 0):
s += '{:+.10g}'.format(coeff) + format_pow(n)
s += ')'
return s
In [412]:
print(print_fraction(num_SCGF, denom_SCGF))
In [415]:
print(print_fraction(num_SCGF, denom_SCGF) + ' + ' +\
print_fraction(num_SCGFbc, denom_SCGFbc))
Note: both of these polynomials have two factors of $(1-c)$ in them; so we can simplify further...
In [421]:
polydiv(polydiv(polydiv(num_SCGFbc,1)[0],1)[0],1)
Out[421]:
In [423]:
polydiv(polydiv(denom_SCGFbc,1)[0],1)
Out[423]:
In [424]:
SCGFbc_func = print_fraction(num_SCGF, denom_SCGF) + ' + ' +\
print_fraction(polydiv(polydiv(num_SCGFbc,1)[0],1)[0],
polydiv(polydiv(denom_SCGFbc,1)[0],1)[0])
print(SCGFbc_func)