``````

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

``````
``````

Populating the interactive namespace from numpy and matplotlib

``````
``````

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]:

0.055142527694733337

``````
``````

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]:

0.055140128732007555

``````
``````

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))

``````
``````

--- 1.79489994049 seconds ---

``````
``````

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]:

array([[ 3.  ,  3.  ,  3.  ],
[ 4.75,  3.  ,  3.  ]])

``````
``````

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]:

2.2567583341910251

``````
``````

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]:

0.28656490534959411

``````
``````

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]:

0.017879942593415472

``````
``````

In [27]:

Unumerical-Uself

``````
``````

Out[27]:

-1.9701934288414309

``````
``````

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)

``````
``````

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-a89bee8b1f96> in <module>()
----> 1 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)

<ipython-input-12-1ff4441fa690> in phi_survey(qi, ri, R, a, m, n)
1 def phi_survey(qi,ri,R,a,m,n):
----> 2     return U_r(qi,ri,R,a,n)+U_m(qi,ri,R,a,m)+U_o(qi,a)

<ipython-input-8-5171beba4dd7> in U_r(qi, ri, R, a, n)
22                 elif np.any(n_L != 0):
23                     count += 1
---> 24                     rijn = np.linalg.norm(rs[a_i]-rs[a_j]+n_L)
25                     result += qi[a_i]*qi[a_j]*erfc(a *rijn)/rijn
26

NameError: global name 'rs' is not defined

``````
``````

In [ ]:

``````