#!/usr/bin/env python # slices wien2k cell in x direction -> each slice is a cut parallel to y-z-plane # calling program: programname nx ny nz x_start x_end # # creates 2d_slices with ny times nz points # spacing between different slices is: length of unit cell in x-direction divided by nx # x_start and x_end are defining the x-range that is sampled, including x_start slice and x_end slice # - this can be used do include mirror symmetrie along sampling axis -> you only have to sample half of the unit cell # - that is important only if you want to call this script several times for the same sturcture to run on several cores # # asking for 200 slides each 10 by 10 pixels: programmname 200 10 10 1 200 # be careful not to create identical slices: if you ask for 200 10 10 0 200 slice 0 and 200 are equal # for parallel execution ask: # programmname 200 10 10 1 50& # programmname 200 10 10 51 100& # programmname 200 10 10 101 150& # programmname 200 10 10 151 200& # Theres the possibility to shift each slice in y and z direction by defining shift vectors (in line ~180) # if you want to sample y from 0.01:1.01 define dy=1 dy_den=100 # this can be important if you want to avoid critical potential values because the sampling exactly hits a core import sys, os, string, array import time class Lapw5: def __init__(self, head = None): self.tails = ["clmsum", "clmval", "vcoul", "vtotal", "r2v"] if head: self.head = head else: self.head = os.path.basename(os.getcwd()) self.shells = [3, 3, 3] self.val = 2 #Wass soll ausgegeben werden? clmsum, clmval, vcoul.. self.in5file = head + ".in5c" def make_ctrl(self, origin, xend, yend, nx, ny): output = open(self.in5file, "w") output.write("%d %d %d %d\n" % (origin[0], origin[1], origin[2], origin[3])) output.write("%d %d %d %d\n" % (xend[0], xend[1], xend[2], xend[3])) output.write("%d %d %d %d\n" % (yend[0], yend[1], yend[2], yend[3])) output.write("%d %d %d\n" % (self.shells[0], self.shells[1], self.shells[2])) output.write("%d %d\n" % (nx, ny)) output.write("RHO \n") if self.val != 0: output.write("ATU VAL NODEBUG\n") else: output.write("ATU TOT NODEBUG\n") output.write("NONORTHO") output.close() def _make_def(self): deffile = "lapw5.def" os.system("x lapw5c -d") if self.val != 1: input = open(deffile, "r") lines = input.readlines() input.close() output = open(deffile, "w") for l in lines: l = string.replace(l, "clmval", self.tails[self.val]) #l = string.replace(l, ".in5", ".in5c") output.write(l) output.close() def run(self): self._make_def() os.system("lapw5c lapw5.def") def _string2angle(s): s = string.strip(s) if s: angle = string.atof(s) if angle == 0.0: angle = 90.0 return angle else: return 90.0 def get_lattice(struct): input = open(struct, "r") line = input.readline() line = input.readline() lat_type = string.strip(line[:4]) line = input.readline() line = input.readline() abc = map(string.atof, [line[0:10],line[10:20],line[20:30]]) angle = map(_string2angle, [line[30:40],line[40:50],line[50:60]]) return abc, angle def collect_charge(head, num): rhofile = head + ".rho" rho = array.array("d") for i in range(num): suffix = ".%04d" % i rhof = rhofile + suffix input = open(rhof, "r") line = input.readline() while 1: line = input.readline() if not line: break for word in string.split(line): rho.append(float(word)) return rho def erase_files(head, num): rhofile = head + ".rho" in5file = head + ".in5c" for i in range(num): suffix = ".%04d" % i os.unlink(rhofile + suffix) os.unlink(in5file + suffix) if __name__ == "__main__": import getopt def print_usage(): usage = '''usage: venus.py [-h] [-v] [-p] nx ny nz start_z end_z -h print help. -v verbose mode. -p preserve intermediate files\n''' sys.stderr.write(usage) try: options, resargs = getopt.getopt(sys.argv[1:], 'hvp') except getopt.error, s: sys.stderr.write(str(s)+"\n") print_usage() sys.exit(1) if len(resargs) < 5: print_usage() sys.exit(1) verbose = 0 preserve = 0 for option in options: if option[0] == "-h": print_usage() sys.exit() elif option[0] == "-v": verbose = 1 elif option[0] == "-p": preserve = 1 nx, ny, nz = map(int, resargs[0:3]) x_start, x_end = map(int, resargs[3:5]) print "using lapw5c for creating slices" print "slicing cell in x-direction -> each slice is parallel to y-z-plane" print "mesh:", nx, ny, nz print "x_start, x_end:", x_start, x_end head = os.path.basename(os.getcwd()) lapw5 = Lapw5(head) x=x_start while x <= x_end : if verbose: print " --> ", x, " ", time.strftime('%d_%m_%Y %H:%M:%S') # no shift # origin = [x , 0*nx , 0*nx , nx] # xend = [x , 1*nx , 0*nx , nx] # yend = [x , 0*nx , 1*nx , nx] #shift in y direction and z direction # if you want to shift y by 0.06 of unit-cell-vector: # 0.06=6/100 -> dy=6 dy_den=100 # if you do not want to shift set: dy=0 dy_den=1 # #I'm not sure if a limit is reached for very small shifts when the total denominater becomes very large #there may be problems on lapw5/lapw5c because this big number may not fit into their datatype #(eg if lapw5 is reading the denominator from case.in5 as integer number dy=0 dy_den=1 dz=0 dz_den=1 dx=0 dx_den=1 origin = [(x*dx_den+dx*nx)*dy_den*dz_den, nx*dx_den*(0+dy)*dz_den, nx*dx_den*dy_den*(0+dz), nx*dx_den*dy_den*dz_den] xend = [(x*dx_den+dx*nx)*dy_den*dz_den, nx*dx_den*(dy_den+dy)*dz_den, nx*dx_den*dy_den*(0+dz), nx*dx_den*dy_den*dz_den] yend = [(x*dx_den+dx*nx)*dy_den*dz_den, nx*dx_den*(0+dy)*dz_den, nx*dx_den*dy_den*(dz_den+dz), nx*dx_den*dy_den*dz_den] lapw5.make_ctrl(origin, xend, yend, ny, nz) lapw5.run() suffix = ".%04d" % x os.rename(head + ".in5c", head + ".in5c" + suffix) rhofile = head + ".rho" x +=1 if os.path.exists(rhofile) and os.path.getsize(rhofile) > 0: os.rename(rhofile, rhofile + suffix) else: print rhofile, "does not exist or empty???" sys.exit(1) alat, angle = get_lattice(head + ".struct") output = open(head + ".rho3d" , "w") output.write("cell\n") output.write("%f %f %f\n" % (alat[1], alat[2], alat[0])) output.write("%f %f %f\n" % ( angle[1], angle[2], angle[0])) output.write("%d %d %d " % ( ny, nz,nx)) output.write("%f %f %f\n" % ( alat[1], alat[2], alat[0])) #rho = collect_charge(head, nz/2) # i = 0 # for z in rho: # output.write("%15.8e " % z) # i = i + 1 # if (i % 5) == 0: # output.write("\n") # output.close() #if not preserve: # erase_files(head, nz-1)