In [1]:
####Implementing the one line poisson solver.
"""This notebook has a code written in Python2.
The code solves the poisson's equation numerically with in a unit cell
using the method described in Tomas Arias lecture.""";
In [3]:
"""Importing all the necessary funtions required for this code"""
# numpy for math calculations
import numpy as np
# For pplotting or visualising
import matplotlib.pylab as pylab
%pylab notebook
pylab.rcParams['figure.figsize'] = (10, 6)
import matplotlib.pyplot as plt
# For performing linear algebra operations using matrices
from numpy import linalg as LA
from numpy.linalg import inv
import operator
# for Generating tuples list.
import itertools
In [28]:
### Creating the setup.m file. All these inputs are setup.
"""Creating the Cell and Geometry
Defining the variables
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
s=[a,b,c] => Defines the number of sample points along the
lattice vectors. (a,b,c) are any integers.
sig=[a,b,c,...] => Defining sigma for each gaussian funtion to be summed up.
coeff=array of {+/-}coefficients. The factor is used for gaussian.
Gaussian => coeff[0].Gaussian(sig[0])+coeff[1].Gaussian(sig[1])+....so on.
"""
#Inputs for Geometry
s = [20,25,30]
R= np.diag([6,6,6])
#Input for Charge distribution
sig=[0.75,0.5]
coeff=[-1,1]
In [5]:
def _M(s):
"""Defining the 'M' matrix.
'M' has coordinates in matrix form for each sample point in real space.
We use M to get list of real space sample points 'r'
r=M(Diag(S))^(-1)R^(Transpose)
Args:
input:
s => number of sample points along the lattice vectors.
output:
M matrix.
"""
M = tuple(itertools.product(*[range(0,s[0]),range(0,s[1]),range(0,s[2])]))
M=[list(i) for i in M]
return M
In [6]:
def _N(s):
"""Defining the 'N' matrix.
'N' has coordinates in matrix form for each sample point in reciprocal
space.
We use N to get list of reciprocal lattice vectos 'G'
G= 2*pi*N*R^(inverse)
Args:
input:
s => number of sample points along the lattice vectors.
output:
N matrix.
"""
M=_M(s)
m1=np.transpose(M)[0]
m2=np.transpose(M)[1]
m3=np.transpose(M)[2]
n1=[ m1[i]-s[0] if m1[i] > s[0]/2 else m1[i] for i in range(len(m1)) ]
n2=[ m2[i]-s[1] if m2[i] > s[1]/2 else m2[i] for i in range(len(m2)) ]
n3=[ m3[i]-s[2] if m3[i] > s[2]/2 else m3[i] for i in range(len(m3)) ]
N=np.transpose([n1,n2,n3])
return N
In [7]:
def _gen_r_G_G2(s,R,p=None):
"""Generates r, G and G2 matrices.
Args: input:
s => number of sample points along the lattice vectors.
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
p = 'r', 'G' or 'G2'(square of G matrix)
Output:
p='r' => Outputs r matrix.
p='G' => Outputs G matrix.
p='2' => Outputs G2 matrix.
Example: s = [15,15,15]
R= np.diag([6,6,6])
_gen_r_G_G2(s,R,"r")
"""
if p is None:
return "Please pass a third argument asking for 'r', 'G' or 'G2'."
if p is "r":
r=np.dot(np.dot(_M(s),inv(np.diag(s))),np.transpose(R))
return r
if p is "G":
G=2*np.pi*np.dot(_N(s),inv(R))
return G
if p is "G2":
G=2*np.pi*np.dot(_N(s),inv(R))
G2=np.sum(np.transpose(G*G),axis=0)
return G2
else:
return "Error: Passed an argument other than r, G or G2."
In [8]:
def gaussian(sig,coeff,x):
""" Defining Gaussian charge distribution at point x in 3D space.
Args:
inputs:
sig => sigmas for each Gaussian.
coeff => coefficients for each gaussian funtion.
Note: All the gaussians will be summed up.
x (np.array()): List of input points
Note: x in 3D space is x=Sqrt(a^2+b^2+c^2) where
{a,b,c} are the cartesian coordinates.
Output: Charge distribution at points x.
"""
g=lambda x,sigma: (np.exp(-(x**2)/(2*sigma**2)))/(2*np.pi*sigma**2)**(3./2)
n=lambda x: np.sum([coeff[i]*g(x,sig[i]) for i in range(len(sig))],axis=0)
if len(sig)!=len(coeff):
return "Error: Length of 'sig' not equal to length of 'coeff'."
else:
return n(x)
In [9]:
def charge_dis(s,R,sig,coeff):
"""Generates the charge distribution at each sample point for the unit cell.
We change the midpoint of the unitcell to origin here.
Args:
input:
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
s=[a,b,c] => Defines the number of sample points along the
lattice vectors. (a,b,c) are any integers.
"""
mid = sum(R,axis=1)/2.
dr=np.sqrt(sum((_gen_r_G_G2(s,R,"r")-mid)**2,axis=1))
return gaussian(sig,coeff,dr)
In [10]:
def Ope(inp,R):
"""Outputs the overlap matrix. In our case here we are specific to plane wave basis.
Args:
input:
inp (numpy array): The overlap matrix operator acts on this input.
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
"""
return np.linalg.det(R)*np.dot(np.identity(len(inp)),inp)
In [11]:
def Lope(inp,s,R):
"""Lope is the L operator as defined in Tomas arias Lecture.
Note:Works only for plane wave basis.
Args:
input:
inp (numpy array): The overlap matrix operator acts on this input.
s=[a,b,c] => Defines the number of sample points along the
lattice vectors. (a,b,c) are any integers.
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
"""
if len(_gen_r_G_G2(s,R,"G2")) == len(inp):
return -np.linalg.det(R)*np.dot(np.diag(_gen_r_G_G2(s,R,"G2")),inp)
else:
return "Error: Input vector size is not same as G2"
In [12]:
def LinvOpe(inp,s,R):
"""LinvOpe is the inverse of L operator as defined in Tomas arias Lecture.
Note:Works only for plane wave basis.
Args:
input:
inp (numpy array): The overlap matrix operator acts on this input.
s=[a,b,c] => Defines the number of sample points along the
lattice vectors. (a,b,c) are any integers.
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
"""
if len(_gen_r_G_G2(s,R,"G2")) == len(inp):
Lele=-np.linalg.det(R)*np.diag(_gen_r_G_G2(s,R,"G2"))
Lele[0][0]=1
Leleinv=inv(Lele)
Leleinv[0][0]=0
return np.dot(Leleinv,inp)
else:
return "Error: Input vector size is not same as G2"
In [13]:
def cI(inp,s,R,p=None):
"""cI is the fourier series expansion coeff Matrix. It is the forward transform matrix for the problem.
Note:Works only for plane wave basis.
Args:
input:
inp (numpy array): The overlap matrix operator acts on this input.
s=[a,b,c] => Defines the number of sample points along the
lattice vectors. (a,b,c) are any integers.
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
p (default): Fast fourier transform. #Written with Wiley's code help.
p='s' => slow fourier transform.
"""
if p is None:
f=np.fft.fftn(inp.reshape(s,order="F")).reshape(np.prod(s),order="F")
elif p is "s":
f=np.dot(np.exp(1j*np.dot(_gen_r_G_G2(s,R,"G"),np.transpose(_gen_r_G_G2(s,R,"r")))),inp)
return f
In [14]:
def cJ(inp,s,R,p=None):
"""cI is the fourier series expansion coeff Matrix. You can get the fast fourier transform too.
Note:Works only for plane wave basis.
Args:
input:
inp (numpy array): The overlap matrix operator acts on this input.
s=[a,b,c] => Defines the number of sample points along the
lattice vectors. (a,b,c) are any integers.
R=[a,b,c] => Defines the lattice vectors. (a,b,c) are any +ve integers.
"""
if p is None:
f=np.fft.ifftn(inp.reshape(s,order="F")).reshape(np.prod(s),order="F")
return f
elif p is "s":
f=np.exp(-1j*np.dot(_gen_r_G_G2(s,R,"G"),np.transpose(_gen_r_G_G2(s,R,"r"))))
return (1./np.prod(s))*np.dot(np.transpose(f.conjugate()),inp)
In [15]:
def phi(inp,s,R):
"""Poisson solver in a single line of code."""
return cI(LinvOpe(-4.0*np.pi*Ope(cJ(inp,s,R),R),s,R),s,R)
In [16]:
Uanalytic=((1/sig[0] + 1/sig[1])/2 -np.sqrt(2)/np.sqrt(sig[0]**2+sig[1]**2))/np.sqrt(np.pi)
Uanalytic
Out[16]:
In [17]:
Unumerical=0.5*np.real(np.dot(np.conjugate(cJ(np.real(phi(charge_dis(s,R,sig,coeff),s,R)),s,R,"s")),Ope(cJ(charge_dis(s,R,sig,coeff),s,R,"s"),R)))
Unumerical
Out[17]:
In [17]:
#################code ends ;
In [18]:
# Just to check how long the code takes to run.
import time
start_time = time.time()
Unumerical=0.5*np.real(np.dot(np.conjugate(cJ(np.real(phi(charge_dis(s,R,sig,coeff),s,R)),s,R)),Ope(cJ(charge_dis(s,R,sig,coeff),s,R),R)))
print("--- %s seconds ---" % (time.time() - start_time))
In [19]:
###Ewald energy Calculator;
In [29]:
## Additions to the setup.m file.
"""Place to put the location of the ionic cores.
stored in a nx3 Matrix. The charges are called Z.
"""
# Defining atomic locations and nuclear charge.
"""Two protons at a distance of 1.75 Bohr. H2 molecule."""
X=np.array([[0,0,0],[1.75,0,0]])
Z=1
In [30]:
# Defining the structure factor.
#Note:In paper he said sum(#,2) => in matlab equivalet to python axis=1
def _sf(s,R,X):
return np.sum(np.exp(1j*np.dot(_gen_r_G_G2(s,R,"G"),np.transpose(X))),axis=1)
In [31]:
# Changing the width of Gaussian. We have only one Gaussian with width 0.25.
sig=[0.25]
coeff=[1]
In [32]:
# As the norm of Gaussian should be Z instead of unity. I'm redifining the funtion.
def gaussian(sig,coeff,Z,s,R,X):
""" Defining Gaussian charge distribution at point x in 3D space.
Args:
inputs:
sig => sigmas for each Gaussian.
coeff => coefficients for each gaussian funtion.
Note: All the gaussians will be summed up.
x (np.array()): List of input points
Note: x in 3D space is x=Sqrt(a^2+b^2+c^2) where
{a,b,c} are the cartesian coordinates.
z => Charge of ions.
#Changes: Added new input Z. So the norm is now Z.
Evaluation of total charge.
Output: Charge distribution at points x.
"""
g=lambda x,sigma: Z*(np.exp(-(x**2)/(2*sigma**2)))/(2*np.pi*sigma**2)**(1./2) # Changed 3./2 to 1./2
mid = sum(R,axis=1)/2.
dr=np.sqrt(sum((_gen_r_G_G2(s,R,"r")-mid)**2,axis=1))
g1=g(np.array([dr[i] for i in range(len(dr))]),np.array(sig))
n=np.real(cI(np.prod([cJ(g1,s,R),_sf(s,R,X)],axis=0),s,R))
if len(sig)!=len(coeff):
return "Error: Length of 'sig' not equal to length of 'coeff'."
else:
return n
In [33]:
Y=X+[3,3,3]
Y
Out[33]:
In [34]:
k=gaussian(sig,coeff,Z,s,R,Y)
t=gaussian(sig,coeff,Z,s,R,X)
#k.real[abs(k.real) < 10e-6] = 0.0
In [35]:
Uself=Z**2/(2*np.sqrt(np.pi))*(1/sig[0])*len(X)
Uself
Out[35]:
In [25]:
Unumerical=0.5*np.real(np.dot(np.conjugate(cJ(np.real(phi(t,s,R)),s,R)),Ope(cJ(t,s,R),R)))
Unumerical
Out[25]:
In [36]:
Unumerical=0.5*np.real(np.dot(np.conjugate(cJ(np.real(phi(k,s,R)),s,R)),Ope(cJ(k,s,R),R)))
Unumerical
Out[36]:
In [27]:
Unumerical-Uself
Out[27]:
In [4]:
# Working on the Survey paper.
from math import erf, erfc, sqrt
import numpy as np
from itertools import product
from math import pi
from cmath import exp
In [8]:
def U_r(qi,ri,R,a,n):
"""
Args: qi => list of all charges.
ri => position vectors
R => vell vectors to find volume
a => alpha parameter.
n => number of terms in sum.
Output: Gives the U^r term in the Survey paper.
"""
result = 0
neighbs = np.array([np.array(i) for i in list(product(range(-n,n+1),repeat=3))])
neighbs = np.dot(R,neighbs.T).T
count = 0
for n_L in neighbs:
for a_i in range(len(qi)):
for a_j in range(len(qi)):
if np.all(n_L == 0) and (a_i != a_j):
count += 1
rijn = np.linalg.norm(ri[a_i]-ri[a_j])
result += qi[a_i]*qi[a_j]*erfc(a *rijn)/rijn
elif np.any(n_L != 0):
count += 1
rijn = np.linalg.norm(rs[a_i]-rs[a_j]+n_L)
result += qi[a_i]*qi[a_j]*erfc(a *rijn)/rijn
return result/2.
In [9]:
def U_m(qi,ri,R,a,m):
"""
Args: qi => list of all charges.
ri => position vectors
R => vell vectors to find volume
a => alpha parameter.
m => reciprocal space vector
Output: Gives the U^m term in the Survey paper.
"""
results = 0
V = np.dot(R[0],np.cross(R[1],R[2]))
k1 = 2*np.pi*np.cross(R[1],R[2])/V
k2 = 2*np.pi*np.cross(R[2],R[0])/V
k3 = 2*np.pi*np.cross(R[0],R[1])/V
K = np.array([k1,k2,k3])
ms = [np.dot(K,np.array(i).T) for i in list(product(list(range(-int(m/2 +1),0)+list(range(int(m/2)))),repeat=3))]
print(len(ms))
for m in ms:
if np.any(m != 0):
for a_i in range(len(qi)):
for a_j in range(len(qi)):
results += qi[a_i]*qi[a_j]* exp(-(pi*pi * np.dot(m,m)/(a*a)) + 2*pi*1j* np.dot(m,ri[a_i]-ri[a_j]))/np.dot(m,m)
return results/(2*pi*V)
In [10]:
def U_o(qi,a):
"""
Args: qi => list of all charges.
a => alpha parameter.
Output: Gives the U^0 term in the Survey paper.
"""
result = 0
for z in qi:
result += z*z
return -a*result/(2*np.sqrt(np.pi))
In [12]:
def phi_survey(qi,ri,R,a,m,n):
return U_r(qi,ri,R,a,n)+U_m(qi,ri,R,a,m)+U_o(qi,a)
In [13]:
phi_survey([1,1],np.array([[0.0,0.0,0.0],[1.75,0.0,0.0]]),[[6,0,0],[0,6,0],[0,0,6]],0.55024769904855042,20,4)
In [ ]: