#! /usr/bin/env python # F. Nattino # July 2013, Leiden import sys import string import math import Constants class Slab: def __init__( self, Packing, Surface, A, Zinterlayer, N, Vacuum , Verbose=False ): if Verbose: # Print Input print " \n Lattice constant: %15.13f " %A print " Generate slab geometry for a ", N[0], " x ", N[1], " unit cell" print " Interlayer distances: ", Zinterlayer print " Vacuum space: %4.3f " %Vacuum print " Packing structure of the metal: ", Packing print " Surface : ", Surface self.Basis = [] self.v1 = [] self.v2 = [] self.Packing = Packing self.Surface = Surface self.A = A self.Zinterlayer = Zinterlayer self.N = N self.Vacuum = Vacuum self.Verbose = Verbose self.Cartesian = [] self.Direct = [] def GenerateBasisVectors( self ): # Define primitive surface unit cell if self.Packing == 'fcc' : if self.Surface == '111' : self.v1 = [ self.A / math.sqrt( 2. ) , 0. ] self.v2 = [ - self.A / ( 2. * math.sqrt( 2. ) ), self.A / 2. * math.sqrt( 3./2. ) ] else : print " No data for ", self.Packing, self.Surface, " surface. Abort. " quit() elif self.Packing == 'bcc' : if self.Surface == '110' : self.v1 = [ self.A * math.sqrt( 2. ) / 2. , - self.A / 2. ] self.v2 = [ self.A * math.sqrt( 2. ) / 2. , self.A / 2. ] elif self.Surface == '100' : self.v1 = [ self.A , 0. ] self.v2 = [ 0. , self.A ] else : print " No data for ", self.Packing, self.Surface, " surface. Abort. " quit() else : print " No data for ", self.Packing, " packing. Abort. " quit() # Define supercell self.Basis.append( [ self.N[0] * self.v1[0] , self.N[0] * self.v1[1] , 0. ] ) self.Basis.append( [ self.N[1] * self.v2[0] , self.N[1] * self.v2[1] , 0. ] ) self.Basis.append( [ 0. , 0. , sum( self.Zinterlayer ) + self.Vacuum ] ) if self.Verbose: print " \n Primitive Surface Unit cell: \n " print " % 15.13f % 15.13f " %tuple(self.v1) print " % 15.13f % 15.13f \n " %tuple(self.v2) print " \n Supercell: \n " print " % 15.13f % 15.13f % 15.13f " %tuple(self.Basis[0]) print " % 15.13f % 15.13f % 15.13f " %tuple(self.Basis[1]) print " % 15.13f % 15.13f % 15.13f \n " %tuple(self.Basis[2]) return self.Basis def GenerateCoordinates( self, CoordinateType='Cartesian' ): # List of the translation for atoms in a layer TranslationsList = [] for i in range( self.N[0] ): for j in range( self.N[1] ): TranslationsList.append( ( float(i), float(j) ) ) # List of the relative translations between layers LayerShiftList = [] if self.Packing == 'fcc' : if self.Surface == '111' : LayerShiftList.append( ( 0. , 0. )) # Shift of the first layer with respect to first LayerShiftList.append( ( self.v1[0] / 2. , self.v2[1] / 3. )) # Shift of the second layer with respect to first LayerShiftList.append( ( self.v1[0] , self.v2[1] * 2. / 3. )) # Shift of the third layer with respect to first ModuleElement = 3 elif self.Packing == 'bcc' : if self.Surface == '110' : LayerShiftList.append( ( 0. , 0. )) LayerShiftList.append( ( self.v1[0] , 0. )) ModuleElement = 2 if self.Surface == '100' : LayerShiftList.append( ( 0. , 0. )) LayerShiftList.append( ( self.v1[0] / 2. , self.v2[1] / 2. )) ModuleElement = 2 Coordinates = [ 0., 0., 0. ] if self.Verbose: print " \n Equilibrium Cartesian Coordinates ( Ang ): \n " for i in range( len( self.Zinterlayer ) + 1 ): j = i % ModuleElement LayerShift = LayerShiftList[j] if i != 0 : Coordinates[2] = math.fabs( Coordinates[2] + self.Basis[2][2] - self.Zinterlayer[ i -1 ] ) % self.Basis[2][2] for Translation in TranslationsList : Position = [ 0. , 0. , 0. ] Position[0] = Coordinates[0] + LayerShift[0] + ( Translation[0]* self.v1[0] + Translation[1] * self.v2[0] ) Position[1] = Coordinates[1] + LayerShift[1] + ( Translation[0]* self.v1[1] + Translation[1] * self.v2[1] ) Position[2] = Coordinates[2] self.Cartesian.append( Position ) if self.Verbose: print " % 15.13f % 15.13f % 15.13f " %tuple(Position) if self.Verbose: print " \n Equilibrium Direct Coordinates: \n " for Atom in self.Cartesian : Position = self.CartesianToDirect( Atom ) self.Direct.append( Position ) if self.Verbose: print " % 15.13f % 15.13f % 15.13f " %tuple( Position ) if self.Verbose: print " \n " if CoordinateType == 'Cartesian': return self.Cartesian elif CoordinateType == 'Direct': return self.Direct def CartesianToDirect( self, Cart ): # Transform cartesian coordinates to direct InvertBasis = self.Invert3x3Matrix( self.Basis ) Dir = [ 0., 0., 0. ] Dir[0] = Cart[0] * InvertBasis[0][0] + Cart[1] * InvertBasis[1][0] + Cart[2] * InvertBasis[2][0] Dir[1] = Cart[0] * InvertBasis[0][1] + Cart[1] * InvertBasis[1][1] + Cart[2] * InvertBasis[2][1] Dir[2] = Cart[0] * InvertBasis[0][2] + Cart[1] * InvertBasis[1][2] + Cart[2] * InvertBasis[2][2] return Dir def InitializeTemperature( self, Temperature, Frequency ): from random import gauss # Compute values for Surface temperature Displacement StDevVel = math.sqrt( BoltzmannConstant * Temperature / ( MassAmu * AmuToKg ) ) StDevPos = [] for i in range( len( Frequency ) ): StDevPosLayer = [] for j in range(3): Omega = Frequency[i][j] * eVToJoule * 1.e-3 / DiracConstant StDevPosLayerQ = math.sqrt( BoltzmannConstant * Temperature / ( MassAmu * AmuToKg * Omega **2.) ) StDevPosLayer.append( StDevPosLayerQ ) StDevPos.append( StDevPosLayer ) if self.Verbose: print " \n Standard deviation for generation of velocities: % 15.13f \n" %StDevVel for i in range( len( Frequency ) ): print " Standard deviation ( m ) for generation of displacements (X,Y,Z) for layer ", i, " : " print " %.6e %.6e %.6e " %tuple( StDevPos[i]) # Print Cartesian Coordinates with temperature print " \n Cartesian Coordinates with temperature ( Ang ): \n " Coordinates = [ 0., 0., 0. ] Potential = 0. for i in range( len( self.Cartesian ) ): LayerNumber = i // ( self.N[0] * self.N[1] ) if LayerNumber <= ( len( Frequency ) - 1 ): for j in range(3): dq = gauss( 0., StDevPos[LayerNumber][j] ) Potential = Potential + 0.5 * dq **2. * MassAmu * AmuToKg * ( Frequency[LayerNumber][j] * eVToJoule * 1.e-3 / DiracConstant ) **2. Coordinates[j] = Cartesian[i][j] + dq * 1.e10 print " % 15.13f % 15.13f % 15.13f " %tuple( Coordinates ) else: print " % 15.13f % 15.13f % 15.13f " %tuple( self.Cartesian[i] ) print " \n Velocities ( Ang / fs ): \n " # Print Velocities Velocities = [ 0. , 0., 0. ] Kinetic = 0. for i in range( len( self.Cartesian ) ): LayerNumber = i // ( self.N[0] * self.N[1] ) if LayerNumber <= ( len( Frequency ) - 1 ): for j in range(3): dVq = gauss( 0., StDevVel ) Kinetic = Kinetic + 0.5 * dVq **2. * MassAmu * AmuToKg Velocities[j] = dVq * 1.e-5 print " % 15.13f % 15.13f % 15.13f " %tuple( Velocities ) else: print " % 15.13f % 15.13f % 15.13f " %tuple( [0.,0.,0.] ) print " \n Potential computed for independent HO ( K ): % 15.10f" %( Potential / ( 0.5 * len( Frequency ) * self.N[0] * self.N[1] * 3. * BoltzmannConstant )) print " Kinetic Energy computed for independent HO ( K ): % 15.10f\n" %( Kinetic / ( 0.5 * len( Frequency ) * self.N[0] * self.N[1] * 3. * BoltzmannConstant )) #******************************************************************************************************************** # Calculate inverse of 3x3 matrix def Invert3x3Matrix( self, Matrix ): # First compute the determinant Det = Matrix[ 0][ 0]*( Matrix[ 1][ 1]*Matrix[ 2][ 2] - Matrix[ 1][ 2]*Matrix[ 2][ 1] ) \ - Matrix[ 0][ 1]*( Matrix[ 1][ 0]*Matrix[ 2][ 2] - Matrix[ 1][ 2]*Matrix[ 2][ 0] ) \ + Matrix[ 0][ 2]*( Matrix[ 1][ 0]*Matrix[ 2][ 1] - Matrix[ 1][ 1]*Matrix[ 2][ 0] ) # The inverse matrix is also 3x3.. Inverse = [ [ 0., 0., 0. ], [ 0., 0., 0. ], [ 0., 0., 0. ] ] # Compute element by element for i1 in range(1,4): for i2 in range (1,4): Inverse[ i2-1][ i1-1] = 1.0 / Det * ( Matrix[ i1 %3 ][ i2 %3 ]* \ Matrix[ (i1 + 1)%3 ][ (i2 + 1)%3 ]- \ Matrix[ i1 %3 ][ (i2 + 1)%3 ]* \ Matrix[ (i1 + 1)%3 ][ i2 %3 ] ) return Inverse