vasp中用于电荷密度计算的脚本（chgsum.pl）

#!/usr/bin/env python #encoding: utf-8 '''aselite is a striped down single file version of ase that retains the following features: atom and atoms objects, some of ase.io and some of ase.constraints.''' # Copyright 2008, 2009 CAMd # (see accompanying license files for details). from math import cos, sin, sqrt import warnings import numpy as np np.seterr(all='raise') import os import copy import sys import time from os.path import isfile import collections def read_any(filename): try: return read_vasp(filename) except: pass try: return read_xyz(filename) except: pass try: return read_con(filename) except: pass raise IOError("Could not read file %s." % filename) def write_jmol(filename, atoms, eigenvalues, eigenvectors): f_xyz = open(filename,'w') for i in range(len(eigenvectors)): mode = eigenvectors[:,i] mode.shape = (len(mode)/3,3) f_xyz.write("%i\n"%len(atoms)) f_xyz.write("%f\n"%eigenvalues[i]) for j,atom in enumerate(atoms): f_xyz.write("%s %f %f %f %f %f %f\n" % (atom.symbol, atom.position[0], atom.position[1], atom.position[2], mode[j,0], mode[j,1], mode[j,2])) f_xyz.close() def get_atomtypes(fname): """Given a file name, get the atomic symbols. The function can get this information from OUTCAR and POTCAR format files. The files can also be compressed with gzip or bzip2. """ atomtypes=[] if fname.find('.gz') != -1: import gzip f = gzip.open(fname) elif fname.find('.bz2') != -1: import bz2 f = bz2.BZ2File(fname) else: f = open(fname) for line in f: if line.find('TITEL') != -1: atomtypes.append(line.split()[3].split('_')[0].split('.')[0]) return atomtypes def atomtypes_outpot(posfname, numsyms): """Try to retreive chemical symbols from OUTCAR or POTCAR If getting atomtypes from the first line in POSCAR/CONTCAR fails, it might be possible to find the data in OUTCAR or POTCAR, if these files exist. posfname -- The filename of the POSCAR/CONTCAR file we're trying to read numsyms -- The number of symbols we must find """ import os.path as op import glob # First check files with exactly same name except POTCAR/OUTCAR instead # of POSCAR/CONTCAR. fnames = [posfname.replace('POSCAR', 'POTCAR').replace('CONTCAR', 'POTCAR')] fnames.append(posfname.replace('POSCAR', 'OUTCAR').replace('CONTCAR', 'OUTCAR')) # Try the same but with compressed files fsc = [] for fn in fnames: fsc.append(fn + '.gz') fsc.append(fn + '.bz2') for f in fsc: fnames.append(f) # Finally try anything with POTCAR or OUTCAR in the name vaspdir = op.dirname(posfname) fs = glob.glob(vaspdir + '*POTCAR*') for f in fs: fnames.append(f) fs = glob.glob(vaspdir + '*OUTCAR*') for f in fs: fnames.append(f) tried = [] files_in_dir = os.listdir('.') for fn in fnames: if fn in files_in_dir: tried.append(fn) at = get_atomtypes(fn) if len(at) == numsyms: return at raise IOError('Could not determine chemical symbols. Tried files ' + str(tried)) def get_atomtypes_from_formula(formula): """Return atom types from chemical formula (optionally prepended with and underscore). """ symbols = string2symbols(formula.split('_')[0]) atomtypes = [symbols[0]] for s in symbols[1:]: if s != atomtypes[-1]: atomtypes.append(s) return atomtypes def read_vasp(filename='CONTCAR'): """Import POSCAR/CONTCAR type file. Reads unitcell, atom positions and constraints from the POSCAR/CONTCAR file and tries to read atom types from POSCAR/CONTCAR header, if this fails the atom types are read from OUTCAR or POTCAR file. """ if isinstance(filename, str): f = open(filename) else: # Assume it's a file-like object f = filename # First line should contain the atom symbols , eg. "Ag Ge" in # the same order # as later in the file (and POTCAR for the full vasp run) atomtypes = f.readline().split() # Sometimes the first line in POSCAR/CONTCAR is of the form # "CoP3_In-3.pos". Check for this case and extract atom types if len(atomtypes) == 1 and '_' in atomtypes[0]: atomtypes = get_atomtypes_from_formula(atomtypes[0]) lattice_constant = float(f.readline().split()[0]) # Now the lattice vectors a = [] for ii in range(3): s = f.readline().split() floatvect = float(s[0]), float(s[1]), float(s[2]) a.append(floatvect) basis_vectors = np.array(a) * lattice_constant # Number of atoms. Again this must be in the same order as # in the first line # or in the POTCAR or OUTCAR file atom_symbols = [] numofatoms = f.readline().split() #vasp5.1 has an additional line which gives the atom types #the following try statement skips this line try: int(numofatoms[0]) except ValueError: numofatoms = f.readline().split() # check for comments in numofatoms line and get rid of them if necessary commentcheck = np.array(['!' in s for s in numofatoms]) if commentcheck.any(): # only keep the elements up to the first including a '!': numofatoms = numofatoms[:np.arange(len(numofatoms))[commentcheck][0]] numsyms = len(numofatoms) if len(atomtypes) < numsyms: # First line in POSCAR/CONTCAR didn't contain enough symbols. atomtypes = atomtypes_outpot(f.name, numsyms) else: try: for atype in atomtypes[:numsyms]: if not atype in chemical_symbols: raise KeyError except KeyError: atomtypes = atomtypes_outpot(f.name, numsyms) for i, num in enumerate(numofatoms): numofatoms[i] = int(num) [atom_symbols.append(atomtypes[i]) for na in range(numofatoms[i])] # Check if Selective dynamics is switched on sdyn = f.readline() selective_dynamics = sdyn[0].lower() == "s" # Check if atom coordinates are cartesian or direct if selective_dynamics: ac_type = f.readline() else: ac_type = sdyn cartesian = ac_type[0].lower() == "c" or ac_type[0].lower() == "k" tot_natoms = sum(numofatoms) atoms_pos = np.empty((tot_natoms, 3)) if selective_dynamics: selective_flags = np.empty((tot_natoms, 3), dtype=bool) for atom in range(tot_natoms): ac = f.readline().split() atoms_pos[atom] = (float(ac[0]), float(ac[1]), float(ac[2])) if selective_dynamics: curflag = [] for flag in ac[3:6]: curflag.append(flag == 'F') selective_flags[atom] = curflag # Done with all reading if type(filename) == str: f.close() if cartesian: atoms_pos *= lattice_constant atoms = Atoms(symbols = atom_symbols, cell = basis_vectors, pbc = True) if cartesian: atoms.set_positions(atoms_pos) else: atoms.set_scaled_positions(atoms_pos) if selective_dynamics: constraints = [] indices = [] for ind, sflags in enumerate(selective_flags): if sflags.any() and not sflags.all(): constraints.append(FixScaled(atoms.get_cell(), ind, sflags)) elif sflags.all(): indices.append(ind) if indices: constraints.append(FixAtoms(indices)) if constraints: atoms.set_constraint(constraints) atoms.format = 'vasp' return atoms def read_vasp_out(filename='OUTCAR',index = 'all'): """Import OUTCAR type file. Reads unitcell, atom positions, energi