#!/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 MassPtamu = 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 ReboundsToConsider = 5 # from 0 to ReboundsToConsider - 1 ########################################################### # Directories where looking for trajectories and name of the files with initial conditions RootDir = "~/VASPN2W110/PBE/Dynamics/" for KinEnergy in KinEnergies: FileOutputEnergyTransferName = "EnergyTransferDistribution-%5.3feV-Sigma60meV.dat" %KinEnergy FileOutputEnergyTransfer = open( FileOutputEnergyTransferName, 'w') FileOutputMinZCoMName = "MinimumZCoMDistribution-%5.3feV.dat" %KinEnergy FileOutputMinZCoM = open( FileOutputMinZCoMName, 'w') # Loop over directories looking for scattered trajectories at a certain collision energy Scattered = [] # 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] == "Scattered" or line[4] == "Backreaction" : traj = Directory + line[2] # The list Scattered contains the directory name where to look for files Scattered.append( traj ) # Ntraj counts the scattered trajectories per number or rebounds Ntraj = [] EnergyTransferValues = np.linspace(-0.5, 2.5, num=600, endpoint=False) DeltaEnergyTransfer = EnergyTransferValues[1] - EnergyTransferValues[0] SigmaEnergyTransfer = 0.060 MinZCoMValues = np.linspace(0.0, 6.0, num=600, endpoint=False) DeltaMinZCoM = MinZCoMValues[1] - MinZCoMValues[0] SigmaMinZCoM = 0.020 EnergyTransferDistributionPerRebound = [] MinZCoMDistributionPerRebound = [] Sum = [] SumSq = [] for i in range( ReboundsToConsider ): Ntraj.append( 0 ) EnergyTransferDistributionPerRebound.append( np.zeros( len( EnergyTransferValues ))) MinZCoMDistributionPerRebound.append( np.zeros( len( MinZCoMValues ))) Sum.append( 0. ) SumSq.append( 0. ) # print a header in the files string = "#EnergyTransfer " for i in range( ReboundsToConsider ): string = string + "#%1d " %i string = string + "\n" FileOutputEnergyTransfer.write( string ) string = "#MinimumZCoM " for i in range( ReboundsToConsider ): string = string + "#%1d " %i string = string + "\n" FileOutputMinZCoM.write( string ) # Now that the trajectories have been identified loop over them for traj in Scattered: # 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. # Loop over timesteps of the single traj and find initial and final total energies as maximum initial and final kinetic energies 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 ZCoM = ( Traj.rCar[i][N1AtomPos][2] + Traj.rCar[i][N2AtomPos][2] ) / 2. vCoMz = ( Traj.vCar[i][N1AtomPos][2] + Traj.vCar[i][N2AtomPos][2] ) / 2. # Total initial energy: max(Ekin) for molecule approaching the surface and Z > 5 Ang if vCoMz < 0. and ZCoM > 5.: if Ekin > EInitial : EInitial = Ekin # Count rebounds if vCoMz > 0. or GoingUp : GoingUp = True if vCoMz < 0.: Rebounds = Rebounds + 1 GoingUp = False # Total final energy: max(Ekin) for molecule leaving the surface and Z > 5 Ang if vCoMz > 0. and ZCoM > 5. and ZCoM < 7.: if Ekin > EFinal : EFinal = Ekin if ZCoM < MinZCoM: MinZCoM = ZCoM if ZCoM > 7.: break Etransfer = EInitial - EFinal # consider scattered trajectories with 0, 1, 2 .. and >= len(Ntraj) - 1 rebounds. if Rebounds >= len( Ntraj ): Rebounds = len( Ntraj ) - 1 Ntraj[ Rebounds ] = Ntraj[ Rebounds ] + 1 Sum[ Rebounds ] = Sum[ Rebounds ] + Etransfer SumSq[ Rebounds ] = SumSq[ Rebounds ] + Etransfer**2. # Calculate distributions for i in range( len( EnergyTransferValues )): EnergyTransfer = EnergyTransferValues[ i ] EnergyTransferDistributionPerRebound[ Rebounds ][i] = EnergyTransferDistributionPerRebound[ Rebounds ][i] + exp( -( EnergyTransfer - Etransfer )**2. / ( 2. * SigmaEnergyTransfer **2. )) for i in range( len( MinZCoMValues )): Value = MinZCoMValues[i] MinZCoMDistributionPerRebound[ Rebounds ][i] = MinZCoMDistributionPerRebound[ Rebounds ][i] + exp( -( Value - MinZCoM ) **2. / ( 2. * SigmaMinZCoM **2. )) del Traj # Write distributions for i in range(len( EnergyTransferValues )): string = "% 8.5f " %EnergyTransferValues[i] for j in range( ReboundsToConsider ): string = string + " % 8.5f " % (EnergyTransferDistributionPerRebound[ j ][ i ] / sum( EnergyTransferDistributionPerRebound[ j ] ) / DeltaEnergyTransfer ) string = string + "\n" FileOutputEnergyTransfer.write( string ) for i in range(len( MinZCoMValues )): string = "% 8.5f " %MinZCoMValues[i] for j in range( ReboundsToConsider ): string = string + " % 8.5f " % (MinZCoMDistributionPerRebound[ j ][ i ] / sum( MinZCoMDistributionPerRebound[ j ] ) / DeltaMinZCoM ) string = string + "\n" FileOutputMinZCoM.write( string ) string = "#NTrajs " for j in range( ReboundsToConsider ): string = string + "#%-3d " %Ntraj[j] string = string + "\n" FileOutputEnergyTransfer.write( string ) # Calculate average values per each rebound (and add condition: more than one rebound) SumMoreThanOneRebound = 0. NtrajMoreThanOneRebound = 0 string = "#Average " for j in range( ReboundsToConsider ): string = string + "#%-8.5f " %( Sum[j] / float( Ntraj[j] )) if j > 0: SumMoreThanOneRebound = SumMoreThanOneRebound + Sum[j] NtrajMoreThanOneRebound = NtrajMoreThanOneRebound + Ntraj[j] string = string + "\n" FileOutputEnergyTransfer.write( string ) SumSqMoreThanOneRebound = 0. string = "#StError " for j in range( ReboundsToConsider ): string = string + "#%-8.5f " %( sqrt( ( SumSq[ j] / float( Ntraj[j] ) - (Sum[j] / float( Ntraj[j] ) )**2.) / float( Ntraj[j] ) ) ) if j > 0: SumSqMoreThanOneRebound = SumSqMoreThanOneRebound + SumSq[ j] string = string + "\n" FileOutputEnergyTransfer.write( string ) Average = SumMoreThanOneRebound/ float( NtrajMoreThanOneRebound ) StErrorMoreThanOneRebound = sqrt( (SumSqMoreThanOneRebound / float( NtrajMoreThanOneRebound ) - (SumMoreThanOneRebound / float( NtrajMoreThanOneRebound ) )**2. ) / float( NtrajMoreThanOneRebound )) string = "#AverageAndErrorMoreThanOneRebound %-8.5f %-8.5f " %(Average, StErrorMoreThanOneRebound) FileOutputEnergyTransfer.write( string ) FileOutputEnergyTransfer.close() FileOutputMinZCoM.close()