import os import sys import shutil import copy import glob import csv import numpy as np from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, pi sys.path.append("c:/Programs/python/lib") sys.path.append("d:/Programs/python/lib") from tkfile import tkFile import tkre from tkutils import IsDir, IsFile, SplitFilePath from tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg from tksci.tksci import Reduce01, Round from tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3 from tkcrystal.tkcif import tkCIF, tkCIFData from tkcrystal.tkcrystal import tkCrystal from tkcrystal.tkvasp import tkVASP """ Convert VASP MD output files to XCrySDen AXSF file """ #================================ # global parameters #================================ debug = 0 # input path CAR_path = 'OUTCAR' xsffile = None #============================= # Treat argments #============================= def usage(): global CAR_path print("") print("Usage:") print(" (a) python {} outcar_path".format(sys.argv[0])) print(" ex: python {} {}".format(sys.argv[0], CAR_path)) def updatevars(): global CAR_path, xsffile CAR_path = getarg(1, CAR_path) dirname, basename, filebody, ext = SplitFilePath(CAR_path) xsffile = os.path.join(dirname, filebody + '.xsf') def WriteXSFAnimationFileHeader(outfile, nStep): xsf = tkFile(outfile, 'w') if not xsf: return None xsf.Write("ANIMSTEPS {}\n".format(nStep)) xsf.Write("CRYSTAL\n") xsf.Close() return 1 def WriteXSFFileFromCrystal(outfile, iCycle, cry): types = cry.AtomTypeList() xsf = tkFile(outfile, 'a') if not xsf: return None aij = cry.LatticeVectors() xsf.Write("PRIMVEC {}\n".format(iCycle)) xsf.Write(" %16.10f %16.10f %16.10f\n" % (aij[0][0], aij[0][1], aij[0][2])) xsf.Write(" %16.10f %16.10f %16.10f\n" % (aij[1][0], aij[1][1], aij[1][2])) xsf.Write(" %16.10f %16.10f %16.10f\n" % (aij[2][0], aij[2][1], aij[2][2])) xsf.Write("CONVVEC {}\n".format(iCycle)) xsf.Write(" %16.10f %16.10f %16.10f\n" % (aij[0][0], aij[0][1], aij[0][2])) xsf.Write(" %16.10f %16.10f %16.10f\n" % (aij[1][0], aij[1][1], aij[1][2])) xsf.Write(" %16.10f %16.10f %16.10f\n" % (aij[2][0], aij[2][1], aij[2][2])) types = cry.AtomTypeList() ntypes = len(types) sites = cry.ExpandedAtomSiteList() nsites = len(sites) xsf.Write("PRIMCOORD {}\n".format(iCycle)) xsf.Write(" %d 1\n" % (nsites)) #$Crystal2 != ''の時、重心の変化を調べる centx1, centy1, centz1 = 0.0, 0.0, 0.0 centx2, centy2, centz2 = 0.0, 0.0, 0.0 for i in range(nsites): atom = sites[i] name = atom.AtomNameOnly() itype = atom.iAtomType() AtomicNumber = types[itype].AtomicNumber() x, y, z = atom.Position(1) fx, fy, fz = atom.Force() if fx is None: fx, fy, fz = atom.Velocity() xc, yc, zc = cry.FractionalToCartesian(x, y, z) if fz is None: fx = fy = fz = 0.0 xsf.Write("%3d %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f\n" % (AtomicNumber, xc, yc, zc, fx, fy, fz)) xsf.Close() return 1 # my $xc = new XCrySDen; # $xc->WriteXSFAnimationFileHeader($XSFFile, $nStep); # $xc->WriteXSFFileFromCrystal($XSFFile, $iStep, $Crystal); def outcar2axsf(): global CAR_path, xsffile vasp = tkVASP() CAR_path = vasp.getdir(CAR_path) print("") INCAR_path = vasp.get_INCAR(CAR_path) POSCAR_path = vasp.get_POSCAR(CAR_path) CONTCAR_path = vasp.get_CONTCAR(CAR_path) OUTCAR_path = vasp.get_OUTCAR(CAR_path) print("CAR dir : ", CAR_path) print("INCAR : ", INCAR_path) print("POSCAR : ", POSCAR_path) print("CONTCAR : ", CONTCAR_path) print("OUTCAR : ", OUTCAR_path) print("AXSF file: ", xsffile) print("") print("*** Read initial crystal structure from [{}]".format(POSCAR_path)) initial_cry = vasp.read_poscar(POSCAR_path) # initial_cry.PrintInf("cell") a, b, c, alpha, beta, gamm = initial_cry.LatticeParameters() print("cell: {:12.8f} {:12.8f} {:12.8f} A {:10.6f} {:10.6f} {:10.6f}".format(a, b, c, alpha, beta, gamm)) print("") print("*** Read final crystal structure from [{}]".format(CONTCAR_path)) final_cry = vasp.read_poscar(CONTCAR_path) # final_cry.PrintInf("cell") a, b, c, alpha, beta, gamm = final_cry.LatticeParameters() print("cell: {:12.8f} {:12.8f} {:12.8f} A {:10.6f} {:10.6f} {:10.6f}".format(a, b, c, alpha, beta, gamm)) print("") print("Read [{}]".format(INCAR_path)) incar_inf = vasp.read_incar_inf(INCAR_path) if not incar_inf: terminate("Error: Can not read [{}]".format(INCAR_path), usage = usage) sample_name = incar_inf['SYSTEM'] if not sample_name: sample_name = 'hogehoge' dt = pfloat(incar_inf['POTIM']) print("") print("sample_name:", sample_name) print("POTIM: {} fs".format(dt)) initial_cifpath = sample_name + '-initial.cif' final_cifpath = sample_name + '-final.cif' historycsvfile = sample_name + '-history.csv' cif = tkCIFData() print("") print("Write to [{}]".format(initial_cifpath)) cif.CreateCIFFileFromCCrystal(initial_cry, initial_cifpath) print("Write to [{}]".format(final_cifpath)) cif = tkCIFData() cif.CreateCIFFileFromCCrystal(final_cry, final_cifpath) nstructures = vasp.read_n_crystalstructures(OUTCAR_path) print("# of crystal structures:", nstructures) if not WriteXSFAnimationFileHeader(xsffile, nstructures): terminate("", usage = usage) print("") print("Read crystal structures from [{}]".format(OUTCAR_path)) fp = tkFile(OUTCAR_path, 'r') if not fp: terminate("Error: Can not read [{}]".format(outcarpath), usage = usage) types = initial_cry.AtomTypeList() ntypes = len(types) sites = initial_cry.ExpandedAtomSiteList() nsites = len(sites) poslist0 = initial_cry.get_position_list(mode = 'all', IsReduce01 = True) atomnamelist = initial_cry.get_atom_name_list(mode = 'all', NameOnly = True) offset = make_matrix2(nsites, 3, defval = 0.0) nsitesfortypes = make_matrix1(ntypes, defval = 0) for i in range(nsites): itype = sites[i].iAtomType() nsitesfortypes[itype] += 1 print("# of atom sites for atom types = ", nsitesfortypes) count = 0 while 1: inf, cry = vasp.read_next_crystalstructure(fp, initial_cry, IsReduce01 = True) if not cry: break WriteXSFFileFromCrystal(xsffile, count, cry) count += 1 fp.Close() terminate(None, usage = usage) def main(): updatevars() print("") print("=============== Convert VASP output files to XCrySDen AXSF file ============") print("") outcar2axsf() if __name__ == "__main__": main()