#!/usr/bin/python ######################################################################################################################## # # ############ # # OVERVIEW # (and what this Python module does...) # ############ # # Module written by Mark Wijzenbroek in October/November 2014. Tested with VASP 5.2.12, should also work with later # versions. # # Read a VASP XML file and compute/obtain positions, velocities and forces and convert them into a readable format. # NOTE: This will probably only work exactly if VASP is linked with stepver.o. # # It should be kept in mind that VASP applies a leap frog integrator (stepver.F) and as such, the velocities that VASP # wants are half a time step off (earlier) compared to the positions. In fact, graphically, what VASP does is this: # # STEP /1\ /2\ /3\ /n\ ... # t/dt 0 1 2 n-1 # v o o o o ... # x o o o o ... # a o o o o ... # # ############################# # # THE POSCAR/CONTCAR FORMAT # (why you need to watch out, and some pesky little details the VASP manual doesn't # ############################# tell you, but are important nonetheless...) # # A POSCAR file consists of four regions: the header, positions, velocities and predictor corrector coordinates. # ( http://cms.mpi.univie.ac.at/vasp/guide/node59.html ) The header and positions are mandatory. # For a new calculation usually only the first two (header and positions) or three regions (also velocities) are # provided. If no velocities are provided, they are initialized from a Maxwell-Boltzmann distribution at T=TEBEG # ( http://cms.mpi.univie.ac.at/vasp/guide/node117.html ). If velocities are provided, they are given at t=-dt/2, # not at t=0, due to the used propagator. # # For a continuation calculation CONTCAR is copied to POSCAR. VASP then uses a rather strange format in which the # predictor corrector coordinates *are* included. The positions given in CONTCAR are of the **last completed step # of the MD**. In other words, at t=(NSW-1)*dt. The velocities however are given at t=(NSW-1/2)*dt. As such, this # is by itself not a valid input for VASP. (The velocities are now given half a timestep after the positions instead # of the other way around, and the positions are wrong.) To solve this, VASP uses the predictor corrector coordinates # to start from the right initial geometry. These predictor corrector coordinates consist again of multiple parts: # the header, the positions, the velocities and another array (unused in the leap frog propagator). The header consists # of two lines ("1" to enable the predictor corrector coordinates and an array of 4 floats related to SMASS and the # Nose algorithm if SMASS>0, see also http://cms.mpi.univie.ac.at/vasp/guide/node118.html ). After that the real # positions (at t=NSW*dt) are given and they overrule(!) the ones given in the regular positions block. The velocities # given in the block below are scaled with respect to timestep and basis vectors, are not used and correspond to # t=(NSW-1)*dt. # # ############################## # # VASP'S LEAP FROG ALGORITHM # (well, idem.) # ############################## # # For the first time step, the positions are given at t=0 and the velocities are given at t=-dt/2. The forces and # potential energy are then calculated for the geometry at t=0. VASP then propagates the velocities, using the forces # from the geometry corresponding to t=0, one whole time step (to t=dt/2). After that, with the newly calculated # velocities, the positions are updated one whole timestep to t=dt. # # If velocity rescaling is requested it is done immediately at t=0 (step 1), and then every t=NBLOCK*dt, and as such # your computed velocities may then differ at these steps from the ones given in POSCAR. The rescaling is done at # whole time steps, by doing a half update of the velocities and then rescaling these to yield the correct temperature. # The computed scaling factor at t=0 is used to rescale velocities at t=-dt/2 and these are then used to propagate the # whole time step (the half step propagation is completely discarded). The printed kinetic energies in VASP correspond # to whole timestep velocities, for all steps. # # This algorithm is in total done NSW times. Due to the algorithm, the way to properly compute the velocities given the # positions and forces (accelerations) is by: # # v_{n+1/2} = ( x_{n+1} - x_{n} ) / dt, therefore # v_{n-1/2} = ( x_{n+1} - x_{n} ) / dt - ( a_{n} * dt ). # # ############ # # SEE ALSO # # ############ # # 1. stepver.F in VASP source code # 2. http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet # 3. http://en.wikipedia.org/wiki/Leapfrog_integration # ######################################################################################################################## import libxml2 import sys import os import math import numpy as np def DirectToCartesian( x, basis ): return [ x[0] * basis[0][0] + x[1] * basis[1][0] + x[2] * basis[2][0], x[0] * basis[0][1] + x[1] * basis[1][1] + x[2] * basis[2][1], x[0] * basis[0][2] + x[1] * basis[1][2] + x[2] * basis[2][2] ] # Find the nearest mirror image for all atoms in list next. def FindMirrorImage( cur, next, basis ): new = [] for i in range( 0, len( cur ) ): shortest = None for j in [ -1, 0, 1 ]: for k in [ -1, 0, 1 ]: for l in [ -1, 0, 1 ]: temp = [ next[i][0] + j * basis[0][0] + k * basis[1][0] + l * basis[2][0], next[i][1] + j * basis[0][1] + k * basis[1][1] + l * basis[2][1], next[i][2] + j * basis[0][2] + k * basis[1][2] + l * basis[2][2] ] r = ( temp[0] - cur[i][0] )**2 + ( temp[1] - cur[i][1] )**2 + ( temp[2] - cur[i][2] )**2 if shortest is None or r < shortest[0]: shortest = [ r, j, k, l ] new.append( [ next[i][0] + shortest[1] * basis[0][0] + shortest[2] * basis[1][0] + shortest[3] * basis[2][0], next[i][1] + shortest[1] * basis[0][1] + shortest[2] * basis[1][1] + shortest[3] * basis[2][1], next[i][2] + shortest[1] * basis[0][2] + shortest[2] * basis[1][2] + shortest[3] * basis[2][2] ] ) return new def FindMirrorImage2( next, curdirect, nextdirect, basis ): new = [] for i in range( 0, len( curdirect ) ): shortest = [ 0.0, 0, 0, 0 ] if curdirect[i][0] < 0.25 and nextdirect[i][0] > 0.75: shortest[1] = -1 elif curdirect[i][0] > 0.75 and nextdirect[i][0] < 0.25: shortest[1] = 1 if curdirect[i][1] < 0.25 and nextdirect[i][1] > 0.75: shortest[2] = -1 elif curdirect[i][1] > 0.75 and nextdirect[i][1] < 0.25: shortest[2] = 1 if curdirect[i][2] < 0.25 and nextdirect[i][2] > 0.75: shortest[3] = -1 elif curdirect[i][2] > 0.75 and nextdirect[i][2] < 0.25: shortest[3] = 1 new.append( [ next[i][0] + shortest[1] * basis[0][0] + shortest[2] * basis[1][0] + shortest[3] * basis[2][0], next[i][1] + shortest[1] * basis[0][1] + shortest[2] * basis[1][1] + shortest[3] * basis[2][1], next[i][2] + shortest[1] * basis[0][2] + shortest[2] * basis[1][2] + shortest[3] * basis[2][2] ] ) return new class VASPXMLFile: xmldoc = None filename = None atoms = [] atomtypes = [] atomcounts = [] atommasses = [] energies_e0 = [] energies_f = [] energies_ekin = [] geometries = [] geometries_cartesian = [] velocities = [] forces = [] accelerations = [] # Note: this array actually contains acceleration * dt basis = [] selective = [] # The VASP " x x x " string selective_int = [] # With bools nsw = None potim = None def __init__( self, filename ): self.atoms = [] self.atomtypes = [] self.atomcounts = [] self.atommasses = [] self.energies_e0 = [] self.energies_f = [] self.energies_ekin = [] self.geometries = [] self.geometries_cartesian = [] self.velocities = [] self.forces = [] self.accelerations = [] self.basis = [] self.selective = [] self.selective_int = [] self.filename = filename self.xmldoc = libxml2.parseFile( filename ) calcs = self.xmldoc.xpathEval( "/modeling/calculation" ) #structs = doc.xpathEval("//calculation/structure") #[not(@name)]") #energies = doc.xpathEval("//calculation/scstep[last()]/energy/i[@name='e_0_energy']") self.nsw = int( self.xmldoc.xpathEval( "/modeling/incar/i[@name='NSW']" )[0].content ) self.potim = float( self.xmldoc.xpathEval( "/modeling/incar/i[@name='POTIM']" )[0].content ) #atomtypes = [] #temp = doc.xpathEval( "/modeling/atominfo/array[@name='atomtypes']/set/rc" ) temp = self.xmldoc.xpathEval( "/modeling/atominfo/array[@name='atoms']/set/rc" ) for i in temp: atomname = i.xpathEval( "c[1]" )[0].content.strip() self.atoms.append( atomname ) temp = self.xmldoc.xpathEval( "/modeling/atominfo/array[@name='atomtypes']/set/rc" ) for i in temp: atomtype = i.xpathEval( "c[2]" )[0].content.strip() atomcount = int( i.xpathEval( "c[1]" )[0].content.strip() ) atommass = float( i.xpathEval( "c[3]" )[0].content.strip() ) self.atomtypes.append( atomtype ) self.atomcounts.append( atomcount ) self.atommasses.append( atommass ) # because VASP does stupid things, as usual, check whether we manually specified a different mass temp = self.xmldoc.xpathEval( "/modeling/incar/v[@name='POMASS']" ) if len(temp) > 0: self.atommasses = [ float(x) for x in temp[0].content.strip().split() ] temp = self.xmldoc.xpathEval( "/modeling/structure[@name='initialpos']/varray[@name='selective']/v" ) for i in temp: selarr = [ x.lower() == "t" for x in i.content.strip().split() ] self.selective.append( i.content.strip() ) self.selective_int.append( selarr ) step = 0 for calc in calcs: f = float( calc.xpathEval( "scstep[last()]/energy/i[@name='e_fr_energy']" )[0].content ) e0 = float( calc.xpathEval( "scstep[last()]/energy/i[@name='e_0_energy']" )[0].content ) ekin = float( calc.xpathEval( "energy/i[@name='kinetic']" )[0].content ) etot = float( calc.xpathEval( "energy/i[@name='total']" )[0].content ) self.energies_e0.append( e0 ) self.energies_f.append( f ) self.energies_ekin.append( ekin ) # energies.append( [ f, e0, ekin, etot ] ) # print step * potim, f, e0, ekin, etot basis = [ [ float( y ) for y in x.content.split() ] for x in calc.xpathEval( "structure/crystal/varray[@name='basis']/v" ) ] positions = [ [ float( y ) for y in x.content.split() ] for x in calc.xpathEval( "structure/varray[@name='positions']/v" ) ] forces = [ [ float( y ) for y in x.content.split() ] for x in calc.xpathEval( "varray[@name='forces']/v" ) ] self.basis.append( basis ) self.geometries.append( positions ) self.forces.append( forces ) positions_cartesian = [] for x in positions: positions_cartesian.append( DirectToCartesian( x, basis ) ) #[ x[0] * basis[0][0] + x[1] * basis[1][0] + x[2] * basis[2][0], # x[0] * basis[0][1] + x[1] * basis[1][1] + x[2] * basis[2][1], # x[0] * basis[0][2] + x[1] * basis[1][2] + x[2] * basis[2][2] ] ) # sys.exit() self.geometries_cartesian.append( positions_cartesian ) step += 1 # now, we do the first structure of the next run, which we need to compute velocities basis = [ [ float( y ) for y in x.content.split() ] for x in self.xmldoc.xpathEval( "/modeling/structure[@name='finalpos']/crystal/varray[@name='basis']/v" ) ] finalstructure = [ [ float( y ) for y in x.content.split() ] for x in self.xmldoc.xpathEval( "/modeling/structure[@name='finalpos']/varray[@name='positions']/v" ) ] self.basis.append( basis ) self.geometries.append( finalstructure ) finalstructure_cartesian = [] for x in finalstructure: finalstructure_cartesian.append( [ x[0] * basis[0][0] + x[1] * basis[1][0] + x[2] * basis[2][0], x[0] * basis[0][1] + x[1] * basis[1][1] + x[2] * basis[2][1], x[0] * basis[0][2] + x[1] * basis[1][2] + x[2] * basis[2][2] ] ) self.geometries_cartesian.append( finalstructure_cartesian ) for i in range( 0, len( self.forces ) ): accel = [ [ self.forces[i][j][k] * 1.60217733e-9 / ( 1.6605402e-27 * self.GetMass( self.atoms[j] ) ) * self.potim * 1.0e-15 * 1.0e-5 for k in range( 0, 3 ) ] for j in range( 0, len(self.forces[i]) ) ] for i in range( 0, len( accel ) ): for j in range( 0, 3 ): if not self.selective_int[i][j]: accel[i][j] = 0.0 self.accelerations.append( accel ) # finally, calculate the velocities based on the positions we have for i in range( 0, len( calcs ) ): velocities_temp = [] for j in range( 0, len( self.geometries_cartesian[i] ) ): next = self.geometries_cartesian[i+1] cur = self.geometries_cartesian[i] nextdirect = self.geometries[i+1] curdirect = self.geometries[i] # now, we find the mirror image which is closest by # new = next # new = FindMirrorImage( cur, next, basis ) new = FindMirrorImage2( next, curdirect, nextdirect, basis ) # velocities_temp.append( [ ( new[j][0] - cur[j][0] ) / self.potim - 0.5 * self.forces[i][j][0] * 1.60217733e-9 / ( 1.6605402e-27 * self.GetMass( self.atoms[j] ) ) * self.potim * 1.0e-15 * 1.0e-5, # ( new[j][1] - cur[j][1] ) / self.potim - 0.5 * self.forces[i][j][1] * 1.60217733e-9 / ( 1.6605402e-27 * self.GetMass( self.atoms[j] ) ) * self.potim * 1.0e-15 * 1.0e-5, # ( new[j][2] - cur[j][2] ) / self.potim - 0.5 * self.forces[i][j][2] * 1.60217733e-9 / ( 1.6605402e-27 * self.GetMass( self.atoms[j] ) ) * self.potim * 1.0e-15 * 1.0e-5 ] ) # velocities_temp.append( [ ( new[j][0] - cur[j][0] ) / self.potim - self.forces[i][j][0] * 1.60217733e-9 / ( 1.6605402e-27 * self.GetMass( self.atoms[j] ) ) * self.potim * 1.0e-15 * 1.0e-5, # ( new[j][1] - cur[j][1] ) / self.potim - self.forces[i][j][1] * 1.60217733e-9 / ( 1.6605402e-27 * self.GetMass( self.atoms[j] ) ) * self.potim * 1.0e-15 * 1.0e-5, # ( new[j][2] - cur[j][2] ) / self.potim - self.forces[i][j][2] * 1.60217733e-9 / ( 1.6605402e-27 * self.GetMass( self.atoms[j] ) ) * self.potim * 1.0e-15 * 1.0e-5 ] ) # THE VELOCITIES AND THE POSITIONS BELONG NOW TO THE SAME TIME STEP, THAT IS x(t) and v(t) velocities_temp.append( [ ( new[j][0] - cur[j][0] ) / self.potim - 0.5 * self.accelerations[i][j][0], ( new[j][1] - cur[j][1] ) / self.potim - 0.5 * self.accelerations[i][j][1], ( new[j][2] - cur[j][2] ) / self.potim - 0.5 * self.accelerations[i][j][2] ] ) # THE VELOCITIES AND THE POSITIONS BELONG NOW TO THE SAME TIME STEP, THAT IS x(t) and v(t - 0.25fs) # velocities_temp.append( [ ( new[j][0] - cur[j][0] ) / self.potim - 0.75 * self.accelerations[i][j][0], # ( new[j][1] - cur[j][1] ) / self.potim - 0.75 * self.accelerations[i][j][1], # ( new[j][2] - cur[j][2] ) / self.potim - 0.75 * self.accelerations[i][j][2] ] ) self.velocities.append( velocities_temp ) self.xmldoc.freeDoc() def GetMass( self, symbol ): for i in range( 0, len( self.atommasses ) ): if self.atomtypes[i] == symbol: return self.atommasses[i] raise Exception( "Could not find mass for element " + symbol ) def ComputeTemperature( self, atommask=None, minushalf=False ): if atommask is None: atommask = range( 0, len( self.atoms ) ) DOF = sum( [ x.count(True) for idx, x in enumerate( self.selective_int ) if idx in atommask ] ) # print DOF # return AmuToKg = 1.6605402e-27 BoltzmannConstant = 1.38065e-23 eVToJoule = 1.60217733e-19 #1.60218e-19 DiracConstant = 1.054571e-34 list = [] list2 = [] for i in range( 0, self.nsw ): Ekin = 0.0 for j in atommask: add = [ 0.0, 0.0, 0.0 ] if minushalf == False: add = [ 0.5 * self.accelerations[i][j][k] for k in range( 0, 3 ) ] atomvels = [ self.velocities[i][j][k] + add[k] for k in range( 0, 3 ) ] Ekin += 0.5 * AmuToKg * self.GetMass( self.atoms[j] ) * ( ( atomvels[0] * 1.0e5 )**2 + ( atomvels[1] * 1.0e5 )**2 + ( atomvels[2] * 1.0e5 )**2 ) # print Ekin / eVToJoule list.append( Ekin / eVToJoule ) list2.append( ( Ekin / eVToJoule ) / ( 0.5 * DOF * 8.617343e-5 ) ) return list, list2 def GetPoscar( self, step ): arr = [] # header arr.append( self.filename + " step " + str(step) ) arr.append( "1.0" ) # basis vectors arr.append( "{0[0]:20.15f} {0[1]:20.15f} {0[2]:20.15f}".format( self.basis[step][0] ) ) arr.append( "{0[0]:20.15f} {0[1]:20.15f} {0[2]:20.15f}".format( self.basis[step][1] ) ) arr.append( "{0[0]:20.15f} {0[1]:20.15f} {0[2]:20.15f}".format( self.basis[step][2] ) ) # arr.append( " ".join( [ str(x) for x in self.basis[step][0] ] ) ) # arr.append( " ".join( [ str(x) for x in self.basis[step][1] ] ) ) # arr.append( " ".join( [ str(x) for x in self.basis[step][2] ] ) ) # atom types arr.append( " ".join( self.atomtypes ) ) arr.append( " ".join( [ str(x) for x in self.atomcounts ] ) ) # positions arr.append( "Selective dynamics" ) arr.append( "Cartesian" ) for i in range( 0, len( self.geometries_cartesian[step] ) ): arr.append( "{0[0]:20.15f} {0[1]:20.15f} {0[2]:20.15f}".format( self.geometries_cartesian[step][i] ) + " " + self.selective[i] ) # arr.append( " ".join( [ str(x) for x in self.geometries_cartesian[step][i] ] ) + " " + self.selective[i] ) # velocities arr.append( " " ) for i in range( 0, len( self.geometries_cartesian[step] ) ): arr.append( "{0[0]:20.15f} {0[1]:20.15f} {0[2]:20.15f}".format( self.velocities[step][i] ) ) # arr.append( " ".join( [ str(x) for x in self.velocities[step][i] ] ) ) return "\n".join( arr ) def GetForces( self, step ): arr = [] for i in range( 0, len( self.forces[step] ) ): arr.append( " ".join( [ str(x) for x in self.forces[step][i] ] ) ) return "\n".join( arr ) def GetMoleculeInfo( self, step ): atom1 = self.geometries_cartesian[step][-2] atom2 = self.geometries_cartesian[step][-1] # pos = FindMirrorImage2( [ atom2 ], [ self.geometries[step][-2] ], [ self.geometries[step][-1] ], self.basis[step] ) pos = FindMirrorImage( [atom1], [atom2], self.basis[step] ) atom2 = pos[0] r = math.sqrt( ( atom1[0] - atom2[0] )**2 + ( atom1[1] - atom2[1] )**2 + ( atom1[2] - atom2[2] )**2 ) theta = math.acos( ( atom2[2] - atom1[2] ) / r ) phi = math.atan2( ( atom2[1] - atom1[1] ) / ( r * math.sin( theta ) ), ( atom2[0] - atom1[0] ) / ( r * math.sin( theta ) ) ) X = ( atom1[0] + atom2[0] ) / 2.0 Y = ( atom1[1] + atom2[1] ) / 2.0 Z = ( atom1[2] + atom2[2] ) / 2.0 AmuToKg = 1.6605402e-27 eVToJoule = 1.60217733e-19 #1.60218e-19 vel1 = [ self.velocities[step][-2][k] + 0.5 * self.accelerations[step][-2][k] for k in range( 0, 3 ) ] vel2 = [ self.velocities[step][-1][k] + 0.5 * self.accelerations[step][-1][k] for k in range( 0, 3 ) ] mass1 = self.GetMass( self.atoms[-2] ) * AmuToKg mass2 = self.GetMass( self.atoms[-1] ) * AmuToKg Ekin1 = [ 0.5 * mass1 * ( vel1[k] * 1.0e5 )**2 / eVToJoule for k in range( 0, 3 ) ] Ekin2 = [ 0.5 * mass2 * ( vel2[k] * 1.0e5 )**2 / eVToJoule for k in range( 0, 3 ) ] # Ekin1 = 0.5 * mass1 * ( ( vel1[0] * 1.0e5 )**2 + ( vel1[1] * 1.0e5 )**2 + ( vel1[2] * 1.0e5 )**2 ) / eVToJoule # Ekin2 = 0.5 * mass2 * ( ( vel2[0] * 1.0e5 )**2 + ( vel2[1] * 1.0e5 )**2 + ( vel2[2] * 1.0e5 )**2 ) / eVToJoule angm = [ ( atom1[2] - ( atom1[2] + atom2[2] ) / 2.0 ) * mass1 * vel1[0] * 1.0e5 - ( atom1[0] - ( atom1[0] + atom2[0] ) / 2.0 ) * mass1 * vel1[2] * 1.0e5 + ( atom2[2] - ( atom1[2] + atom2[2] ) / 2.0 ) * mass2 * vel2[0] * 1.0e5 - ( atom2[0] - ( atom1[0] + atom2[0] ) / 2.0 ) * mass2 * vel2[2] * 1.0e5, ( atom1[0] - ( atom1[0] + atom2[0] ) / 2.0 ) * mass1 * vel1[1] * 1.0e5 - ( atom1[1] - ( atom1[1] + atom2[1] ) / 2.0 ) * mass1 * vel1[0] * 1.0e5 + ( atom2[0] - ( atom1[0] + atom2[0] ) / 2.0 ) * mass2 * vel2[1] * 1.0e5 - ( atom2[1] - ( atom1[1] + atom2[1] ) / 2.0 ) * mass2 * vel2[0] * 1.0e5, ( atom1[1] - ( atom1[1] + atom2[1] ) / 2.0 ) * mass1 * vel1[2] * 1.0e5 - ( atom1[2] - ( atom1[2] + atom2[2] ) / 2.0 ) * mass1 * vel1[1] * 1.0e5 + ( atom2[1] - ( atom1[1] + atom2[1] ) / 2.0 ) * mass2 * vel2[2] * 1.0e5 - ( atom2[2] - ( atom1[2] + atom2[2] ) / 2.0 ) * mass2 * vel2[1] * 1.0e5 ] redmass = ( mass1 * mass2 ) / ( mass1 + mass2 ) Ekinrot = ( angm[0]**2 + angm[1]**2 + angm[2]**2 ) / ( 2 * redmass * r**2 ) / eVToJoule Ekintrans = [ 0.5 * ( mass1 + mass2 ) * ( ( vel1[k] + vel2[k] ) / 2.0 * 1.0e5 )**2 / eVToJoule for k in range( 0, 3 ) ] # Ekintrans = 0.5 * ( mass1 + mass2 ) * ( ( vel1[0] + vel2[0] ) / 2.0 * 1.0e5 )**2 + ( ( vel1[1] + vel2[1] ) / 2.0 * 1.0e5 )**2 + ( ( vel1[2] + vel2[2] ) / 2.0 * 1.0e5 )**2 ) / eVToJoule return [ X, Y, Z, r, math.degrees( theta ), math.degrees( phi ), Ekin1, Ekin2, Ekintrans, Ekinrot ] def AtomNumber( symbol ): if symbol == "H": return 1 elif symbol == "He": return 2 elif symbol == "C": return 6 elif symbol == "O": return 8 elif symbol == "Ru": return 44 def PrintXSF( vaspfiles ): # we skip the last one here, it's part of the next trajectory steps = [ len( x.geometries ) - 1 for x in vaspfiles ] print "ANIMSTEPS " + str( sum( steps ) ) print "CRYSTAL" print "PRIMVEC" print " ".join( [ str(x) for x in vaspfiles[0].basis[0][0] ] ) print " ".join( [ str(x) for x in vaspfiles[0].basis[0][1] ] ) print " ".join( [ str(x) for x in vaspfiles[0].basis[0][2] ] ) geom = 0 for k in vaspfiles: for i in range( 0, len( k.geometries ) - 1 ): geom += 1 print "PRIMCOORD " + str(geom) print str(len( k.geometries[i] )) + " 1" for j in range( 0, len( k.geometries[i] ) ): print str( AtomNumber( k.atoms[j] ) ) + " " + str( k.geometries_cartesian[i][j][0] ) + " " + str( k.geometries_cartesian[i][j][1] ) + " " + str( k.geometries_cartesian[i][j][2] ) #print str( AtomNumber( k.atoms[j] ) ) + " " + str( k.velocities[i][j][0] ) + " " + str( k.velocities[i][j][1] ) + " " + str( k.velocities[i][j][2] ) def WriteXDATCAR( vaspfiles, basename ): f = open( basename + ".XDATCAR", "w" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].basis[0][0] ] ) + "\n" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].basis[0][1] ] ) + "\n" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].basis[0][2] ] ) + "\n" ) f.write( " ".join( vaspfiles[0].atomtypes ) + "\n" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].atomcounts ] ) + "\n" ) tstep = 0 for k in vaspfiles: for i in range( 0, k.nsw ): f.write( "Direct configuration " + str( tstep ) + "\n" ) tstep += 1 for j in range( 0, len( k.atoms ) ): f.write( " ".join( [ str(x) for x in k.geometries[i][j] ] ) + "\n" ) f.close() f = open( basename + ".POSCAR", "w" ) f.write( "Dummy for VMD\n" ) f.write( "1.0\n" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].basis[0][0] ] ) + "\n" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].basis[0][1] ] ) + "\n" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].basis[0][2] ] ) + "\n" ) f.write( " ".join( vaspfiles[0].atomtypes ) + "\n" ) f.write( " ".join( [ str(x) for x in vaspfiles[0].atomcounts ] ) + "\n" ) f.close() def WriteEnergies( vaspfiles, basename ): f = open( basename + ".dat", "w" ) tfinal = 0.0 for k in vaspfiles: for i in range( 0, k.nsw ): f.write( "%s %s %s\n" % ( tfinal + i * k.potim, k.energies_e0[i], k.energies_ekin[i] ) ) tfinal += k.nsw * k.potim f.close() def PrintAverage( vaspfiles, basename ): e0list = [] eklist = [] for k in vaspfiles: e0list += k.energies_e0 eklist += k.energies_ekin e0avg = np.mean( e0list ) ekavg = np.mean( eklist ) e0std = np.std( e0list ) ekstd = np.std( eklist ) print basename + ": " + " = " + ( "%.4f" % e0avg ) + " +/- " + ( "%.4f" % e0std ) + " eV" + ", " + " = " + ( "%.4f" % ekavg ) + " +/- " + ( "%.4f" % ekstd ) + " eV" def ArrayDifference( array1, array2 ): for i in range( len(array1) ): for j in range( len(array1[i]) ): if abs( array1[i][j] - array2[i][j] ) > 1e-6: print "Array 1 and 2 are different, element: ", i, j return print "Array 1 and 2 are the same" #print velocities[498] #print velocities[499] #print geometries[0] #print geometries[499] #print geometries[500] files = [] #basename = sys.argv[1] for i in sys.argv[1:]: if os.path.isfile( i ): file = VASPXMLFile( i ) files.append( file ) if len( files ) < 1: raise Exception( "No files given" ) for idx, i in enumerate( files[0].velocities ): print files[0].GetPoscar( idx ) #ekin, temperature = files[0].ComputeTemperature( atommask=[0,1,2,3,4,5,6,7,8,9,10,11] ) #for i in range( 0, len( ekin ) ): # print i, ekin[i], files[0].energies_e0[i] - -.48817816E+02, temperature[i] #for idx, i in enumerate( files[0].velocities ): # print "VELOCITY AT STEP ", idx # print files[0].GetPoscar( idx ) #PrintXSF( files ) #files[0].ComputeTemperature( [ 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 ] ) #geom0 = files[0].geometries_cartesian[0] #geom1 = files[0].geometries_cartesian[1] #geom2 = files[0].geometries_cartesian[2] #velo0 = files[0].velocities[0] #velo1 = files[0].velocities[1] #force0 = files[0].forces[0] #force1 = files[0].forces[1] #basis = [[4.7692608377837704, 0.0000000000000000, 0.0000000000000000], [2.3846304188918865, 4.1303010427949998, 0.0000000000000000], [0.0000000000000000, 0.0000000000000000, 24.7266484043800006]] #xin = [] #xin.append( [0.0000000000000000, 0.0000000000000000, 0.0000000000000000] ) #xin.append( [0.3333333300000000, 0.3333333300000000, 0.0000000000000000] ) #xin.append( [0.6666666700000000, 0.6666666700000000, 0.0000000000000000] ) #xin.append( [0.9995993229707619, 0.0104254528002154, 0.1758762944763794] ) #xin.append( [0.3329089412792469, 0.9991914062553362, 0.0852085174919084] ) #xin.append( [0.3203349046503709, 0.3493018602598180, 0.1702770029746930] ) #xin.append( [0.6724423362747426, 0.3387351792856449, 0.0888385973934579] ) #xin.append( [0.6666530297274932, 0.6772749129192408, 0.1721612313649886] ) #xin.append( [0.0025876768973611, 0.6825319192754531, 0.0838436500023350] ) #xin.append( [0.9989286623516918, 0.0028168736031176, 0.3560446705773232] ) #xin.append( [0.3468865407450996, 0.0101195128856467, 0.2597139934165441] ) #xin.append( [0.3388791515113049, 0.3386605320124267, 0.3487954551439195] ) #xin.append( [0.6648843343552205, 0.3550512069862028, 0.2618363566878614] ) #xin.append( [0.6567615377349245, 0.6815053756400236, 0.3446713332598961] ) #xin.append( [0.0094104163033680, 0.6665193222220812, 0.2658669749670140] ) #xin.append( [0.9205404306381292, 0.0330499410618753, 0.4305214704436808] ) #xin.append( [0.8886153393869977, 0.0232054165272665, 0.4771058985949266] ) #xin = [ DirectToCartesian( x, files[0].basis[0] ) for x in xin ] #vin = [] #vin.append( [0.00000000E+00, 0.00000000E+00, 0.00000000E+00] ) #vin.append( [0.00000000E+00, 0.00000000E+00, 0.00000000E+00] ) #vin.append( [0.00000000E+00, 0.00000000E+00, 0.00000000E+00] ) #vin.append( [-0.57610545E-03, -0.13469326E-02, -0.85473546E-03] ) #vin.append( [-0.23303169E-02, 0.20938568E-02, -0.10496250E-02] ) #vin.append( [-0.12032429E-02, -0.13701917E-02, -0.46645956E-03] ) #vin.append( [0.22551811E-02, 0.64890924E-03, -0.25812406E-03] ) #vin.append( [0.14546777E-02, -0.14059357E-02, 0.23679875E-02] ) #vin.append( [-0.12905740E-02, 0.17032949E-02, 0.36663380E-03] ) #vin.append( [0.11144495E-02, 0.10292183E-04, 0.43102331E-03] ) #vin.append( [-0.12954504E-02, -0.63242078E-03, -0.31933018E-03] ) #vin.append( [-0.11809582E-02, 0.90490784E-03, 0.64873246E-04] ) #vin.append( [0.10474547E-02, -0.24743738E-02, -0.18721890E-02] ) #vin.append( [0.17984530E-02, 0.23159869E-03, 0.21333710E-02] ) #vin.append( [0.43303324E-03, -0.55856289E-03, -0.29806680E-03] ) #vin.append( [-0.45854883E-02, -0.15235388E-02, -0.14129254E-02] ) #vin.append( [-0.22589256E-02, 0.58178645E-04, -0.28334214E-02] ) #print geom0 #print geom1 #print files[0].geometries[0] #print xin #print velo0 #print files[0].velocities[-1] #print vin #sys.exit() ## step 0: compute accelerations (multiplied with timestep) #accel0 = [] #for i in range( 0, len(force0) ): # accel0.append( [ force0[i][j] * 1.60217733e-9 / ( 1.6605402e-27 * files[0].GetMass( files[0].atoms[i] ) ) * files[0].potim * 1.0e-15 * 1.0e-5 for j in range( 0, 3 ) ] ) #accel1 = [] #for i in range( 0, len(force1) ): # accel1.append( [ force1[i][j] * 1.60217733e-9 / ( 1.6605402e-27 * files[0].GetMass( files[0].atoms[i] ) ) * files[0].potim * 1.0e-15 * 1.0e-5 for j in range( 0, 3 ) ] ) # #print "Now doing regular velocity verlet..." #print "------------------------------------" ## velocity verlet: step 1 ## v(t + .5*dt) = v + .5*a(t)*dt #velohalf = [] #for i in range( 0, len(velo0) ): # velohalf.append( [ velo0[i][j] + 0.5 * accel0[i][j] for j in range( 0, 3 ) ] ) # ## velocity verlet: step 2 ## x(t + dt) = x(t) + v(t + .5dt)*dt #newgeom = [] #for i in range( 0, len(geom0) ): # newgeom.append( [ geom0[i][j] + velohalf[i][j] * files[0].potim for j in range( 0, 3 ) ] ) #print #print "Testing geometries..." #ArrayDifference( geom1[3:], newgeom[3:] ) # ## velocity verlet: step 3 ## already done :) # ## velocity verlet: step 4 ## v(t + dt) = v(t + .5dt) + .5*a(t + dt)*dt #newvelo = [] #for i in range( 0, len(velohalf) ): # newvelo.append( [ velohalf[i][j] + .5 * accel1[i][j] for j in range( 0, 3 ) ] ) #print #print "Testing velocities..." #ArrayDifference( velo1[3:], newvelo[3:] ) ##print velo0 ##print velo1 ##print newvelo # #print #print "Now doing leapfrog..." #print "---------------------" ## leapfrog: step 1 ## x(t) = x(t - dt) + v(t - .5*dt)*dt #newgeom = [] #for i in range( 0, len(geom0) ): # newgeom.append( [ geom0[i][j] + velo0[i][j] * files[0].potim for j in range( 0, 3 ) ] ) ##print geom0 ##print geom1 ##print newgeom #print #print "Testing geometries..." #print geom1 #print newgeom #ArrayDifference( geom1[3:], newgeom[3:] ) # ## leapfrog: step 2 ## v(t + .5dt) = v(t - .5dt) + a(t)*dt #newvelo = [] #for i in range( 0, len(velo0) ): # newvelo.append( [ velo0[i][j] + accel1[i][j] for j in range( 0, 3 ) ] ) #print #print "Testing velocities..." #print velo1 #print newvelo #ArrayDifference( velo1[3:], newvelo[3:] ) # #print #print "Now doing VASP leapfrog..." #print "--------------------------" # #print 0.01 + accel0[0][0] #print 4.1 + (0.01+accel0[0][0]) * files[0].potim #newvelo = [] #for i in range( 0, len(velo0) ): # newvelo.append( [ velo0[i][j] + accel1[i][j] for j in range( 0, 3 ) ] ) #print #print "Testing velocities..." ##ArrayDifference( velo1[3:], newvelo[3:] ) #print velo1 #print newvelo ##for i in velo0: ## print i # #newgeom = [] #for i in range( 0, len(geom0 ) ): # newgeom.append( [ geom1[i][j] + newvelo[i][j] * files[0].potim for j in range( 0, 3 ) ] ) #print #print "Testing geometries..." #print newgeom #print geom2 ##print geom1 # #print "TESTING STUFF" #print velo0 #print accel0 #print np.mean( files[0].ComputeTemperature( atommask=[9,11,13] ) ) #PrintAverage( files[2:], basename ) #WriteEnergies( files, basename ) #WriteXDATCAR( files, basename ) #print files[0].velocities[0] #print files[0].velocities[499]