• R/O
  • HTTP
  • SSH
  • HTTPS

pwdft: Commit

Source code and sample for Educational-PWDFT


Commit MetaInfo

Revision222a7b36ea78813caff0622ef392bcd20069a5a4 (tree)
Zeit2018-12-27 18:31:48
Autormitsuaki1987 <kawamitsuaki@gmai...>
Commitermitsuaki1987

Log Message

Add script to generate QE input from cif file

Ändern Zusammenfassung

Diff

--- a/.gitignore
+++ b/.gitignore
@@ -1,7 +1,8 @@
1-*.o
2-*.mod
3-*.x
4-make.inc
5-*~
6-doc/html/
7-*.a
1+*.o
2+*.mod
3+*.x
4+make.inc
5+*~
6+doc/html/
7+*.a
8+tool/__pycache__
--- /dev/null
+++ b/doc/index.html
@@ -0,0 +1,110 @@
1+<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
2+<html lang="ja">
3+ <head>
4+ <meta name="Content-Language" content="ja">
5+ <meta http-equiv="Content-Style-Type" content="text/css">
6+ <meta http-equiv="Content-Script-Type" content="text/javascript">
7+ <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
8+ <meta http-equiv="Pragma" content="no-cache">
9+ <meta http-equiv="Cache-Control" content="no-cache">
10+ <title>Web page for Numerical Analysis for Material Science 2 (DFT part)</title>
11+ </head>
12+
13+ <body bgcolor="CCFFCC">
14+ <h1>
15+ <p>Numerical Analysis for Material Science 2 (DFT part)</p>
16+ <p>物質科学のための計算数理 2 (DFT part)</p>
17+ </h1>
18+ <hr>
19+
20+ <h2><a href="./html/index.html">Source browser</a></h2>
21+ <h2><a href="https://osdn.net/projects/educational-pwdft/scm/git/pwdft/">Git browser</a></h2>
22+ <h2>
23+ <code>
24+ $ git clone git://git.osdn.net/gitroot/educational-pwdft/pwdft.git
25+ </code>
26+ </h2>
27+
28+ <ul>
29+ <li>
30+ <h2>8th : Basics of DFT calculation
31+ (<a href="2018W-8.pptx">PPTX</a> / <a href="2018W-8.pdf">PDF</a>)</h2>
32+ <h3>
33+ <ul>
34+ <li>Hohenberg-Kohn theorem</li>
35+ <li>Kohn-Sham method</li>
36+ <li>Local density approximation</li>
37+ <li>Bloch theorem</li>
38+ <li>Hands-on of QuantumESPRESSO : MgB<sub>2</sub></li>
39+ </ul>
40+ </h3>
41+ </li>
42+ <li>
43+ <h2>
44+ 9th : Solution of Kohn-Sham eq.
45+ (<a href="2018W-9.pptx">PPTX</a> / <a href="2018W-9.pdf">PDF</a>)
46+ </h2>
47+ <h3>
48+ <ul>
49+ <li>Plane wave representation of Kohn-Sham eq.</li>
50+ <li>Direct diagonalization</li>
51+ <li>Iterative method : LOBPCG</li>
52+ </ul>
53+ </h3>
54+ </li>
55+ <li>
56+ <h2>
57+ 10th : DFT self-consistent loop
58+ (<a href="2018W-10.pptx">PPTX</a> / <a href="2018W-10.pdf">PDF</a>)
59+ </h2>
60+ <h3>
61+ <ul>
62+ <li>SCF loop</li>
63+ <li>Kohn-Sham potential</li>
64+ <li>Broyden method</li>
65+ <li>Plot volumetric (grid) data</li>
66+ </ul>
67+ </h3>
68+ </li>
69+ <li>
70+ <h2>
71+ 11th Total energy
72+ (<a href="2018W-11.pptx">PPTX</a> / <a href="2018W-11.pdf">PDF</a>)
73+ </h2>
74+ <h3>
75+ <ul>
76+ <li>Total energy</li>
77+ <li>Ewald sum</li>
78+ <li>Tetrahedron method</li>
79+ <li>Fermi surface</li>
80+ </ul>
81+ </h3>
82+ </li>
83+ <li>
84+ <h2>
85+ 12th : Advanced topics
86+ (<a href="2018W-12.pptx">PPTX</a> / <a href="2018W-12.pdf">PDF</a>)
87+ </h2>
88+ <h3>
89+ <ul>
90+ <li>Non-local pseudopotentials</li>
91+ <li>Generalized gradient approximation</li>
92+ <li>Spin density functional theory</li>
93+ </ul>
94+ </h3>
95+ </li>
96+ <li>
97+ <h2>
98+ 13th : Practical calculation
99+ (<a href="2018W-13.pptx">PPTX</a> / <a href="2018W-13.pdf">PDF</a>)
100+ </h2>
101+ <h3>
102+ <ul>
103+ <li>CIF file and structure database</li>
104+ </ul>
105+ </h3>
106+ </li>
107+ </ul>
108+
109+ </body>
110+</html>
--- a/src/energy.F90
+++ b/src/energy.F90
@@ -155,7 +155,7 @@ contains
155155 norm = sqrt(dot_product(bvec(1:3,ii), bvec(1:3,ii)))
156156 eta = eta + norm
157157 end do
158- eta = eta / 3.0d0
158+ eta = eta / 3.0d0 / sqrt(2.0d0) / sqrt(2.0d0*pi)
159159 write(*,*) " Ewald eta [Bohr^-1] : ", eta
160160 !
161161 ! Reciprocal space
@@ -232,7 +232,7 @@ contains
232232 !
233233 end do
234234 !
235- write(*,*) " Eward energy : ", Eew*htr2ev
235+ write(*,*) " Ewald energy : ", Eew*htr2ev
236236 Etot = Etot + Eew
237237 !
238238 end subroutine ewald
--- /dev/null
+++ b/tool/cif2qe.py
@@ -0,0 +1,280 @@
1+#!/usr/bin/python3
2+import sys
3+import pymatgen
4+import seekpath
5+import numpy
6+from pymatgen.core.periodic_table import get_el_sp
7+import json
8+
9+
10+def write_atom(f, avec, typ, nat, pos, atom):
11+ print("CELL_PARAMETERS angstrom", file=f)
12+ for ii in range(3):
13+ print(" %f %f %f" % (avec[ii, 0], avec[ii, 1], avec[ii, 2]), file=f)
14+ print("ATOMIC_SPECIES", file=f)
15+ for ityp in typ:
16+ print(" %s %f %s" % (ityp, pymatgen.Element(ityp).atomic_mass, sssp[str(ityp)]["filename"]), file=f)
17+ print("ATOMIC_POSITIONS crystal", file=f)
18+ for iat in range(nat):
19+ print(" %s %f %f %f" % (
20+ atom[iat], pos[iat][0], pos[iat][1], pos[iat][2]), file=f)
21+
22+
23+def write_middle(f, nat, ntyp, ecutwfc, ecutrho, psdir):
24+ print(" pseudo_dir = \'%s/\'" % psdir, file=f)
25+ print("/", file=f)
26+ print("&SYSTEM", file=f)
27+ print(" ibrav = 0", file=f)
28+ print(" nat = %d" % nat, file=f)
29+ print(" ntyp = %d" % ntyp, file=f)
30+ print(" ecutwfc = %f" % ecutwfc, file=f)
31+ print(" ecutrho = %f" % ecutrho, file=f)
32+
33+
34+if __name__ == '__main__':
35+
36+ args = sys.argv
37+ if len(args) < 2:
38+ print("Usage:")
39+ print("$ cif2qe.py cif-file pseudo-directory [dk_path] [dk_grid]")
40+ print("Default:")
41+ print("$ cif2qe.py cif-file pseudo-directory 0.1 0.2")
42+ exit(0)
43+ #
44+ # CIF parser
45+ #
46+ structure = pymatgen.Structure.from_file(args[1])
47+ #
48+ # Default value
49+ #
50+ dk_path = 0.1
51+ dk_grid = 0.2
52+ #
53+ psdir = args[2]
54+ f = open(psdir + '/sssp_efficiency.json', 'r')
55+ sssp = json.load(f)
56+ #
57+ if len(args) > 3:
58+ dk_path = float(args[3])
59+ if len(args) > 4:
60+ dk_grid = float(args[4])
61+ #
62+ print(" dk for band : {0}".format(dk_path))
63+ print(" dk for grid : {0}".format(dk_grid))
64+
65+ structure.remove_oxidation_states()
66+
67+ #
68+ # Band path and primitive lattice
69+ #
70+ skp = seekpath.get_explicit_k_path((structure.lattice.matrix, structure.frac_coords,
71+ [pymatgen.Element(str(spc)).number for spc in structure.species]),
72+ reference_distance=dk_path)
73+ #
74+ # Lattice information
75+ #
76+ avec = skp["primitive_lattice"]
77+ bvec = skp["reciprocal_primitive_lattice"]
78+ pos = skp["primitive_positions"]
79+ nat = len(skp["primitive_types"])
80+ atom = [str(get_el_sp(iat)) for iat in skp["primitive_types"]]
81+ typ = set(atom)
82+ ntyp = len(typ)
83+ typ = set(atom)
84+ #
85+ # WFC and Rho cutoff
86+ #
87+ ecutwfc = 0.0
88+ ecutrho = 0.0
89+ for ityp in typ:
90+ if ecutwfc < sssp[str(ityp)]["cutoff"]:
91+ ecutwfc = sssp[str(ityp)]["cutoff"]
92+ if ecutrho < sssp[str(ityp)]["cutoff"] * sssp[str(ityp)]['dual']:
93+ ecutrho = sssp[str(ityp)]["cutoff"] * sssp[str(ityp)]['dual']
94+ #
95+ # k and q grid
96+ #
97+ nk = numpy.zeros(3, numpy.int_)
98+ for ii in range(3):
99+ norm = numpy.sqrt(numpy.dot(bvec[ii][:], bvec[ii][:]))
100+ nk[ii] = round(norm / dk_grid)
101+ print(norm)
102+ print("SCF k-grid : ", nk[0], nk[1], nk[2])
103+ #
104+ # Band path
105+ #
106+ print("Band path")
107+ for ipath in range(len(skp["path"])):
108+ start = skp["explicit_segments"][ipath][0]
109+ final = skp["explicit_segments"][ipath][1] - 1
110+ print("%5d %8s %10.5f %10.5f %10.5f %8s %10.5f %10.5f %10.5f" % (
111+ final - start + 1,
112+ skp["explicit_kpoints_labels"][start],
113+ skp["explicit_kpoints_rel"][start][0],
114+ skp["explicit_kpoints_rel"][start][1],
115+ skp["explicit_kpoints_rel"][start][2],
116+ skp["explicit_kpoints_labels"][final],
117+ skp["explicit_kpoints_rel"][final][0],
118+ skp["explicit_kpoints_rel"][final][1],
119+ skp["explicit_kpoints_rel"][final][2]))
120+ #
121+ # Number of electrons
122+ #
123+ nbnd = len(atom)*10
124+ #
125+ # scf.in : Charge density
126+ #
127+ with open("scf.in", 'w') as f:
128+ print("&CONTROL", file=f)
129+ print(" calculation = \'scf\'", file=f)
130+ write_middle(f, nat, ntyp, ecutwfc, ecutrho, psdir)
131+ print(" occupations = \'tetrahedra_opt\'", file=f)
132+ print("/", file=f)
133+ #
134+ print("&ELECTRONS", file=f)
135+ print(" mixing_beta = 0.3", file=f)
136+ print(" conv_thr = %e" % (float(nat)*1.0e-10), file=f)
137+ print("/", file=f)
138+ write_atom(f, avec, typ, nat, pos, atom)
139+ print("K_POINTS automatic", file=f)
140+ print(" %d %d %d 0 0 0" % (nk[0], nk[1], nk[2]), file=f)
141+ #
142+ # nscf.in : Dense k grid
143+ #
144+ with open("nscf.in", 'w') as f:
145+ print("&CONTROL", file=f)
146+ print(" calculation = \'nscf\'", file=f)
147+ write_middle(f, nat, ntyp, ecutwfc, ecutrho, psdir)
148+ print(" occupations = \'tetrahedra_opt\'", file=f)
149+ print(" nbnd = %d" % nbnd, file=f)
150+ print(" la2f = .true.", file=f)
151+ print("/", file=f)
152+ print("&ELECTRONS", file=f)
153+ print("/", file=f)
154+ write_atom(f, avec, typ, nat, pos, atom)
155+ print("K_POINTS automatic", file=f)
156+ print(" %d %d %d 0 0 0" % (nk[0]*2, nk[1]*2, nk[2]*2), file=f)
157+ #
158+ # band.in : Plot band
159+ #
160+ with open("band.in", 'w') as f:
161+ print("&CONTROL", file=f)
162+ print(" calculation = \'bands\'", file=f)
163+ write_middle(f, nat, ntyp, ecutwfc, ecutrho, psdir)
164+ print(" nbnd = %d" % nbnd, file=f)
165+ print("/", file=f)
166+ print("&ELECTRONS", file=f)
167+ print("/", file=f)
168+ write_atom(f, avec, typ, nat, pos, atom)
169+ print("K_POINTS crystal", file=f)
170+ print(len(skp["explicit_kpoints_rel"]), file=f)
171+ for ik in range(len(skp["explicit_kpoints_rel"])):
172+ print(" %f %f %f 1.0" % (
173+ skp["explicit_kpoints_rel"][ik][0],
174+ skp["explicit_kpoints_rel"][ik][1],
175+ skp["explicit_kpoints_rel"][ik][2]),
176+ file=f)
177+ #
178+ # bands.in : Read by bands.x
179+ #
180+ with open("bands.in", 'w') as f:
181+ print("&BANDS", file=f)
182+ print(" lsym = .false.", file=f)
183+ print("/", file=f)
184+ #
185+ # proj.in : Read by projwfc.x
186+ #
187+ with open("proj.in", 'w') as f:
188+ print("&PROJWFC", file=f)
189+ print(" emin = ", file=f)
190+ print(" emax = ", file=f)
191+ print(" deltae = 0.1", file=f)
192+ print("/", file=f)
193+ #
194+ # pp.in : Plot Kohn-Sham orbitals
195+ #
196+ with open("pp.in", 'w') as f:
197+ print("&INPUTPP ", file=f)
198+ print(" plot_num = 7", file=f)
199+ print(" kpoint = 1", file=f)
200+ print(" kband(1) = %d" % 1, file=f)
201+ print(" kband(2) = %d" % nbnd, file=f)
202+ print(" lsign = .true.", file=f)
203+ print("/", file=f)
204+ print("&PLOT ", file=f)
205+ print(" iflag = 3", file=f)
206+ print(" output_format = 5", file=f)
207+ print(" fileout = \".xsf\"", file=f)
208+ print("/", file=f)
209+ #
210+ # band.gp : Gnuplot script
211+ #
212+ with open("band.gp", 'w') as f:
213+ print("set terminal pdf color enhanced \\", file=f)
214+ print("dashed dl 0.5 size 8.0cm, 6.0cm", file=f)
215+ print("set output \"band.pdf\"", file=f)
216+ print("#", file=f)
217+ n_sym_points = 1
218+ final = 0
219+ x0 = numpy.linalg.norm(avec[0, :]) * 0.5 / numpy.pi
220+ print("x%d = %f" % (n_sym_points, x0*skp["explicit_kpoints_linearcoord"][final]), file=f)
221+ for ipath in range(len(skp["path"])):
222+ start = skp["explicit_segments"][ipath][0]
223+ if start != final:
224+ n_sym_points += 1
225+ print("x%d = %f" % (n_sym_points, x0*skp["explicit_kpoints_linearcoord"][start]), file=f)
226+ n_sym_points += 1
227+ final = skp["explicit_segments"][ipath][1] - 1
228+ print("x%d = %f" % (n_sym_points, x0*skp["explicit_kpoints_linearcoord"][final]), file=f)
229+ print("#", file=f)
230+ print("set border lw 2", file=f)
231+ print("#", file=f)
232+ print("set style line 1 lt 1 lw 2 lc 0 dashtype 2", file=f)
233+ print("set style line 2 lt 1 lw 2 lc 0", file=f)
234+ print("set style line 3 lt 1 lw 1 lc 1", file=f)
235+ print("#", file=f)
236+ print("set ytics scale 3.0, -0.5 font \'Cmr10,18\'", file=f)
237+ print("set xtics( \\", file=f)
238+ n_sym_points = 1
239+ final = 0
240+ label_f = skp["explicit_kpoints_labels"][final]
241+ if label_f == "GAMMA":
242+ label_f = "\\241"
243+ print("\"%s\" x%d" % (label_f, n_sym_points), end="", file=f)
244+ for ipath in range(len(skp["path"])):
245+ start = skp["explicit_segments"][ipath][0]
246+ label_s = skp["explicit_kpoints_labels"][start]
247+ if label_s == "GAMMA":
248+ label_s = "\\241"
249+ label_f = skp["explicit_kpoints_labels"][final]
250+ if label_f == "GAMMA":
251+ label_f = "\\241"
252+ if start != final:
253+ n_sym_points += 1
254+ print(", \\\n\"%s%s\" x%d" % (label_f, label_s, n_sym_points), end="", file=f)
255+ n_sym_points += 1
256+ final = skp["explicit_segments"][ipath][1] - 1
257+ label_f = skp["explicit_kpoints_labels"][final]
258+ if label_f == "GAMMA":
259+ label_f = "\\241"
260+ print(", \\\n\"%s\" x%d" % (label_f, n_sym_points), end="", file=f)
261+ print(") \\\noffset 0.0, 0.0 font \'Cmr10,18\'", file=f)
262+ print("#", file=f)
263+ print("set grid xtics ls 2 front", file=f)
264+ print("#", file=f)
265+ print("unset key", file=f)
266+ print("#", file=f)
267+ print("set xzeroaxis ls 1", file=f)
268+ print("#", file=f)
269+ print("set ylabel \"Energy from {/Cmmi10 \\042}_F [eV]\" offset - 0.5, 0.0 font \'Cmr10,18\'", file=f)
270+ print("#", file=f)
271+ n_sym_points = 1
272+ final = 0
273+ for ipath in range(len(skp["path"])):
274+ start = skp["explicit_segments"][ipath][0]
275+ if start == final:
276+ n_sym_points += 1
277+ final = skp["explicit_segments"][ipath][1] - 1
278+ else:
279+ break
280+ print("plot[:][emin:emax] \"bands.out.gnu\" u 1:($2-ef) w l ls 3", file=f)
Show on old repository browser