In [5]:

    

%load_ext autoreload
%autoreload 2


    

In [1]:

    

import numpy as np
from numpy.linalg import det
from numpy.linalg import inv
from numpy.fft import fftn, ifftn
from itertools import product
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


    

In [2]:

    

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from pydft.crystal import *
from pydft.bases.planewave import *


    

In [32]:

    

sigma1 = 0.75
sigma2 = 0.50

lat_const = 6
r = [lat_const,0,0]
S = [16,16,16]
R = makeR(r)
origin = 1./2*np.sum(R,1)
crystal = Crystal(R,S)
g(crystal)

# Create vector of basis expansion coefficients for the electron density.
n = -gauss_charge_dist(crystal._r,origin,sigma1) + (
    gauss_charge_dist(crystal._r,origin,sigma2))


    

In [37]:

    

# Plane refers to one of the pages of the n and r matrices.
plane=0

npl = n[plane*S[1]*S[2]:(plane+1)*S[1]*S[2]]
rpl = crystal._r[plane*S[1]*S[2]:(plane+1)*S[1]*S[2]]
xpl = [rpl[i,0] for i in range(len(rpl))]
ypl = [rpl[i,1] for i in range(len(rpl))]

fig = plt.figure()
ax = plt.gca(projection='3d')
ax.scatter(xpl,ypl,list(npl))
plt.title("charge density distribution")
plt.show()
plt.close()


    

In [14]:

    

ϕ=-4*np.pi*I(crystal, Linv(crystal, O(crystal, J(crystal, n))))
UNumeric = 1./2*np.real(np.dot(J(crystal, ϕ),O(crystal, J(crystal, n))))
UNumeric


    

    




/Users/jeremy/Codes/dft/pydft/bases/planewave.py:74: RuntimeWarning: divide by zero encountered in double_scalars
  np.linalg.norm(gi)**-2 for gi in supercell._g])*v
/Users/jeremy/Codes/dft/pydft/bases/planewave.py:74: RuntimeWarning: invalid value encountered in multiply
  np.linalg.norm(gi)**-2 for gi in supercell._g])*v






    
Out[14]:





0.055140235057336791

In [39]:

    

# Define the Gaussian broadening parameters
sigma1 = 0.75
sigma2 = 0.50

# Find the analytic solution.
U_analytic = ((1/sigma1 + 1/sigma2)/2 - np.sqrt(2)/
                np.sqrt(sigma1**2 + sigma2**2))/np.sqrt(np.pi)
Uerrors = []
ss = range(8,15)
for si in ss:
    # Create integer, sampling matrices.
    S = [si,si,si]

    # Define Coordinates.
    lat_const = 6
    r = [lat_const,0,0]
    R = makeR(r)
    origin = 1./2*np.sum(R,1)
    
    # Create the crystal
    crystal = Crystal(R, S)
    g(crystal)
    
    # Create vector of basis expansion coefficients for the electron density.
    n = -gauss_charge_dist(crystal._r, origin,sigma1) + (
        gauss_charge_dist(crystal._r,origin,sigma2))

    ϕ = -4*np.pi*I(crystal, Linv(crystal, O(crystal, J(crystal, n))))
    U_numeric = 1./2*np.real(np.dot(J(crystal, ϕ),O(crystal, J(crystal, n))))


    Uerrors.append(abs(U_numeric-U_analytic))
    del crystal
plt.semilogy(ss, Uerrors)
plt.xlabel("Sampling points per edge")
plt.ylabel("Error")
plt.show()
plt.close()


    

    




/Users/jeremy/Codes/dft/pydft/bases/planewave.py:74: RuntimeWarning: divide by zero encountered in double_scalars
  np.linalg.norm(gi)**-2 for gi in crystal._g])*v
/Users/jeremy/Codes/dft/pydft/bases/planewave.py:74: RuntimeWarning: invalid value encountered in multiply
  np.linalg.norm(gi)**-2 for gi in crystal._g])*v






    

In [6]:

    

from pydft.crystal import *
from pydft.energies.ewald import *

import numpy as np

# We should get something close to what Arias got with similar parameters.
LC = 6
r = [LC,0,0]
R = makeR(r)
S = [15,15,15]
crystal = Crystal(R, S, LC)
crystal._X = [[0., 0., 0.]]
energy = ewald_energy(crystal)
print(energy)


    

    




-0.32695557836

In [9]:

    

error = abs(energy - (-1./3))
print(error)


    

    




0.0063777549733

In [99]:

    

from pydft.energies.schrodinger import *
from pydft.crystal import *


    

In [104]:

    

LC = 1.
R = makeR([LC, 0., 0.])
S = [2,2,2]
crystal = Crystal(R, S, LC)
basis = "planewave"
potential = "sho"
w = 1.
E = getE(crystal, basis, potential, [w])


    

In [105]:

    

E


    

    
Out[105]:





(238.37050562614488-5.0115250491145069e-13j)