Easy Ab initio calculation with ASE-Siesta-Pyscf

No installation necessary, just download a ready to go container for any system, or run it into the cloud

We first import the necessary libraries and define the system using ASE



In [ ]:

    
# import libraries and set up the molecule geometry

from ase.units import Ry, eV, Ha
from ase.calculators.siesta import Siesta
from ase import Atoms
import numpy as np
import matplotlib.pyplot as plt

from ase.build import molecule

CH4 = molecule("CH3COOH")

# visualization of the particle
from ase.visualize import view
view(CH4, viewer='x3d')

We can then run the DFT calculation using Siesta



In [ ]:

    
# enter siesta input and run siesta
siesta = Siesta(
    mesh_cutoff=150 * Ry,
    basis_set='DZP',
    pseudo_qualifier='lda',
    energy_shift=(10 * 10**-3) * eV,
    fdf_arguments={
        'SCFMustConverge': False,
        'COOP.Write': True,
        'WriteDenchar': True,
        'PAO.BasisType': 'split',
        'DM.Tolerance': 1e-4,
        'DM.MixingWeight': 0.1,
        'MaxSCFIterations': 300,
        'DM.NumberPulay': 4,
        'XML.Write': True})

CH4.set_calculator(siesta)
e = CH4.get_potential_energy()

The TDDFT calculations with PySCF-NAO



In [ ]:

    
# compute polarizability using pyscf-nao

freq = np.arange(0.0, 15.0, 0.05)
siesta.pyscf_tddft(label="siesta", jcutoff=7, iter_broadening=0.15/Ha,
        xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7, freq = freq)



In [ ]:

    
# plot polarizability with matplotlib
%matplotlib inline

fig = plt.figure(1)
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.plot(siesta.results["freq range"], siesta.results["polarizability nonin"][:, 0, 0].imag)
ax2.plot(siesta.results["freq range"], siesta.results["polarizability inter"][:, 0, 0].imag)

ax1.set_xlabel(r"$\omega$ (eV)")
ax2.set_xlabel(r"$\omega$ (eV)")

ax1.set_ylabel(r"Im($P_{xx}$) (au)")
ax2.set_ylabel(r"Im($P_{xx}$) (au)")

ax1.set_title(r"Non interacting")
ax2.set_title(r"Interacting")

fig.tight_layout()



In [ ]:

    
import plotly.offline as py
import plotly.graph_objs as go

py.init_notebook_mode(connected=True)

trace = go.Scatter(x=siesta.results["freq range"], y=siesta.results["polarizability inter"][:, 0, 0].imag)

print(CH4.get_positions())
py.iplot([trace])

Compute the spatial distributoin of the density change at resonance frequency



In [ ]:

    
res = 10.5/Ha 
lim = 20.0 # Bohr
box = np.array([[-lim, lim],
                [-lim, lim],
                [-lim, lim]])
from pyscf.nao.m_comp_spatial_distributions import spatial_distribution

spd = spatial_distribution(siesta.results["density change inter"], freq/Ha, box, label="siesta")
spd.get_spatial_density(8.25/Ha)



In [ ]:

    
center = np.array([spd.dn_spatial.shape[0]/2, spd.dn_spatial.shape[1]/2, spd.dn_spatial.shape[2]/2], dtype=int)



In [ ]:

    
fig2 = plt.figure(2, figsize=(15, 12))

cmap="seismic"
ax1 = fig2.add_subplot(1, 3, 1)
vmax = np.max(abs(spd.dn_spatial[center[0], :, :].imag))
vmin = -vmax
ax1.imshow(spd.dn_spatial[center[0], :, :].imag, interpolation="bicubic", vmin=vmin, vmax=vmax, cmap=cmap, extent=[spd.mesh[1][0], spd.mesh[1][spd.mesh[1].shape[0]-1], spd.mesh[2][0], spd.mesh[2][spd.mesh[2].shape[0]-1]])

ax2 = fig2.add_subplot(1, 3, 2)
vmax = np.max(abs(spd.dn_spatial[:, center[1], :].imag))
vmin = -vmax
ax2.imshow(spd.dn_spatial[:, center[1], :].imag, interpolation="bicubic", vmin=vmin, vmax=vmax, cmap=cmap, extent=[spd.mesh[0][0], spd.mesh[0][spd.mesh[0].shape[0]-1], spd.mesh[2][0], spd.mesh[2][spd.mesh[2].shape[0]-1]])

ax3 = fig2.add_subplot(1, 3, 3)
vmax = np.max(abs(spd.dn_spatial[:, :, center[2]].imag))
vmin = -vmax
ax3.imshow(spd.dn_spatial[:, :, center[2]].imag, interpolation="bicubic", vmin=vmin, vmax=vmax, cmap=cmap, extent=[spd.mesh[0][0], spd.mesh[0][spd.mesh[0].shape[0]-1], spd.mesh[1][0], spd.mesh[1][spd.mesh[1].shape[0]-1]])

ax1.set_xlabel(r"y (Bohr)")
ax2.set_xlabel(r"x (Bohr)")
ax3.set_xlabel(r"x (Bohr)")

ax1.set_ylabel(r"z (Bohr)")
ax2.set_ylabel(r"z (Bohr)")
ax3.set_ylabel(r"y (Bohr)")

ax1.set_title(r"Im($\delta n$) in the $x$ plane")
ax2.set_title(r"Im($\delta n$) in the $y$ plane")
ax3.set_title(r"Im($\delta n$) in the $z$ plane")

fig2.savefig("density_chang_im_CH3COOH.pdf", format="pdf")



In [ ]: