add dft-scripts and new images

.gitignore

@@ -24,3 +24,8 @@ molecules
 ltxpnggh-pages
 gh-pages
 ltxpng
+/auto/dft.el
+/dft.brf
+/dft.tex
+/mk4ht.cfg
+/topics-to-discuss.org

README.org

@@ -10,3 +10,8 @@ You are welcome to contribute examples of how to run calculations.
 Bug/typo/error reports should be sent to:
 
 John Kitchin (jkitchin@andrew.cmu.edu) or reported as an [[https://github.com/jkitchin/dft-book/issues][issue]].
+
+* TODO get this article for dft course
+
+
+http://pubs.rsc.org/en/Content/ArticleLanding/2014/GC/C3GC41368C#!divAbstract

dft-scripts/script-100.py

@@ -21,38 +21,24 @@
                Atom('O', (0.5 + u) * a1 + (0.5 - u) * a2 + 0.5 * a3),
                Atom('O', (0.5 - u) * a1 + (0.5 + u) * a2 + 0.5 * a3)],
               cell=[a1, a2, a3])
-v0 = atoms.get_volume()
-cell0 = atoms.get_cell()
-factors = [0.9, 0.95, 1.0, 1.05, 1.1] #to change volume by
-energies, volumes = [], []
+KPOINTS = [2, 3, 4, 5, 6, 7, 8]
+energies = []
 ready = True
-for f in factors:
-    v1 = f*v0
-    cell_factor = (v1 / v0)**(1. / 3.)
-    atoms.set_cell(cell0 * cell_factor, scale_atoms=True)
-    with jasp('bulk/tio2/step1-{0:1.2f}'.format(f),
+for k in KPOINTS:
+    with jasp('bulk/tio2/kpts-{0}'.format(k),
               encut=520,
-              kpts=(5,5,5),
-              isif=2, # relax internal degrees of freedom
-              ibrion=1,
-              nsw=50,
+              kpts=(k, k, k),
               xc='PBE',
               sigma=0.05,
               atoms=atoms) as calc:
         try:
             energies.append(atoms.get_potential_energy())
-            volumes.append(atoms.get_volume())
         except (VaspSubmitted, VaspQueued):
             ready = False
 if not ready:
     import sys; sys.exit()
-plt.plot(volumes, energies)
-plt.xlabel('Vol. ($\AA^3)$')
+plt.plot(KPOINTS, energies)
+plt.xlabel('number of k-points in each vector')
 plt.ylabel('Total energy (eV)')
-plt.savefig('images/tio2-step1.png')
-print '#+tblname: tio2-vol-ene'
-print '#+caption: Total energy of TiO_{2} vs. volume.'
-print '| Volume ($\AA^3$) | Energy (eV) |'
-print '|-'
-for v, e in zip(volumes, energies):
-    print '| {0} | {1} |'.format(v, e)
+plt.savefig('images/tio2-kpt-convergence.png')
+plt.show()

dft-scripts/script-101.py

@@ -1,29 +1,58 @@
+# step 1 frozen atoms and shape at different volumes
+from ase import Atom, Atoms
+import numpy as np
 from jasp import *
+import matplotlib.pyplot as plt
+'''
+create a TiO2 structure from the lattice vectors at
+http://cst-www.nrl.navy.mil/lattice/struk/c4.html
+'''
+a = 4.59 # experimental degrees of freedom.
+c = 2.96
+u = 0.3 #internal degree of freedom!
+#primitive vectors
+a1 = a*np.array([1.0, 0.0, 0.0])
+a2 = a*np.array([0.0, 1.0, 0.0])
+a3 = c*np.array([0.0, 0.0, 1.0])
+atoms = Atoms([Atom('Ti', [0., 0., 0.]),
+               Atom('Ti', 0.5 * a1 + 0.5 * a2 + 0.5 * a3),
+               Atom('O', u * a1 + u * a2),
+               Atom('O', -u * a1 - u * a2),
+               Atom('O', (0.5 + u) * a1 + (0.5 - u) * a2 + 0.5 * a3),
+               Atom('O', (0.5 - u) * a1 + (0.5 + u) * a2 + 0.5 * a3)],
+              cell=[a1, a2, a3])
+v0 = atoms.get_volume()
+cell0 = atoms.get_cell()
 factors = [0.9, 0.95, 1.0, 1.05, 1.1] #to change volume by
-energies1, volumes1 = [], [] # from step 1
-energies, volumes = [], [] # for step 2
+energies, volumes = [], []
 ready = True
 for f in factors:
-    with jasp('bulk/tio2/step1-{0:1.2f}'.format(f)) as calc:
-        atoms = calc.get_atoms()
-        energies1.append(atoms.get_potential_energy())
-        volumes1.append(atoms.get_volume())
-        calc.clone('bulk/tio2/step2-{0:1.2f}'.format(f))
-    # now set ISIF=4 and run
-    with jasp('bulk/tio2/step2-{0:1.2f}'.format(f),
-              isif=4) as calc:
-        atoms = calc.get_atoms()
+    v1 = f*v0
+    cell_factor = (v1 / v0)**(1. / 3.)
+    atoms.set_cell(cell0 * cell_factor, scale_atoms=True)
+    with jasp('bulk/tio2/step1-{0:1.2f}'.format(f),
+              encut=520,
+              kpts=(5,5,5),
+              isif=2, # relax internal degrees of freedom
+              ibrion=1,
+              nsw=50,
+              xc='PBE',
+              sigma=0.05,
+              atoms=atoms) as calc:
         try:
             energies.append(atoms.get_potential_energy())
             volumes.append(atoms.get_volume())
         except (VaspSubmitted, VaspQueued):
             ready = False
 if not ready:
     import sys; sys.exit()
-import matplotlib.pyplot as plt
-plt.plot(volumes1, energies1, volumes, energies)
+plt.plot(volumes, energies)
 plt.xlabel('Vol. ($\AA^3)$')
 plt.ylabel('Total energy (eV)')
-plt.legend(['step 1', 'step 2'], loc='best')
-plt.savefig('images/tio2-step2.png')
-plt.show()
+plt.savefig('images/tio2-step1.png')
+print '#+tblname: tio2-vol-ene'
+print '#+caption: Total energy of TiO_{2} vs. volume.'
+print '| Volume ($\AA^3$) | Energy (eV) |'
+print '|-'
+for v, e in zip(volumes, energies):
+    print '| {0} | {1} |'.format(v, e)

dft-scripts/script-102.py

@@ -1,10 +1,29 @@
 from jasp import *
-with jasp('bulk/tio2/step2-1.05') as calc:
-    calc.clone('bulk/tio2/step3')
-with jasp('bulk/tio2/step3',
-          isif=3) as calc:
-    calc.calculate()
-    atoms = calc.get_atoms()
-    print calc
-from pyspglib import spglib
-print '\nThe spacegroup is {0}'.format(spglib.get_spacegroup(atoms))
+factors = [0.9, 0.95, 1.0, 1.05, 1.1] #to change volume by
+energies1, volumes1 = [], [] # from step 1
+energies, volumes = [], [] # for step 2
+ready = True
+for f in factors:
+    with jasp('bulk/tio2/step1-{0:1.2f}'.format(f)) as calc:
+        atoms = calc.get_atoms()
+        energies1.append(atoms.get_potential_energy())
+        volumes1.append(atoms.get_volume())
+        calc.clone('bulk/tio2/step2-{0:1.2f}'.format(f))
+    # now set ISIF=4 and run
+    with jasp('bulk/tio2/step2-{0:1.2f}'.format(f),
+              isif=4) as calc:
+        atoms = calc.get_atoms()
+        try:
+            energies.append(atoms.get_potential_energy())
+            volumes.append(atoms.get_volume())
+        except (VaspSubmitted, VaspQueued):
+            ready = False
+if not ready:
+    import sys; sys.exit()
+import matplotlib.pyplot as plt
+plt.plot(volumes1, energies1, volumes, energies)
+plt.xlabel('Vol. ($\AA^3)$')
+plt.ylabel('Total energy (eV)')
+plt.legend(['step 1', 'step 2'], loc='best')
+plt.savefig('images/tio2-step2.png')
+plt.show()

dft-scripts/script-103.py

@@ -1,10 +1,10 @@
 from jasp import *
-with jasp('bulk/tio2/step3') as calc:
+with jasp('bulk/tio2/step2-1.05') as calc:
+    calc.clone('bulk/tio2/step3')
+with jasp('bulk/tio2/step3',
+          isif=3) as calc:
+    calc.calculate()
     atoms = calc.get_atoms()
-    print 'default ismear: ',atoms.get_potential_energy()
-    calc.clone('bulk/tio2/step4')
-with jasp('bulk/tio2/step4',
-          ismear=-5,
-          nsw=0) as calc:
-    atoms = calc.get_atoms()
-    print 'ismear=-5:      ',atoms.get_potential_energy()
+    print calc
+from pyspglib import spglib
+print '\nThe spacegroup is {0}'.format(spglib.get_spacegroup(atoms))

dft-scripts/script-104.py

@@ -1,25 +1,10 @@
 from jasp import *
-from ase import Atom, Atoms
-from ase.utils.eos import EquationOfState
-LC = [3.75, 3.80, 3.85, 3.90, 3.95, 4.0, 4.05, 4.1]
-volumes, energies = [],[]
-for a in LC:
-    atoms = Atoms([Atom('Pd', (0, 0, 0))],
-                  cell=0.5 * a*np.array([[1.0, 1.0, 0.0],
-                                         [0.0, 1.0, 1.0],
-                                         [1.0, 0.0, 1.0]]))
-    with jasp('bulk/Pd-LDA-{0}'.format(a),
-              encut=350,
-              kpts=(12,12,12),
-              xc='LDA',
-              atoms=atoms):
-        try:
-            e = atoms.get_potential_energy()
-            energies.append(e)
-            volumes.append(atoms.get_volume())
-        except (VaspSubmitted, VaspQueued):
-            pass
-if len(energies) == len(LC):
-    eos = EquationOfState(volumes, energies)
-    v0, e0, B = eos.fit()
-    print 'LDA lattice constant is {0:1.3f} Ang^3'.format((4*v0)**(1./3.))
+with jasp('bulk/tio2/step3') as calc:
+    atoms = calc.get_atoms()
+    print 'default ismear: ',atoms.get_potential_energy()
+    calc.clone('bulk/tio2/step4')
+with jasp('bulk/tio2/step4',
+          ismear=-5,
+          nsw=0) as calc:
+    atoms = calc.get_atoms()
+    print 'ismear=-5:      ',atoms.get_potential_energy()

dft-scripts/script-105.py

@@ -2,34 +2,24 @@
 from ase import Atom, Atoms
 from ase.utils.eos import EquationOfState
 LC = [3.75, 3.80, 3.85, 3.90, 3.95, 4.0, 4.05, 4.1]
-GGA = {'AM':'AM05',
-       'PE':'PBE',
-       'PS':'PBEsol',
-       'RP':'RPBE'}
-for key in GGA:
-    volumes, energies = [],[]
-    for a in LC:
-        atoms = Atoms([Atom('Pd', (0, 0, 0))],
-                      cell=0.5 * a*np.array([[1.0, 1.0, 0.0],
-                                             [0.0, 1.0, 1.0],
-                                             [1.0, 0.0, 1.0]]))
-        with jasp('bulk/Pd-GGA-{1}-{0}'.format(a,key),
-                  encut=350,
-                  kpts=(12,12,12),
-                  xc='LDA',
-                  gga=key,
-                  atoms=atoms):
-            try:
-                e = atoms.get_potential_energy()
-                energies.append(e)
-                volumes.append(atoms.get_volume())
-            except (VaspSubmitted, VaspQueued):
-                pass
-    if len(energies) == len(LC):
-        eos = EquationOfState(volumes, energies)
-        v0, e0, B = eos.fit()
-        print '{1:6s} lattice constant is {0:1.3f} Ang^3'.format((4*v0)**(1./3.),
-                                                             GGA[key])
-    else:
-        print energies, LC
-        print '{0} is not ready'.format(GGA[key])
+volumes, energies = [],[]
+for a in LC:
+    atoms = Atoms([Atom('Pd', (0, 0, 0))],
+                  cell=0.5 * a*np.array([[1.0, 1.0, 0.0],
+                                         [0.0, 1.0, 1.0],
+                                         [1.0, 0.0, 1.0]]))
+    with jasp('bulk/Pd-LDA-{0}'.format(a),
+              encut=350,
+              kpts=(12,12,12),
+              xc='LDA',
+              atoms=atoms):
+        try:
+            e = atoms.get_potential_energy()
+            energies.append(e)
+            volumes.append(atoms.get_volume())
+        except (VaspSubmitted, VaspQueued):
+            pass
+if len(energies) == len(LC):
+    eos = EquationOfState(volumes, energies)
+    v0, e0, B = eos.fit()
+    print 'LDA lattice constant is {0:1.3f} Ang^3'.format((4*v0)**(1./3.))

dft-scripts/script-106.py

@@ -1,29 +1,35 @@
 from jasp import *
-from ase.lattice.cubic import FaceCenteredCubic
-from ase import Atoms, Atom
-# bulk system
-atoms = FaceCenteredCubic(directions=[[0,1,1],
-                                      [1,0,1],
-                                      [1,1,0]],
-                                      size=(1,1,1),
-                                      symbol='Rh')
-with jasp('bulk/bulk-rh',
-          xc='PBE',
-          encut=350,
-          kpts=(4,4,4),
-          isif=3,
-          ibrion=2,
-          nsw=10,
-          atoms=atoms) as calc:
-    bulk_energy = atoms.get_potential_energy()
-# atomic system
-atoms = Atoms([Atom('Rh',[5, 5, 5])],
-              cell=(7, 8, 9))
-with jasp('bulk/atomic-rh',
-          xc='PBE',
-          encut=350,
-          kpts=(1, 1, 1),
-          atoms=atoms) as calc:
-    atomic_energy = atoms.get_potential_energy()
-cohesive_energy = atomic_energy - bulk_energy
-print 'The cohesive energy is {0:1.3f} eV'.format(cohesive_energy)
+from ase import Atom, Atoms
+from ase.utils.eos import EquationOfState
+LC = [3.75, 3.80, 3.85, 3.90, 3.95, 4.0, 4.05, 4.1]
+GGA = {'AM':'AM05',
+       'PE':'PBE',
+       'PS':'PBEsol',
+       'RP':'RPBE'}
+for key in GGA:
+    volumes, energies = [],[]
+    for a in LC:
+        atoms = Atoms([Atom('Pd', (0, 0, 0))],
+                      cell=0.5 * a*np.array([[1.0, 1.0, 0.0],
+                                             [0.0, 1.0, 1.0],
+                                             [1.0, 0.0, 1.0]]))
+        with jasp('bulk/Pd-GGA-{1}-{0}'.format(a,key),
+                  encut=350,
+                  kpts=(12,12,12),
+                  xc='LDA',
+                  gga=key,
+                  atoms=atoms):
+            try:
+                e = atoms.get_potential_energy()
+                energies.append(e)
+                volumes.append(atoms.get_volume())
+            except (VaspSubmitted, VaspQueued):
+                pass
+    if len(energies) == len(LC):
+        eos = EquationOfState(volumes, energies)
+        v0, e0, B = eos.fit()
+        print '{1:6s} lattice constant is {0:1.3f} Ang^3'.format((4*v0)**(1./3.),
+                                                             GGA[key])
+    else:
+        print energies, LC
+        print '{0} is not ready'.format(GGA[key])