新脚本&新功能

VASP.py

@@ -26,8 +26,7 @@ def read_total_atoms():
         return int(number)
 
 
-
-def grep_OUTCAR(command):
+def grep_outcar(command):
     pipe = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
     (content, error) = (pipe.stdout.readlines(), pipe.stderr.read())
     if content:
@@ -39,3 +38,115 @@ def grep_OUTCAR(command):
         print ""
         exit(0)
     return []
+
+
+# this function can only read file in VASP 5.0 or later
+def read_poscar(file_name):
+    space = re.compile(r'\s+')
+    with open(file_name) as input_file:
+        content = input_file.readlines()
+
+        lattice = float(content[1].strip())
+
+        basis = []
+        for i in xrange(2, 5):
+            line = space.split(content[i].strip())
+            basis.append([float(line[0]), float(line[1]), float(line[2])])
+
+        elements = space.split(content[5].strip())
+        num_atoms = space.split(content[6].strip())
+        for i in xrange(0, len(num_atoms)):
+            num_atoms[i] = int(num_atoms[i])
+
+        index = 0
+        if re.search(r'^[Ss]', content[7]) is None:
+            selectiveflag = ''
+            index = 7
+        else:
+            selectiveflag = content[7].strip()
+            index = 8
+
+        if re.search(r'^[Dd]', content[index]) is None:
+            coordinate_type = 'Cartesian'
+        else:
+            coordinate_type = 'Direct'
+
+        coordinates = []
+        selective = []
+        start = index + 1
+        end = start + sum(num_atoms)
+        if selectiveflag == '':
+            for i in xrange(start, end):
+                line = space.split(content[i].strip())
+                coordinates.append([float(line[0]), float(line[1]), float(line[2])])
+        else:
+            for i in xrange(start, end):
+                line = space.split(content[i].strip())
+                coordinates.append([float(line[0]), float(line[1]), float(line[2])])
+                selective.append([line[3], line[4], line[5]])
+
+    if re.search(r'^[Cc]', coordinate_type) is None:
+        coordinate_type = 'Cartesian'
+        coordinates = dirkar(basis, coordinates)
+
+    return lattice, basis, elements, num_atoms, selectiveflag, coordinate_type, coordinates, selective
+
+
+def write_poscar(file_name, lattice, basis, elements, num_atoms, selectiveflag, coordinate_type, coordinates, selective):
+    with open(file_name) as output_file:
+        description = ""
+        for atom in elements:
+            description += "%s  " % atom
+        output_file.write(description.rstrip() + '\n')
+        output_file.write("  %15.10f\n" % lattice)
+        for i in xrange(0, 3):
+            output_file.write("  %15.10f  %15.10f  %15.10f\n" % (basis[i][0], basis[i][1], basis[i][2]))
+        output_file.write(description.rstrip() + '\n')
+        num_atom = ''
+        for num in num_atoms:
+            num_atom += "%d  " % num
+        output_file.write(num_atom.rstrip() + '\n')
+        if selectiveflag != '':
+            output_file.write(selectiveflag + '\n')
+        output_file.write(coordinate_type + '\n')
+
+        if re.search(r'^[Dd]', coordinate_type):
+            coordinates = kardir(basis, coordinates)
+
+        for i in xrange(0, len(coordinates)):
+            output_file.write("%16.10f  %16.10f  %16.10f %s %s %s" % (coordinates[i][0], coordinates[i][1], coordinates[i][2],
+                                                                      selective[i][0], selective[i][1], selective[i][2]))
+
+
+def dirkar(basis, coordinates):
+    for i in xrange(0, len(coordinates)):
+        v1 = coordinates[i][0] * basis[0][0] + coordinates[i][1] * basis[1][0] + coordinates[i][2] * basis[2][0]
+        v2 = coordinates[i][0] * basis[0][1] + coordinates[i][1] * basis[1][1] + coordinates[i][2] * basis[2][1]
+        v3 = coordinates[i][0] * basis[0][2] + coordinates[i][1] * basis[1][2] + coordinates[i][2] * basis[2][2]
+        coordinates[i][0] = v1
+        coordinates[i][1] = v2
+        coordinates[i][2] = v3
+    return coordinates
+
+
+def kardir(basis, coordinates):
+    inverse = [[basis[1][1]*basis[2][2]-basis[2][1]*basis[1][2], basis[2][1]*basis[0][2]-basis[0][1]*basis[2][2], basis[0][1]*basis[1][2]-basis[1][1]*basis[0][2]],
+               [basis[2][0]*basis[1][2]-basis[1][0]*basis[2][2], basis[0][0]*basis[2][2]-basis[2][0]*basis[0][2], basis[1][0]*basis[0][2]-basis[0][0]*basis[1][2]],
+               [basis[1][0]*basis[2][1]-basis[2][0]*basis[1][1], basis[2][0]*basis[0][1]-basis[0][0]*basis[2][1], basis[0][0]*basis[1][1]-basis[1][0]*basis[0][1]]]
+    omega = basis[0][0]*basis[1][1]*basis[2][2] + basis[0][1]*basis[1][2]*basis[2][0] + basis[0][2]*basis[1][0]*basis[2][1] - \
+            basis[0][2]*basis[1][1]*basis[2][0] + basis[1][2]*basis[2][1]*basis[0][0] + basis[2][2]*basis[0][1]*basis[1][0]
+
+    inverse = [[inverse[0][0]/omega, inverse[0][1]/omega, inverse[0][2]/omega],
+               [inverse[1][0]/omega, inverse[1][1]/omega, inverse[1][2]/omega],
+               [inverse[2][0]/omega, inverse[2][1]/omega, inverse[2][2]/omega]]
+
+    for i in xrange(0, len(coordinates)):
+        v1 = coordinates[i][0] * inverse[0][0] + coordinates[i][1] * inverse[1][0] + coordinates[i][2] * inverse[2][0]
+        v2 = coordinates[i][0] * inverse[0][1] + coordinates[i][1] * inverse[1][1] + coordinates[i][2] * inverse[2][1]
+        v3 = coordinates[i][0] * inverse[0][2] + coordinates[i][1] * inverse[1][2] + coordinates[i][2] * inverse[2][2]
+
+        # move atoms to primative cell
+        coordinates[i][0] = v1 + 60 - int(v1 + 60)
+        coordinates[i][1] = v2 + 60 - int(v2 + 60)
+        coordinates[i][2] = v3 + 60 - int(v3 + 60)
+    return coordinates

excoor.py

@@ -5,14 +5,14 @@
 import sys
 import os
 from VASP import read_total_atoms
-from VASP import grep_OUTCAR
+from VASP import grep_outcar
 
 
 print ""
 number_of_atoms = read_total_atoms()
 
 string = "grep \"POSITION\" -A %d OUTCAR" % (int(number_of_atoms) + 15)
-content = grep_OUTCAR(string)
+content = grep_outcar(string)
 
 space = re.compile(r'\s+')
 position = re.compile(r'POSITION')

posdiff.py

@@ -0,0 +1,73 @@
+#!/bin/env python
+
+import sys
+from VASP import read_poscar
+
+threshold = 0.001
+if len(sys.argv) == 4:
+    threshold = float(sys.argv[3])
+elif len(sys.argv) != 3:
+    print ""
+    print "usage: poscar.py POSCAR1 POSCAR2 threshold"
+    print "try again!"
+    exit(0)
+
+lattice1, basis1, elements1, num_atoms1, selectiveflag1, coordinate_type1, coordinates1, selective1 = read_poscar(sys.argv[1])
+lattice2, basis2, elements2, num_atoms2, selectiveflag2, coordinate_type2, coordinates2, selective2 = read_poscar(sys.argv[2])
+
+
+if abs(lattice1 - lattice2) > threshold:
+    print ""
+    print "Lattices are different:"
+    print "%s: %10.5f" % (sys.argv[1], lattice1)
+    print "%s: %10.5f" % (sys.argv[2], lattice2)
+    print ""
+
+basis_diff = []
+for i in xrange(0, 3):
+    for j in xrange(0, 3):
+        if abs(basis1[i][j] - basis2[i][j]) > threshold:
+            basis_diff.append([i, j, basis1[i][j] - basis2[i][j]])
+
+if basis_diff:
+    print ""
+    print "Basis are differernt:"
+    print "    basis_id(within 3*3 array)            difference"
+    for i in basis_diff:
+        print "        (%4d,%4d)                  %14.10f" % (i[0] + 1, i[1] + 1, i[2])
+    print ""
+
+if elements1 != elements2:
+    print ""
+    print "Elements are differernt:"
+    print "%s: " % sys.argv[1], elements1
+    print "%s: " % sys.argv[2], elements2
+    print ""
+    exit(0)
+
+if num_atoms1 != num_atoms2:
+    print ""
+    print "Number of atoms are different:"
+    print "%s: " % sys.argv[1], num_atoms1
+    print "%s: " % sys.argv[2], num_atoms2
+    print ""
+    exit(0)
+
+coordinate_diff = []
+for i in xrange(0, len(coordinates1)):
+    for j in xrange(0, 3):
+        if abs(coordinates1[i][j] - coordinates2[i][j]) > threshold:
+            coordinate_diff.append([i, j, coordinates1[i][j] - coordinates2[i][j]])
+
+if coordinate_diff:
+    print ""
+    print "coordinates are different: "
+    print "      atoms_id              difference"
+    for i in coordinate_diff:
+        print "    (%4d,%4d)          %14.10f" % (i[0] + 1, i[1] + 1, i[2])
+    print ""
+
+
+
+
+