#!/usr/bin/env python ################################################# from os import path # from numpy import array, dot, degrees, average, argmin, sqrt # from numpy import abs as npabs # from numpy.linalg import norm, inv # from molecule import com # from geom import angle_between as ab # ################################################# ##################################################################################################### ### P A R A M E T E R S ############################################################################# ##################################################################################################### R_CRIT_min = 2.0 R_CRIT_max = 3.0 Z_CRIT = 6.5 tstep = 0.4 # [ fs ] #T_MAX = 2000 # [ fs ] T_MAX = 1000 # [ fs ] MPt = 106.420 # atomic masses MC = 12.000 # MH = 1.008 # MD = 2.014 # M = MC + MH + 3*MD # methane mass surface_atoms = 36 # ######################################### # au2Ang = 0.52917721092 # conversion constants amu2au = 1822.88839 # fs2au = 41.341373337 # au2eV = 27.21138505 # deg2rad = 0.0174532925 # J2eV = 6.24181 * 10**18 # h = 6.6260695729 * 10**-34 # # Kb = 8.6173857e-05 # [ eV/K ] Boltzmann amu2kg = 1.6605402e-27 # ang2m = 1e-10 # fs2s = 1e-15 # eV2J = 1.60217733e-19 # # ######################################### ##################################################################################################### ### restart number ################################################################################## runs=-1 while path.exists('XDATCAR_%d' % (runs+1)): runs += 1 if runs == 0: poscar = 'POSCAR' else: poscar = 'POSCAR_0' if path.exists('CONTCAR'): contcar = 'CONTCAR' else: contcar = 'CONTCAR_%d' % runs ##################################################################################################### POSlines = open(poscar).readlines() # read file a = [float(x) for x in POSlines[2].split()] # store cell parameters b = [float(x) for x in POSlines[3].split()] # c = [float(x) for x in POSlines[4].split()] # Cell = array([a, b, c]) # tansformation matrix Cell_= inv(Cell) ### Positions ######################################################################################## C = [ ] # postition of the atoms H = [ ] # for each timestep D1 = [ ] # D2 = [ ] # D3 = [ ] # header = 7 # number of lines of the XDATCAR header block = surface_atoms + 5 + 1 # surf atoms + methane + headline #----------------------------------------------------------------------------------------------------- def XDAT_analyze(inxdat, block): XDATlines = open(inxdat).readlines() steps = int(XDATlines[-(surface_atoms+6)].split()[2]) for i in range(steps): C.append ( [ float(x) for x in XDATlines[( header + i*block + surface_atoms + 1 )].split() ] ) H.append ( [ float(x) for x in XDATlines[( header + i*block + surface_atoms + 2 )].split() ] ) D1.append( [ float(x) for x in XDATlines[( header + i*block + surface_atoms + 3 )].split() ] ) D2.append( [ float(x) for x in XDATlines[( header + i*block + surface_atoms + 4 )].split() ] ) D3.append( [ float(x) for x in XDATlines[( header + i*block + surface_atoms + 5 )].split() ] ) #----------------------------------------------------------------------------------------------------- for i in range(1, (runs+1)): XDAT_analyze('XDATCAR_%d' % i, block) # if 'XDATCAR' is present it means that the calculation if path.exists('XDATCAR'): XDAT_analyze('XDATCAR', block) # is running at the moment. Therefore the steps computed steps = len(C) #for i in range(steps): #print "%10.8f %10.8f %10.8f" % tuple(C[86]) #print "%10.8f %10.8f %10.8f" % tuple(H[86]) #print "%10.8f %10.8f %10.8f" % tuple(D1[86]) #print "%10.8f %10.8f %10.8f" % tuple(D2[86]) #print "%10.8f %10.8f %10.8f" % tuple(D3[86]) #print ### Real Traj ################################################################################# ### Find the real trajectory of the atmos considering the eventual supercell ### supercell border crossing #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def Real_Trajectory(atom, ref, Real): '''atom is an array containing, for each step, [x, y, z] of the atom in direct coordinates ref is the reference atom in the molecule. atom coordinates will be changed so that the real position of the atom will be considered and not necessarly the periodic image in the [0,0,0] cell. The reference atoms is needed to find out which cell to start from''' ref = array(ref) atom = array(atom) r_min = 99999 # absolute distance between the reference and the atom diff = 3*[0.] # componentwise distance for i in range(-1, 2): # find out in which cell is the atom for j in range(-1, 2): # in respect of the reference one dist = ref[0] - (atom[0] + array([i,j,0])) # dist = norm(dist) # this is the distance between the reference atom # and the periodic image in cell [i,j,0] of the tested atom if dist < r_min: # store the periodic image closest to the reference atom r_min = dist # at the beginning of the trajectory a_min = i # b_min = j # # cell = array([a_min , b_min, 0]) # the cell considered is the real [0, 0, 0] + cell[a_min, b_min, 0] Real.append( atom[0] + cell ) # initial step according to the for i in range(1, len(ref)): # assume that the next step will Real.append(atom[i]+cell) # be in the same cell # for j in range(3): # check if that's actually the case diff[j] = Real[i][j] - Real[(i-1)][j] # if diff[j] < -0.5: # if it's not update the position cell[j] += 1. # and the reference cell accordingly Real[i][j] += 1. # elif diff[j] > 0.5: # cell[j] -= 1. # Real[i][j] -= 1. # #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def distance(A,B): "return the distance between two 3D points" return norm(array(A)-array(B)) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C_r = [] # real postition of the atoms H_r = [] # for each timestep D1_r = [] # D2_r = [] # D3_r = [] # Real_Trajectory( C, C, C_r ) Real_Trajectory( H, C, H_r ) Real_Trajectory( D1, C, D1_r ) Real_Trajectory( D2, C, D2_r ) Real_Trajectory( D3, C, D3_r ) #testC = open('TESTC.dat', 'w') #testCr = open('TESTCr.dat', 'w') #DIFF = [ ] #for i in range(len(C)): # # print '%9.6f %9.6f' % (C[i][0], C_r[i][0]) # DIFF.append(norm(abs(array(C[i])-array(C_r[i])))) # testC.write('%9.6f %9.6f %9.6f\n' % tuple(C[i])) # testCr.write('%9.6f %9.6f %9.6f\n' % tuple(C_r[i])) #print max(DIFF), min(DIFF) #testC.close() #testCr.close() ### Direct to Cartesian ############################################################################ C_c = dot(C_r , Cell) # real postition of the atoms H_c = dot(H_r , Cell) # for each timestep D1_c = dot(D1_r , Cell) # in cartesian coordinates D2_c = dot(D2_r , Cell) # D3_c = dot(D3_r , Cell) # ### Velocities #################################################################################### #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def velocities(positions, dt, velocities): steps = len(positions) v = 3 * [0.0] for i in range(1, (steps-1)): for j in range(3): v[j] = ( positions[(i+1)][j] - positions[(i-1)][j] ) / (2*dt) velocities.append( [float(x) for x in v] ) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def kinetic(v, M): v_tot = norm( v ) Ki = 0.5 * M * ( v_tot**2 ) * amu2kg * ang2m**2 / fs2s**2 / eV2J return Ki #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ vC = [] vH = [] vD1 = [] vD2 = [] vD3 = [] #head #pos&vel surface #posCH4 #n vC.append( [float(x) for x in open(poscar).readlines()[ ( 9 + 2*surface_atoms + 5 + 1 ) ].split()] ) # initial velocity is vH.append( [float(x) for x in open(poscar).readlines()[ ( 9 + 2*surface_atoms + 5 + 2 ) ].split()] ) # taken from the POSCAR file vD1.append( [float(x) for x in open(poscar).readlines()[ ( 9 + 2*surface_atoms + 5 + 3 ) ].split()] ) # vD2.append( [float(x) for x in open(poscar).readlines()[ ( 9 + 2*surface_atoms + 5 + 4 ) ].split()] ) # vD3.append( [float(x) for x in open(poscar).readlines()[ ( 9 + 2*surface_atoms + 5 + 5 ) ].split()] ) # velocities( C_c, tstep, vC ) # compute velocities velocities( H_c, tstep, vH ) # velocities( D1_c, tstep, vD1 ) # velocities( D2_c, tstep, vD2 ) # velocities( D3_c, tstep, vD3 ) # vC.append( [float(x) for x in open(contcar).readlines()[ ( 9 + 2*surface_atoms + 5 + 1 ) ].split()] ) # save the last step velocity vH.append( [float(x) for x in open(contcar).readlines()[ ( 9 + 2*surface_atoms + 5 + 2 ) ].split()] ) # form the most recent vD1.append( [float(x) for x in open(contcar).readlines()[ ( 9 + 2*surface_atoms + 5 + 3 ) ].split()] ) # CONTCAR vD2.append( [float(x) for x in open(contcar).readlines()[ ( 9 + 2*surface_atoms + 5 + 4 ) ].split()] ) # vD3.append( [float(x) for x in open(contcar).readlines()[ ( 9 + 2*surface_atoms + 5 + 5 ) ].split()] ) # vC = array(vC) vH = array(vH) vD1 = array(vD1) vD2 = array(vD2) vD3 = array(vD3) ### COM ######################################################################################## COM = [] # position of the CoM vCOM = [] # velocity of the CoM KCOM = [] # kinetic E of the CoM for i in range(steps): CenterOfMass_p = MC * C_c[i] + MH * H_c[i] + MD * ( D1_c[i] + D2_c[i] + D3_c[i] ) CenterOfMass_v = MC * vC[i] + MH * vH[i] + MD * ( vD1[i] + vD2[i] + vD3[i] ) COM.append( CenterOfMass_p/M ) vCOM.append( CenterOfMass_v/M ) k_ics = kinetic(vCOM[-1][0], M) k_ypsilon = kinetic(vCOM[-1][1], M) k_zed = kinetic(vCOM[-1][2], M) KCOM.append( [ k_ics , k_ypsilon , k_zed ] ) ### Bonds Lenght ################################################################################## #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def rVSt(A1, A2, r): printWARNING = False steps = len(A1) for k in range(steps): r.append( norm(A1[k]-A2[k]) ) if r[-1] > 3.5: printWARNING = True if printWARNING: open('BONDveryLONG', 'w').write("\n") #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ rCH = [] rCD1 = [] rCD2 = [] rCD3 = [] rVSt(C_c, H_c, rCH) rVSt(C_c, D1_c, rCD1) rVSt(C_c, D2_c, rCD2) rVSt(C_c, D3_c, rCD3) ### Evaluate Outcome ################################################################################# #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def reacted(radius, positions, label, R_CRIT_min, R_CRIT_max): '''given a vector containing the radius (bond distance) at each timestep, it returns True if the bond has been broken. the bond is considered broken if it is larger than R_CRIT_max or if it stays larger than R_CRIT_min for TRESHOLD timesteps''' BACK = False # True if the molecule back reacts: r > R_CRIT_min fot t < TRESHOLD REACT = False TRESHOLD = int(float(100)/tstep) # if reacted propagate for 100 fs more steps = len(radius) d_steps = 0 react_t = 9999. REF = array(C_c) for k in range(steps): # compute dist from the atom to all the first images of C alldist = [ ] for m in range(-1,2): for n in range(-1,2): alldist.append(distance(positions[k], REF[k] + m*Cell[0] + n*Cell[1])) if d_steps > 0 and radius[k] < R_CRIT_min: BACK = True d_steps = 0 elif radius[k] > R_CRIT_max: REACT = True elif alldist[4] > min(alldist): REACT = True elif radius[k] > R_CRIT_min and d_steps < TRESHOLD: if d_steps == 0: react_t = k * tstep d_steps += 1 elif d_steps > TRESHOLD: REACT = True return label, REACT, BACK, react_t #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ reaction = [] reaction.append( reacted( rCH , H_c, 'H' , R_CRIT_min, R_CRIT_max ) ) reaction.append( reacted( rCD1, D1_c, 'D1' , R_CRIT_min, R_CRIT_max ) ) reaction.append( reacted( rCD2, D2_c, 'D2' , R_CRIT_min, R_CRIT_max ) ) reaction.append( reacted( rCD3, D3_c, 'D3' , R_CRIT_min, R_CRIT_max ) ) SCATTERING = False BACK_REACTION = False REACTION = False TRAPPING = False if steps*tstep > T_MAX: TRAPPING = True else: for i in range(steps): if COM[i][2] > Z_CRIT and vCOM[i][2] > 0 and SCATTERING == False: SCATTERING = True scat_t = i*tstep N_bond_broken = 0 for i in range(4): if reaction[i][1] == True: N_bond_broken += 1 bond_broken = reaction[i][0] react_t = reaction[i][3] REACTION = True if reaction[i][2] == True: BACK_REACTION = True ### Angles ########################################################################################### def Angles( DISS, H1, H2, H3, C ): "reurns beta and theta in respect of the dissociating atom DISS" DISS, H1, H2, H3, C = array(DISS), array(H1), array(H2), array(H3), array(C) THETA, BETA, GAMMA = [ ], [ ], [ ] for j in range( steps ): CH_bond = DISS[j] - C[j] # dissociating bond baricentro = array(com( [ H1[j], H2[j], H3[j] ], [1., 1., 1.])) # umbrella geometric center == com if the masses are all the same umbrella = C[j] - baricentro THETA.append( degrees(ab(CH_bond , [0,0,1])) ) BETA.append( degrees(ab(umbrella , [0,0,1])) ) GAMMA.append( degrees(ab(CH_bond , umbrella)) ) return THETA, BETA, GAMMA # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if REACTION and bond_broken == 'H' : theta, beta, gamma = Angles(H_c, D1_c, D2_c, D3_c, C_c) if REACTION and bond_broken == 'D1': theta, beta, gamma = Angles(D1_c, H_c, D2_c, D3_c, C_c) if REACTION and bond_broken == 'D2': theta, beta, gamma = Angles(D2_c, H_c, D1_c, D3_c, C_c) if REACTION and bond_broken == 'D3': theta, beta, gamma = Angles(D3_c, H_c, D2_c, D1_c, C_c) if SCATTERING or TRAPPING : theta, beta, gamma = Angles(H_c, D1_c, D2_c, D3_c, C_c) if REACTION or SCATTERING or TRAPPING: output_angles = open( 'angles.dat', 'w' ) output_angles.write('# theta beta gamma\n') for i in range( steps ): output_angles.write( ' %7.4f %7.4f %7.4f\n' % ( theta[i], beta[i], gamma[i] ) ) output_angles.close() ### Outputs ########################################################################################## output_r = open('analysis_r.dat', 'w') output_c = open('analysis_com.dat', 'w') output_r.write('#t r_CH r_CD1 r_CD2 r_CD3#\n') output_c.write('#t Z X Y vZ Kz\n') for i in range(steps): output_r.write('%4.1f %12.8f %12.8f %12.8f %12.8f\n' % ((i*tstep), rCH[i], rCD1[i], rCD2[i], rCD3[i])) output_c.write('%4.1f %12.8f %12.8f %12.8f %12.8f %12.8f\n' % ((i*tstep), COM[i][2], COM[i][0], COM[i][1], vCOM[i][2], KCOM[i][2])) output_r.close() output_c.close() ### Reaction outcome ################################################################################ #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def CountBounces(vel, FinalStep=steps): ''' Computes the number of bounces given a vector with the velocity of the CoM. One bounce is counted when the vel changes sign twice''' first_bounce = 0 sign_changes = 0 for i in range(1, FinalStep): if vel[i]*vel[i-1] < 0.0 : if sign_changes == 0: first_bounce = i sign_changes += 1 return (sign_changes/2), first_bounce #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def XyDist(pos,FinalStep=steps): ''' Computes the distance tavelled in the XY plane given the COM positions ''' dist_tot, dist_ics, dist_why = 0., 0., 0. for i in range(1,FinalStep): dist_tot += norm( array(pos[i][0:2])-array(pos[i-1][0:2]) ) dist_ics += npabs(pos[i][0]-pos[i-1][0]) dist_why += npabs(pos[i][1]-pos[i-1][1]) if FinalStep < len(pos): displace = norm( array(pos[FinalStep][0:2]) - array(pos[0][0:2]) ) else: displace = norm( array(pos[-1][0:2]) - array(pos[0][0:2]) ) return dist_tot, dist_ics, dist_why, displace #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ def Dist2Surf(pos, FinalStep): ''' Compute the distance from methane to the closest top, fcc, hcp and bridge site at the time of the reaction ''' PosCart = pos[FinalStep] surface1layer = array( [ [ 0.0000000000000, 0.0000000000000, 0. ], [ 0.0000000000000, 0.3333333333333, 0. ], [ 0.0000000000000, 0.6666666666667, 0. ], [ 0.0000000000000, 1.0000000000000, 0. ], [ 0.3333333333333, 0.0000000000000, 0. ], [ 0.3333333333333, 0.3333333333333, 0. ], [ 0.3333333333333, 0.6666666666667, 0. ], [ 0.3333333333333, 1.0000000000000, 0. ], [ 0.6666666666667, 0.0000000000000, 0. ], [ 0.6666666666667, 0.3333333333333, 0. ], [ 0.6666666666667, 0.6666666666667, 0. ], [ 0.6666666666667, 1.0000000000000, 0. ], [ 1.0000000000000, 0.0000000000000, 0. ], [ 1.0000000000000, 0.3333333333333, 0. ], [ 1.0000000000000, 0.6666666666667, 0. ], [ 1.0000000000000, 1.0000000000000, 0. ] ] ) surface2layer = array( [ [ 0.1111111111111, 0.1111111111111, 0. ], [ 0.1111111111111, 0.4444444444444, 0. ], [ 0.1111111111111, 0.7777777777778, 0. ], [ 0.1111111111111, 1.1111111111111, 0. ], [ 0.4444444444444, 0.1111111111111, 0. ], [ 0.4444444444444, 0.4444444444444, 0. ], [ 0.4444444444444, 0.7777777777778, 0. ], [ 0.4444444444444, 1.1111111111111, 0. ], [ 0.7777777777778, 0.1111111111111, 0. ], [ 0.7777777777778, 0.4444444444444, 0. ], [ 0.7777777777778, 0.7777777777778, 0. ], [ 0.7777777777778, 1.1111111111111, 0. ], [ 1.1111111111111, 0.1111111111111, 0. ], [ 1.1111111111111, 0.4444444444444, 0. ], [ 1.1111111111111, 0.7777777777778, 0. ], [ 1.1111111111111, 1.1111111111111, 0. ] ] ) surface3layer = array( [ [ -0.1111111111111, 0.2222222222222, 0. ], [ -0.1111111111111, 0.5555555555556, 0. ], [ -0.1111111111111, 0.8888888888889, 0. ], [ -0.1111111111111, -0.1111111111111, 0. ], [ 0.2222222222222, 0.2222222222222, 0. ], [ 0.2222222222222, 0.5555555555556, 0. ], [ 0.2222222222222, 0.8888888888889, 0. ], [ 0.2222222222222, -0.1111111111111, 0. ], [ 0.5555555555556, 0.2222222222222, 0. ], [ 0.5555555555556, 0.5555555555556, 0. ], [ 0.5555555555556, 0.8888888888889, 0. ], [ 0.5555555555556, -0.1111111111111, 0. ], [ 0.8888888888889, 0.2222222222222, 0. ], [ 0.8888888888889, 0.5555555555556, 0. ], [ 0.8888888888889, 0.8888888888889, 0. ], [ 0.8888888888889, -0.1111111111111, 0. ] ] ) surface4layer = array( [ [ 0.1666666666667, 0.0000000000000, 0. ], [ 0.1666666666667, 0.3333333333333, 0. ], [ 0.1666666666667, 0.6666666666667, 0. ], [ 0.1666666666667, 1.0000000000000, 0. ], [ 0.5000000000000, 0.0000000000000, 0. ], [ 0.5000000000000, 0.3333333333333, 0. ], [ 0.5000000000000, 0.6666666666667, 0. ], [ 0.5000000000000, 1.0000000000000, 0. ], [ 0.8333333333333, 0.0000000000000, 0. ], [ 0.8333333333333, 0.3333333333333, 0. ], [ 0.8333333333333, 0.6666666666667, 0. ], [ 0.8333333333333, 1.0000000000000, 0. ], [ 0.0000000000000, 0.1666666666667, 0. ], [ 0.0000000000000, 0.5000000000000, 0. ], [ 0.0000000000000, 0.8333333333333, 0. ], [ 0.3333333333333, 0.1666666666667, 0. ], [ 0.3333333333333, 0.5000000000000, 0. ], [ 0.3333333333333, 0.8333333333333, 0. ], [ 0.6666666666667, 0.1666666666667, 0. ], [ 0.6666666666667, 0.5000000000000, 0. ], [ 0.6666666666667, 0.8333333333333, 0. ], [ 1.0000000000000, 0.1666666666667, 0. ], [ 1.0000000000000, 0.5000000000000, 0. ], [ 1.0000000000000, 0.8333333333333, 0. ], [ 0.1666666666667, 0.1666666666667, 0. ], [ 0.1666666666667, 0.5000000000000, 0. ], [ 0.1666666666667, 0.8333333333333, 0. ], [ 0.5000000000000, 0.1666666666667, 0. ], [ 0.5000000000000, 0.5000000000000, 0. ], [ 0.5000000000000, 0.8333333333333, 0. ], [ 0.8333333333333, 0.1666666666667, 0. ], [ 0.8333333333333, 0.5000000000000, 0. ], [ 0.8333333333333, 0.8333333333333, 0. ] ] ) surface1layer = [ x[0:2] for x in dot(surface1layer, Cell) ] surface2layer = [ x[0:2] for x in dot(surface2layer, Cell) ] surface3layer = [ x[0:2] for x in dot(surface3layer, Cell) ] surface4layer = [ x[0:2] for x in dot(surface4layer, Cell) ] # report COM in the 1st unit cell PosDir = dot(PosCart, Cell_) PosDir %= 1 PosDir %= 1 PosCart = dot(PosDir,Cell) dist_top = 1.e+10 for i in range(len(surface1layer)): newdist = distance(PosCart[0:2],surface1layer[i]) if newdist < dist_top : dist_top = newdist dist_fcc = 1.e+10 for i in range(len(surface2layer)): newdist = distance(PosCart[0:2],surface2layer[i]) if newdist < dist_fcc : dist_fcc = newdist dist_hcp = 1.e+10 for i in range(len(surface3layer)): newdist = distance(PosCart[0:2],surface3layer[i]) if newdist < dist_hcp : dist_hcp = newdist dist_bridge = 1.e+10 for i in range(len(surface4layer)): newdist = distance(PosCart[0:2],surface4layer[i]) if newdist < dist_bridge : dist_bridge = newdist if dist_top < dist_fcc and dist_top < dist_hcp and dist_top < dist_bridge: site = 'top' elif dist_fcc < dist_top and dist_fcc < dist_hcp and dist_fcc < dist_bridge: site = 'fcc' elif dist_hcp < dist_top and dist_hcp < dist_fcc and dist_hcp < dist_bridge: site = 'hcp' else: site = 'bridge' return dist_top, dist_fcc, dist_hcp, dist_bridge, site #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if REACTION: # step when r[i] == r[TransitionState] rTS_CH = 1.61 # TS value for breaking CH TSstep =-1 difference = [ 9999., 9999., 9999., 9999. ] while all( [ x > 4.e-2 for x in difference ] ): TSstep += 1 difference[0] = npabs(rCH[TSstep] - rTS_CH) difference[1] = npabs(rCD1[TSstep] - rTS_CH) difference[2] = npabs(rCD2[TSstep] - rTS_CH) difference[3] = npabs(rCD3[TSstep] - rTS_CH) if bond_broken == 'CH' and argmin(difference) != 0: open('ANALYSIS_PROBLEM', 'w').write('Cleavage mismatch\n') if bond_broken == 'CD1' and argmin(difference) != 1: open('ANALYSIS_PROBLEM', 'w').write('Cleavage mismatch\n') if bond_broken == 'CD2' and argmin(difference) != 2: open('ANALYSIS_PROBLEM', 'w').write('Cleavage mismatch\n') if bond_broken == 'CD3' and argmin(difference) != 3: open('ANALYSIS_PROBLEM', 'w').write('Cleavage mismatch\n') xy_dist, x_dist, y_dist, disp = XyDist(COM, TSstep) bounces, firstbounce = CountBounces([ x[2] for x in vCOM ], TSstep ) dist2topTS, dist2fccTS, dist2hcpTS, dist2bridgeTS, siteTS = Dist2Surf(COM, TSstep) directCOM = dot(COM[TSstep], Cell_) directCOM %= 1 directCOM %= 1 reactCOM = dot(directCOM, Cell) steering = distance( COM[0][0:2], COM[TSstep][0:2] ) else: xy_dist, x_dist, y_dist, disp = XyDist(COM) bounces, firstbounce = CountBounces([ x[2] for x in vCOM ]) if TRAPPING: Vx_y = [ x[0:2] for x in vCOM[firstbounce:] ] Vxy = [ ] Vx = [ ] for i in range(len(Vx_y)): Vxy.append( sqrt(Vx_y[i][0]**2 + Vx_y[i][1]**2) ) for i in range(len(Vx_y)): Vx.append( sqrt(Vx_y[i][0]**2) ) avgVx = average(Vx) avgVxy = average(Vxy) if TRAPPING or SCATTERING: for i in range(1, steps-1): if COM[i-1][2] > COM[i][2] < COM[i+1][2]: closest_step = i break DistOnClosestImpact = distance(COM[0][0:2],COM[closest_step][0:2]) steering = distance( COM[0][0:2], COM[closest_step][0:2] ) if SCATTERING: avgVx = 9999. avgVxy = 9999. dist2topI, dist2fccI, dist2hcpI, dist2bridgeI, siteI = Dist2Surf(COM, 0) initCOM = dot(COM[0], Cell_) initCOM %= 1 initCOM %= 1 initCOM = dot(initCOM, Cell) ################################################################################################### if REACTION or TRAPPING or SCATTERING: output_b = open('PostAnalysis.dat', 'w') output_b.write('-------------------------------------------------\n') output_b.write('Trajectory Post Analysis ------------------------\n') output_b.write('-------------------------------------------------\n') output_b.write('\n%d bounces' % bounces) if TRAPPING or SCATTERING: output_b.write('\ttrapping : %9.4f | %9.4f | xy dist to closest approach %10.4f' % ( avgVxy, avgVx, DistOnClosestImpact ) ) output_b.write('\n\nDistance travelled ----------------------------\n' ) output_b.write('%10.4f Ang in xy\n' % xy_dist) output_b.write('%10.4f Ang in x\n' % x_dist) output_b.write('%10.4f Ang in y\n' % y_dist) output_b.write('Total shift -------------------------------------\n' ) output_b.write('%10.4f Ang in xy\n' % disp ) output_b.write('\ndistance from closest top site') output_b.write('\n %9.4f Initial' % dist2topI ) if REACTION: output_b.write('\t %9.4f TS' % dist2topTS ) output_b.write('\ndistance from closest fcc site') output_b.write('\n %9.4f Initial' % dist2fccI ) if REACTION: output_b.write('\t %9.4f TS' % dist2fccTS ) output_b.write('\ndistance from closest hcp site') output_b.write('\n %9.4f Initial' % dist2hcpI ) if REACTION: output_b.write('\t %9.4f TS' % dist2hcpTS ) output_b.write('\ndistance from closest bridge site') output_b.write('\n %9.4f Initial' % dist2bridgeI ) if REACTION: output_b.write('\t %9.4f TS' % dist2bridgeTS ) output_b.write('\n Closest initial site: %s' % siteI ) if REACTION: output_b.write('\n Closest TS site: %s' % siteTS ) output_b.write('\n\nTHETA ---------------------------------------' ) output_b.write('\n %9.4f Initial' % theta[0] ) if REACTION: output_b.write('\t %9.4f TS' % theta[TSstep] ) output_b.write('\nBETA ------------------------------------------' ) output_b.write('\n %9.4f Initial' % beta[0] ) if REACTION: output_b.write('\t %9.4f TS' % beta[TSstep] ) output_b.write('\nGAMMA -----------------------------------------' ) output_b.write('\n %9.4f Initial' % gamma[0] ) if REACTION: output_b.write('\t %9.4f TS' % gamma[TSstep] ) output_b.write('\n\nCOM initial %12.8f %12.8f\n' % ( initCOM[0], initCOM[1])) if REACTION: output_b.write('COM reaction %12.8f %12.8f\n'% (reactCOM[0], reactCOM[1])) output_b.write('steering in xy %10.4f Ang \n' % steering ) output_b.close() if [ REACTION, TRAPPING, SCATTERING ].count(True) > 1: open('STOPCAR', 'w').write('LSTOP = .TRUE.') open('outcome', 'w').write('WARINIG: MULTIPLE OUTCOMES DETECTED\n') open('WARNING', 'a').write('WARINIG: MULTIPLE OUTCOMES DETECTED\n') open('WARNING', 'a').write('REACT, BACK, TRAP, SCAT == %s, %s, %s, %s\n' % tuple([ REACTION, BACK_REACTION, TRAPPING, SCATTERING ] )) elif N_bond_broken > 1: open('STOPCAR', 'w').write('LSTOP = .TRUE.') open('outcome', 'w').write('WARINIG: %d critic bond lengths\n' % N_bond_broken ) open('WARNING', 'a').write('WARINIG: %d critic bond lengths\n' % N_bond_broken ) elif REACTION and N_bond_broken == 1: open('STOPCAR', 'w').write('LSTOP = .TRUE.') if BACK_REACTION: open('outcome', 'w').write('REACTION - %s\nafter %4.1f fs\n-BACK-\n' % (bond_broken, react_t)) else: open('outcome', 'w').write('REACTION - %s\nafter %4.1f fs\n' % (bond_broken, react_t)) elif TRAPPING: open('STOPCAR', 'w').write('LSTOP = .TRUE.') open('outcome', 'w').write('TRAPPING\ntotal propagation time: %4.1f\n' % (steps*tstep) ) elif SCATTERING: open('STOPCAR', 'w').write('LSTOP = .TRUE.') if BACK_REACTION: open('outcome', 'w').write('SCATTERING - BACK REACTION\nafter %4.1f fs\n' % (scat_t)) else : open('outcome', 'w').write('SCATTERING\nafter %4.1f fs\n' % (scat_t)) else: open('outcome', 'w').write('UNCLEAR \n')