#!/usr/bin/env python ################################################################################ # # # # ################################################################################ # # # F. Nattino # ################################################################################ # # ################################################################################ import sys, os, gzip import numpy as np from ReadTrajectoryVASP import * from math import * def ReadFile( FileName ): if (os.path.splitext( FileName )[1] == ".tgz"): lines = gzip.open( FileName ).readlines() else: lines = open( FileName ).readlines() return lines # Some constants ########################################## amu2au = 1822.88839 Angstrom2au = 1.889725989 fs2au = 41.3413732914 au2K = 315775.04 au2eV = 27.21138505 kb = 8.6173324E-5 MassNamu = 14.001 MassWamu = 183.850 N1AtomPos = 20 N2AtomPos = 21 TimeStep = 1.0 # Input ################################################### KinEnergies = [ 0.9, 1.3, 1.7, 2.287] # eV Ntrajs = 400 MaximumTime = 4500.0 # fs AtomsPerLayer = 4 NumLayers = 4 # Moving layers ########################################################### # Directories where looking for trajectories and name of the files with initial conditions RootDir = "/home/francesco/VASPN2W110/PBE/Dynamics/" for KinEnergy in KinEnergies: FileOutputEnergyTransferName = "EnergyN2VsTime-%5.3feV.dat" %KinEnergy FileOutputEnergyTransfer = open( FileOutputEnergyTransferName, 'w') # Loop over directories looking for reacted trajectories at a certain collision energy Reacted = [] # Directory of set of trajectories Directory = RootDir + "%5.3feV/" %KinEnergy # File containing outcome of trajs FileAnalysisTraj = Directory + 'AnalysisTraj-PeriodicImageInteractionRequired.out' linesAnalysisTraj = open( FileAnalysisTraj, 'r').readlines() for i in range( Ntrajs ): line = linesAnalysisTraj[i].split() if line[4] == "Reaction" : traj = Directory + line[2] # The list Reacted contains the directory name where to look for files Reacted.append( traj ) # Ntraj counts the reacted trajectories TimeValues = np.arange(0.0, MaximumTime) AverageKinEnergyN2Direct = np.zeros( len( TimeValues )) AverageKinEnergyN2Trapping = np.zeros( len( TimeValues )) AverageKinEnergySurfAtomsPerLayer = [] for j in range( NumLayers ): AverageKinEnergySurfAtomsPerLayer.append( np.zeros( len( TimeValues )) ) # print a header in the files string = "#Time #AverageKinEnergyN2Direct\n" FileOutputEnergyTransfer.write( string ) ShortestTimeDirect = MaximumTime ShortestTimeTrapping = MaximumTime NtrajDirect = 0 NtrajTrapping = 0 # Now that the trajectories have been identified loop over them for traj in Reacted: # Read Trajectory Traj = Trajectory( Verbosity=0 ) File = traj + 'POSCAR-1' if os.path.isfile( File ): Traj.ReadPOSCAR( Filename = File ) iFile = 1 while True : File = traj + 'XDATCAR-' + str(iFile) if os.path.isfile( File ): Traj.ReadXDATCAR(Filename = File) iFile = iFile + 1 else: break else: File = traj + 'POSCAR' Traj.ReadPOSCAR( Filename = File ) File = traj + 'XDATCAR' Traj.ReadXDATCAR(Filename = File) Traj.ExtractRealPositionsAndVelocities( TimeStep ) File = traj + 'CONTCAR' Traj.ReadCONTCAR(Filename = File) EInitial = 0. EFinal = 0. Rebounds = 0 GoingUp = False MinZCoM = 6. KinEnergyN2 = np.zeros( Traj.maxnt * TimeStep ) EkinLayer = [] for j in range( NumLayers ) : EkinLayer.append( np.zeros( Traj.maxnt * TimeStep ) ) # Loop over timesteps of the single traj for i in range( Traj.maxnt ): # Define time Time = float( i ) * TimeStep vx = Traj.vCar[i][N1AtomPos][0] vy = Traj.vCar[i][N1AtomPos][1] vz = Traj.vCar[i][N1AtomPos][2] EkinAtom1 = 0.5 * MassNamu * ( vx ** 2. + vy ** 2. + vz ** 2. ) * amu2au * Angstrom2au **2. / fs2au **2. * au2eV vx = Traj.vCar[i][N2AtomPos][0] vy = Traj.vCar[i][N2AtomPos][1] vz = Traj.vCar[i][N2AtomPos][2] EkinAtom2 = 0.5 * MassNamu * ( vx ** 2. + vy ** 2. + vz ** 2. ) * amu2au * Angstrom2au **2. / fs2au **2. * au2eV Ekin = EkinAtom1 + EkinAtom2 KinEnergyN2[i] = Ekin for j in range( NumLayers) : for k in range( AtomsPerLayer ): AtomNumber = j * AtomsPerLayer + k vx = Traj.vCar[i][AtomNumber][0] vy = Traj.vCar[i][AtomNumber][1] vz = Traj.vCar[i][AtomNumber][2] EkinSurfAtom = 0.5 * MassWamu * ( vx ** 2. + vy ** 2. + vz ** 2. ) * amu2au * Angstrom2au **2. / fs2au **2. * au2eV EkinLayer[j][i] = EkinLayer[j][i] + EkinSurfAtom #/ float( AtomsPerLayer ) vCoMz = ( Traj.vCar[i][N1AtomPos][2] + Traj.vCar[i][N2AtomPos][2] ) / 2. # Count the number of rebounds if vCoMz > 0. or GoingUp : GoingUp = True if vCoMz < 0.: Rebounds = Rebounds + 1 GoingUp = False N2Closest = Traj.FindClosestImage( i, N1AtomPos, N2AtomPos ) r = sqrt( ( Traj.rCar[i][N1AtomPos][0] - Traj.rCar[i][N2AtomPos][0] ) ** 2. + \ ( Traj.rCar[i][N1AtomPos][1] - Traj.rCar[i][N2AtomPos][1] ) ** 2. + \ ( Traj.rCar[i][N1AtomPos][2] - Traj.rCar[i][N2AtomPos][2] ) ** 2. ) rsm = sqrt( ( Traj.rCar[i][N1AtomPos][0] - N2Closest[0] ) ** 2. + \ ( Traj.rCar[i][N1AtomPos][1] - N2Closest[1] ) ** 2. + \ ( Traj.rCar[i][N1AtomPos][2] - N2Closest[2] ) ** 2. ) if r > 2. and ( r - rsm ) >= 0.01 : break if Rebounds == 0 : if Time <= ShortestTimeDirect: ShortestTimeDirect = Time AverageKinEnergyN2Direct[0:int( Time )] = AverageKinEnergyN2Direct[0:int( Time )] + KinEnergyN2[0:int( Time )] NtrajDirect = NtrajDirect + 1 if Time > 800. : if Time <= ShortestTimeTrapping: ShortestTimeTrapping = Time AverageKinEnergyN2Trapping[0:int( Time )] = AverageKinEnergyN2Trapping[0:int( Time )] + KinEnergyN2[0:int( Time )] for j in range( NumLayers ): AverageKinEnergySurfAtomsPerLayer[j][0:int( Time )] = AverageKinEnergySurfAtomsPerLayer[j][0:int( Time )] + EkinLayer[j][0:int( Time )] NtrajTrapping = NtrajTrapping + 1 del Traj for i in range( int( ShortestTimeDirect )): string = "% 8.5f % 8.5f\n" %( TimeValues[i], AverageKinEnergyN2Direct[i] / float( NtrajDirect ) ) FileOutputEnergyTransfer.write( string ) string = "#Direct" + str( NtrajDirect ) + "\n\n" FileOutputEnergyTransfer.write( string ) string = "#Time #AverageKinEnergyN2Trapping " for j in range( NumLayers ): string = string + "#KinEnSurfLayer" + str( j ) + " " string = string + "#AverageSurfAtoms\n" FileOutputEnergyTransfer.write( string ) for i in range( int( ShortestTimeTrapping )): string = "% 8.5f % 8.5f" %( TimeValues[i], AverageKinEnergyN2Trapping[i] / float(NtrajTrapping) ) AverageKinEnergySurfAtoms = 0. for j in range( NumLayers ): string = string + "% 8.5f " %( AverageKinEnergySurfAtomsPerLayer[j][i] / float(NtrajTrapping)) AverageKinEnergySurfAtoms = AverageKinEnergySurfAtoms + AverageKinEnergySurfAtomsPerLayer[j][i] / float(NtrajTrapping) #/ float( NumLayers ) string = string + "% 8.5f \n" %AverageKinEnergySurfAtoms FileOutputEnergyTransfer.write( string ) string = "#Trapping" + str( NtrajTrapping ) + "\n" FileOutputEnergyTransfer.write( string ) FileOutputEnergyTransfer.close()