''' filter out atoms outside the interval [middle atoms] for passivation use all atoms for counting bonds problematic if they don't fit in unit cell? ''' #--------------------------------------------------------------------------------- def passivate_mol(configuration,CH_len,nearest_neighbor_distance,exclude=[]): ''' Function to passivate (carbon) atoms by adding missing hydrogen atoms. Finds all atoms with only 2 nearest neighbors and attaches a hydrogen atom in such a way that the 4 atoms lie in the same plane, and the C-H bond points in the opposite direction compared to the average of the two nearest neighbor bond vectors. Parameters: configuration : The molecule configuration to be passivated CH_len : The C-H bond length nearest_neighbor_distance : The distance beyond which an atom is not considered a nearest neighbor exclude : A list of atoms (indices) not to be passivated Returns: The passivated molecule configuration ''' import numpy # Extract the original elements and coordinates (peel off units) elements = configuration.elements() c = configuration.cartesianCoordinates().inUnitsOf(Ang) coordinates = configuration.cartesianCoordinates().tolist() # Compute the distance from each atom to all others (order N) distances = [] for i in c: distances.append(numpy.sqrt(((i-c)**2).sum(axis=1))) distances = numpy.array(distances) # Determine which atoms are nearest neighbors, by the given criterion nearest_neighbors = numpy.where((distances<1.1*nearest_neighbor_distance) & (distances>0.9*nearest_neighbor_distance)) # Count how many nearest neighbors each atoms has num_neighbors = [0,]*len(coordinates) for i in range(len(nearest_neighbors[0])): num_neighbors[nearest_neighbors[0][i]] += 1 # Filter out those atoms with only 2 nearest neighbors low_coordination_atoms = numpy.where((numpy.array(num_neighbors)<3))[0].tolist() for i in low_coordination_atoms: # Exclude atoms that are in the "exclude" list if i in exclude or elements[i]!=Carbon: continue # Find neighbors of the low coordination atoms neighbors = nearest_neighbors[1][numpy.where(nearest_neighbors[0]==i)] if len(neighbors)==2: # For atoms with 2 neighbors... # Make neighbor vectors v1 = c[neighbors[0]]-c[i] v2 = c[neighbors[1]]-c[i] # The direction of the C-H bond is the negative of the sum of the two neighbor bonds H_vec = -(v1+v2) / numpy.sqrt(((v1+v2)**2).sum()) # The position of the H atom is in this direction, from the position # of the atom, distance = C-H bond length H_pos = c[i] + H_vec * CH_len # Add new H atom and its coordinates to the atom list coordinates.append(H_pos.tolist()) elements.append(Hydrogen) else: # For atoms with 1 neighbor only, assume that that neighbor has 3 neighbors in total # Find a neighbor of the neighbor to define the bonding plane second_neighbors = nearest_neighbors[1][numpy.where(nearest_neighbors[0]==neighbors[0])].tolist() # Remove the atom "i" itself from the neighbor list to find the two OTHER neighbors second_neighbors.remove(i) for j in range(len(second_neighbors)): # Use these bonds to make the 2 new hydrogen atoms v = c[neighbors[0]]-c[second_neighbors[j]] H_pos = c[i] + v * CH_len / numpy.sqrt((v**2).sum()) # Add new H atom and its coordinates to the atom list coordinates.append(H_pos.tolist()) elements.append(Hydrogen) # Create new molecule and return return MoleculeConfiguration(elements,coordinates*Ang) #--------------------------------------------------------------------------------- def passivate_bulk(configuration,CH_len,nearest_neighbor_distance,periodicA,periodicB,periodicC,exclude=[]): ''' Function to passivate (carbon) atoms by adding missing hydrogen atoms. Finds all atoms with only 2 nearest neighbors and attaches a hydrogen atom in such a way that the 4 atoms lie in the same plane, and the C-H bond points in the opposite direction compared to the average of the two nearest neighbor bond vectors. Parameters: configuration : The molecule configuration to be passivated CH_len : The C-H bond length nearest_neighbor_distance : The distance beyond which an atom is not considered a nearest neighbor exclude : A list of atoms (indices) not to be passivated Returns: The passivated molecule configuration ''' import numpy # Extract the original elements and coordinates (peel off units) elements = numpy.array(configuration.elements()).tolist() num_atoms = len(elements) unitcell = configuration.bravaisLattice() c = configuration.cartesianCoordinates().inUnitsOf(Ang) coordinates = configuration.cartesianCoordinates().tolist() NA,NB,NC = (1,1,1) if periodicA: NA = 3 if periodicB: NB = 3 if periodicC: NC = 3 supercell = configuration.repeat(NA,NB,NC) offset1 = NA*NB*NC/2 offset2 = offset1+1 s = supercell.cartesianCoordinates().inUnitsOf(Ang) super_coordinates = supercell.cartesianCoordinates().tolist() super_elements = configuration.elements() # Compute the distance from each atom to all others (order N) distances = [] for i in s: distances.append(numpy.sqrt(((i-s)**2).sum(axis=1))) distances = numpy.array(distances) # Determine which atoms are nearest neighbors, by the given criterion nearest_neighbors = numpy.where((distances<1.1*nearest_neighbor_distance) & (distances>1e-4)) # Count how many nearest neighbors each atoms has num_neighbors = [0,]*len(super_coordinates) for i in range(len(nearest_neighbors[0])): num_neighbors[nearest_neighbors[0][i]] += 1 # Filter out those atoms with only 2 nearest neighbors all_low_coordination_atoms = numpy.where((numpy.array(num_neighbors)<3))[0] in_original_cell = numpy.where((all_low_coordination_atoms>=num_atoms*offset1) & (all_low_coordination_atoms