# -------------------------------------------------------------
# dope the edges of a GNR ribbbon with specified atom
# -------------------------------------------------------------
def dopeEdges(atomsToInpect, coords, atom, ymax_, ymin_):
# get atoms on top (+y) edge of this electrode
edgeAtoms = []
for ind in range(len(atomsToInpect)):
if (coords[ind])[1] == ymax_: # is it an edge atom?
# we can quickly sort here
newInd = 0
for ind2 in range(len(edgeAtoms)):
# check if this new z coord is > each one in list, stop when not true and insert
if (coords[ind])[2] > (coords[edgeAtoms[ind2]])[2]:
pass
else:
# insert here
newInd = ind2
break
# now insert
#print ymax_, "edge", coords[ind], "inserting at", newInd, str(edgeAtoms)
edgeAtoms.insert(newInd, ind) # save this index
carbonCnt = 0
for ind in range(len(edgeAtoms)):
if atomsToInpect[edgeAtoms[ind]] == Carbon: # we know it is though
carbonCnt += 1
if carbonCnt % dopingLvl == 0:
atomsToInpect[edgeAtoms[ind]] = atom
# get atoms on bottom (-y) edge of this electrode
edgeAtoms = []
for ind in range(len(atomsToInpect)):
if (coords[ind])[1] == ymin_: # is it an edge atom?
# we can quickly sort here
newInd = 0
for ind2 in range(len(edgeAtoms)):
# check if this new z coord is > each one in list, stop when not true and insert
if (coords[ind])[2] > (coords[edgeAtoms[ind2]])[2]:
pass
else:
# insert here
newInd = ind2
break
# now insert
#print ymax_, "edge", coords[ind], "inserting at", newInd, str(edgeAtoms)
edgeAtoms.insert(newInd, ind) # save this index
carbonCnt = 0
for ind in range(len(edgeAtoms)):
if atomsToInpect[edgeAtoms[ind]] == Carbon: # we know it is though
carbonCnt += 1
if carbonCnt % dopingLvl == 0:
atomsToInpect[edgeAtoms[ind]] = atom
return atomsToInpect
example usage: set variable dopingLvl to something like 3, meaning replace every 3 Carbons with whatever you specify in the function call .. first loop through to get max and min y for the Carbon atoms in the structure, then call as below
# find min and max y for Carbon atoms (edge atoms)
yminCarbon = 10000*Angstrom
ymaxCarbon = 0*Angstrom
for ind in range(len(central_region_coordinates)):
if central_region_elements[ind] == Carbon:
if (central_region_coordinates[ind])[1] < yminCarbon:
yminCarbon = (central_region_coordinates[ind])[1]
if (central_region_coordinates[ind])[1] > ymaxCarbon:
ymaxCarbon = (central_region_coordinates[ind])[1]
# NOW DOPE
left_electrode_elements = dopeEdges(left_electrode_elements, left_electrode_coordinates, Nitrogen, ymaxCarbon, yminCarbon)
right_electrode_elements = dopeEdges(right_electrode_elements, right_electrode_coordinates, Nitrogen, ymaxCarbon, yminCarbon)
central_region_elements = dopeEdges(central_region_elements, central_region_coordinates, Boron, ymaxCarbon, yminCarbon)
Oh yes then you also need this part:
# now make sure electrodes are matched inside bulk
print "len(left_electrode_elements)", len(left_electrode_elements)
print "len(right_electrode_elements)", len(right_electrode_elements)
for ind in range(len(left_electrode_elements)):
central_region_elements[ind] = left_electrode_elements[ind]
for ind in range(len(right_electrode_elements)):
central_region_elements[\
len(central_region_elements) - len(right_electrode_elements) + ind] = right_electrode_elements[ind]
;D