bugfix

tools/vasp/wisdom-sampling-2.py

@@ -0,0 +1,83 @@
+import random, math, os, copy, time
+random.seed(time.time())
+
+f = open("POSCAR", "r")
+lines = f.readlines()
+f.close()
+#line[0] comments
+#line[1] scale factor
+#line[2-4] pbc
+#line[5] elements
+#line[6] natoms
+#line[7] Selective dynamics
+#line[8] label for coordination type
+#line[9] coordination
+
+head = []
+coords = []
+if lines[7].startswith("Selective"): 
+    for i in range(9):
+        head.append(lines[i])
+natoms = 0
+for i in lines[6].strip().split():
+    natoms += int(i)
+for i in range(9,9+natoms):
+    coords.append(lines[i])
+
+x = [0.0, 1.0, 0.0, 1.0, 0.21, 0.23]
+#nsamples = 200
+#nx =  int(math.sqrt(nsamples)) + 1
+nx = 4
+ny = 4
+dx = x[1] - x[0]
+dy = x[3] - x[2]
+dz = x[5] - x[4]
+dx = dx/nx
+dy = dy/nx
+
+shift = 0.125
+
+coords_new = []
+for i in range(4):
+    for j in range(4):
+        xnew = 0.25*i + shift + (0.04 - random.random()*0.4)
+        ynew = 0.25*j + (0.04 - random.random()*0.4)
+        znew = 0.22 + (0.02 - random.random()*0.02)
+        coords_new.append([xnew, ynew, znew])
+
+for i in range(4):
+    for j in range(4):
+        xnew = 0.25*i + (0.04 - random.random()*0.4)
+        ynew = 0.25*j + shift + (0.04 - random.random()*0.4)
+        znew = 0.22 + (0.02 - random.random()*0.02)
+        coords_new.append([xnew, ynew, znew])
+
+for i in range(len(coords_new)):
+    folder = "r_%03d"%i
+    if not os.path.exists(folder):
+        os.mkdir(folder)
+    os.chdir(folder)
+    o = open("POSCAR", "w")
+    new_head = copy.copy(head)
+    new_head[5] = new_head[5].strip() + "    H\n"
+    new_head[6] = new_head[6].strip() + "    1\n"
+    newline = ""
+    for j in coords_new[i]:
+        newline += "%18.8f"%j
+    newline += " T T T\n"
+    for i in new_head:
+        o.write(i)
+    for i in coords:
+        line = i.replace("T", "F")
+        o.write(line)
+    o.write(newline)
+    o.close()
+    os.chdir("..")
+
+if not os.path.exists("ref"):
+    os.mkdir("ref")
+os.chdir("ref")
+o = open("POSCAR", "w")
+for i in lines:
+    o.write(i)
+o.close()

tools/vasp/z-dist/cal-hist.py

@@ -0,0 +1,27 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+data = np.loadtxt("all.dat")
+data = data.transpose()
+
+hist_data = []
+n_total = 0
+for i in range(len(data[0])):
+    ndata = int(data[1][i])
+    n_total += ndata
+    for j in range(ndata):
+        hist_data.append(data[0][i])
+
+hist, bin_edges = np.histogram(hist_data, 12)
+dx = bin_edges[1] - bin_edges[0]
+o = open("all-hist.dat", "w")
+for i in range(len(hist)):
+    o.write("%12.8f%12.4f\n"%(bin_edges[i]+0.5*dx, 
+                              hist[i]*36.0/n_total))
+o.close()
+
+plt.plot(bin_edges[:-1], hist*36.0/n_total)
+plt.show()
+
+
+

tools/vasp/z-dist/gaussian-fit.py

@@ -0,0 +1,35 @@
+import sys
+import numpy as np
+from numpy import exp, loadtxt, pi, sqrt
+import matplotlib.pyplot as plt
+from lmfit import Model
+
+def gaussian(x, amp, cen, wid):
+    """1-d gaussian: gaussian(x, amp, cen, wid)"""
+    return (amp / (sqrt(2*pi) * wid)) * exp(-(x-cen)**2 / (2*wid**2))
+
+
+def gaussian_fit(fname):
+    fname = sys.argv[1]
+
+    data = np.loadtxt(fname)
+    data = data.transpose()
+    x = data[0]
+    y = data[1]
+
+    gmodel = Model(gaussian)
+    result = gmodel.fit(y, x=x, amp=120, cen=9.1, wid=0.2)
+    print(result.fit_report())
+    #plt.plot(data[0], data[1])
+    #plt.plot(x, result.init_fit, 'k--')
+    #plt.plot(x, result.best_fit, 'r-')
+    #plt.show()
+
+def main():
+    if len(sys.argv) > 1:
+        fname = sys.argv[1]
+        gaussian_fit(fname)
+
+if __name__ == "__main__":
+    main()
+

tools/vasp/z-dist/peak-fit.py

@@ -0,0 +1,36 @@
+"""
+@ref: https://plot.ly/python/peak-fitting/
+"""
+import sys
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy
+import peakutils
+
+from scipy import signal
+
+def gaussian(x, mu, sig):
+    return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
+
+def peak_fit(fname):
+    data = np.loadtxt(fname)
+    data = data.transpose()
+
+    indices = peakutils.indexes(data[1], thres=0.0, min_dist=0.05)
+
+    peaks = [[],[]]
+    for i in range(len(indices)):
+        peaks[0].append(data[0][indices[i]])
+        peaks[1].append(data[1][indices[i]])
+
+    if 0:
+        plt.plot(data[0], data[1])
+        plt.scatter(peaks[0], peaks[1]) 
+        plt.show()
+
+def main():
+    if len(sys.argv) > 1:
+        fname = peak_fit(sys.argv[1])
+
+if __name__ == "__main__":
+    main()

tools/vasp/z-dist/xyz2z-profiles.py

@@ -0,0 +1,35 @@
+import numpy as np
+
+f = open("movie.xyz", "r")
+lines = f.readlines()
+f.close()
+
+n_atom = int(lines[0])
+n_frame = len(lines)/(n_atom + 2)
+
+o_list = []
+for i in range(n_atom+2):
+    tokens = lines[i].strip().split()
+    if len(tokens) == 4:
+        if tokens[0] == "O":
+            o_list.append(i)
+
+o_atoms = []
+for i in range(len(o_list)):
+    o_atoms.append([])
+
+for i in range(n_frame):
+    for j in range(len(o_list)):
+        tokens = lines[i*(n_atom + 2)+o_list[j]].strip().split()
+        z = float(tokens[-1])
+        o_atoms[j].append(z)
+
+for i in range(len(o_atoms)):
+    data = np.array(o_atoms[i])
+    hist, bin_edges = np.histogram(data, 20)
+    dx = bin_edges[1] - bin_edges[0]
+    o = open("hist-o-%03d.dat"%i, "w")
+    for i in range(len(hist)):
+        o.write("%12.6f%8d\n"%(bin_edges[i] + 0.5*dx, hist[i]))
+    o.close()
+