Homec4science
Diffusion lammps (efficient_neuronet) Browse

lammps/tools/polybond/lmpsdata.py
3f7d79ab82abefficient_neuronet

Last Change
Raw File

lmpsdata.py
View Options

#
#
# lmpsdata.py
#
# For reading, writing and manipulating lammps data files
# For calculation of certain properities using lammps data files
# For creating VMD input text files using lammps data files
# All x,y,z calculations assume the information includes image flags
class Lmpsdata:
def __init__(self,file,atomtype):
"""initiates lammps data structures"""
self.atomtype=atomtype
self.keywords=[]
self.atoms=[]
self.angles=[]
self.bonds=[]
self.dihedrals=[]
self.dipoles=[]
self.impropers=[]
self.masses=[]
self.shapes=[]
self.velocities=[]
self.anglecoef=[]
self.bondcoef=[]
self.dihedralcoef=[]
self.impropercoef=[]
self.paircoef=[]
self.read(file)
def read(self,file):
"""Reads in lammps data file
Skips the first line of the file (Comment line)
First reads the header portion of the file (sectflag=1)
blank lines are skipped
header keywords delineate an assignment
if no header keyword is found, body portion begins
Second reads the body portion of the file (sectflag=2)
first line of a section has only a keyword
next line is skipped
remaining lines contain values
a blank line signifies the end of that section
if no value is listed on a line than a keyword must be used
File is read until the end
The keywords read in are stored in self.keywords"""
if file=='':
print 'no file is given. Will have to build keywords and class structures manually'
return
sectflag=1
f=open(file,'r')
f.readline()
for line in f:
row=line.split()
if sectflag==1:
if len(row)==0:
#skip line the line is blank
pass
elif len(row)==1:
# Set Sectflag to 2, assume line is a keyword
sectflag=2
checkkey=1 #ensures keyword will be checked in body portion
keyword=row
elif len(row)==2:
if row[1]=='atoms':
self.atomnum=row[0]
self.keywords.append('atoms')
elif row[1]=='bonds':
self.bondnum=row[0]
self.keywords.append('bonds')
elif row[1]=='angles':
self.anglenum=row[0]
self.keywords.append('angles')
elif row[1]=='dihedrals':
self.dihedralnum=row[0]
self.keywords.append('dihedrals')
elif row[1]=='impropers':
self.impropernum=row[0]
self.keywords.append('impropers')
else:
# Set Sectflag to 2, assume line is a keyword
sectflag=2
checkkey=1 #ensures keyword will be checked in body portion
keyword=row
elif len(row)==3:
if row[1]=='atom' and row[2]=='types':
self.atomtypenum=row[0]
self.keywords.append('atom types')
elif row[1]=='bond' and row[2]=='types':
self.bondtypenum=row[0]
self.keywords.append('bond types')
elif row[1]=='angle' and row[2]=='types':
self.angletypenum=row[0]
self.keywords.append('angle types')
elif row[1]=='dihedral' and row[2]=='types':
self.dihedraltypenum=row[0]
self.keywords.append('dihedral types')
elif row[1]=='improper' and row[2]=='types':
self.impropertypenum=row[0]
self.keywords.append('improper types')
else:
# Set Sectflag to 2, assume line is a keyword
sectflag=2
checkkey=1 #ensures keyword will be checked in body portion
keyword=row
elif len(row)==4:
if row[2]=='xlo' and row[3]=='xhi':
self.xdimension=[row[0], row[1]]
self.keywords.append('xlo xhi')
elif row[2]=='ylo' and row[3]=='yhi':
self.ydimension=[row[0], row[1]]
self.keywords.append('ylo yhi')
elif row[2]=='zlo' and row[3]=='zhi':
self.zdimension=[row[0], row[1]]
self.keywords.append('zlo zhi')
else:
# Set Sectflag to 2, assume line is a keyword
sectflag=2
checkkey=1 #ensures keyword will be checked in body portion
keyword=row
elif len(row)==5:
if row[1]=='extra' and row[2]=='bond' and row[3]=='per' and row[4]=='atom':
self.extrabonds=row[0]
self.keywords.append('extra bond per atom')
else:
# Set Sectflag to 2, assume line is a keyword
sectflag=2
checkkey=1 #ensures keyword will be checked in body portion
keyword=row
elif len(row)==6:
if row[3]=='xy' and row[4]=='xz' and row[5]=='yz':
self.tilt=[row[0], row[1], row[2]]
self.keywords.append('xy xz yz')
else:
# Set Sectflag to 2, assume line is a keyword
sectflag=2
checkkey=1 #ensures keyword will be checked in body portion
keyword=row
else:
# set sectflag to 2, assume line is a keyword
sectflag=2
checkkey=1 #ensures keyword will be checked in body portion
keyword=row
elif sectflag==2:
if checkkey==1:
if len(keyword)==1:
if keyword[0]=='Atoms' or keyword[0]=='Velocities' or keyword[0]=='Masses' or\
keyword[0]=='Shapes' or keyword[0]=='Dipoles' or keyword[0]=='Bonds' or\
keyword[0]=='Angles' or keyword[0]=='Dihedrals' or keyword[0]=='Impropers':
bodyflag=1
blanknum=0
self.keywords.append(keyword[0])
checkkey=0
else:
bodyflag=0
checkkey=0
elif len(keyword)==2:
if row[1]=='Coeffs' and (row[0]=='Pair' or row[0]=='Bond' or row[0]=='Angle' or\
row[0]=='Dihedral' or row[0]=='Improper'):#class 2 force field keywords not included
bodyflag=1
blanknum=0
self.keywords.append('{0} {1}'.format(keyword[0],keyword[1]))
checkkey=0
else:
bodyflag=0
checkkey=0
else:
#egnore line and set bodyflag=0
bodyflag=0
checkkey=0
if bodyflag==0: #bodyflag 0 means no body keyword has been found
if len(row)==1:
if row[0]=='Atoms' or row[0]=='Velocities' or row[0]=='Masses' or\
row[0]=='Shapes' or row[0]=='Dipoles' or row[0]=='Bonds' or\
row[0]=='Angles' or row[0]=='Dihedrals' or row[0]=='Impropers':
bodyflag=1
blanknum=0
keyword=row
self.keywords.append(keyword[0])
else:
#egnore line
pass
elif len(row)==2:
if row[1]=='Coeffs' and (row[0]=='Pair' or row[0]=='Bond' or row[0]=='Angle' or\
row[0]=='Dihedral' or row[0]=='Improper'):#class 2 force field keywords not included
bodyflag=1
blanknum=0
keyword=row
self.keywords.append('{0} {1}'.format(keyword[0],keyword[1]))
else:
#egnore line
pass
else:
#egnore line
pass
elif bodyflag==1: #currently assumes 1 or more blank lines are between body data keywords
if len(row)==0:
blanknum+=1
if blanknum>1:
bodyflag=0
elif len(keyword)==1:
if keyword[0]=='Atoms':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.atoms.append(row)
elif keyword[0]=='Velocities':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.velocities.append(row)
elif keyword[0]=='Masses':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.masses.append(row)
elif keyword[0]=='Shapes':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.shapes.append(row)
elif keyword[0]=='Dipoles':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.dipoles.append(row)
elif keyword[0]=='Bonds':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.bonds.append(row)
elif keyword[0]=='Angles':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.angles.append(row)
elif keyword[0]=='Dihedrals':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.dihedrals.append(row)
elif keyword[0]=='Impropers':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.impropers.append(row)
else:
#egnore line and change bodyflag to 0
bodyflag=0
elif len(keyword)==2:
if keyword[0]=='Pair' and keyword[1]=='Coeffs':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.paircoef.append(row)
elif keyword[0]=='Bond' and keyword[1]=='Coeffs':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.bondcoef.append(row)
elif keyword[0]=='Angle' and keyword[1]=='Coeffs':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.anglecoef.append(row)
elif keyword[0]=='Dihedral' and keyword[1]=='Coeffs':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.dihedralcoef.append(row)
elif keyword[0]=='Improper' and keyword[1]=='Coeffs':
try:
int(row[0])
except ValueError:
keyword=row
checkkey=1
else:
self.impropercoef.append(row)
else:
#egnore line and change bodyflag to 0
bodyflag=0
else:
#egnore line and change bodyflag to 0
bodyflag=0
f.close()
def write(self,file,modflag):
"""Write lammps data files using the lammps keywords and lammpsdata structures
writes first line of the file as a Comment line
if no modifications to any of the lammpsdata structures (modflag=0)
Use Keywords to write lammpsdata structures directly
if modifications to any of the lammpsdata structures (modflag=1)
Key Lammpsdata structures like atom numbers, coeficient numbers
need to be modified to match the other modified lammpsdata structures
This section will still use the keywords to write lammpsdata structures.
For all modflags, the code will write data to the file until all of the
keyword's data structures have been finished writing.
The keywords are stored in self.keywords"""
f=open(file,'w')
f.write('polymer data file\n')
for row in self.keywords:
if row=='atoms':
if modflag==0:
f.write('{0} {1}\n'.format(self.atomnum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.atoms),row))
elif row=='bonds':
if modflag==0:
f.write('{0} {1}\n'.format(self.bondnum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.bonds),row))
elif row=='angles':
if modflag==0:
f.write('{0} {1}\n'.format(self.anglenum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.angles),row))
elif row=='dihedrals':
if modflag==0:
f.write('{0} {1}\n'.format(self.dihedralnum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.dihedrals),row))
elif row=='impropers':
if modflag==0:
f.write('{0} {1}\n'.format(self.impropernum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.impropers),row))
elif row=='atom types':
if modflag==0:
f.write('{0} {1}\n'.format(self.atomtypenum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.masses),row))
elif row=='bond types':
if modflag==0:
f.write('{0} {1}\n'.format(self.bondtypenum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.bondcoef),row))
elif row=='angle types':
if modflag==0:
f.write('{0} {1}\n'.format(self.angletypenum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.anglecoef),row))
elif row=='dihedral types':
if modflag==0:
f.write('{0} {1}\n'.format(self.dihedraltypenum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.dihedralcoef),row))
elif row=='improper types':
if modflag==0:
f.write('{0} {1}\n'.format(self.impropertypenum,row))
elif modflag==1:
f.write('{0} {1}\n'.format(len(self.impropercoef),row))
elif row=='xlo xhi':
f.write('{0} {1} {2}\n'.format(self.xdimension[0],self.xdimension[1],row))
elif row=='ylo yhi':
f.write('{0} {1} {2}\n'.format(self.ydimension[0],self.ydimension[1],row))
elif row=='zlo zhi':
f.write('{0} {1} {2}\n'.format(self.zdimension[0],self.zdimension[1],row))
elif row=='extra bond per atom':
f.write('{0} {1}\n'.format(self.extrabonds,row))
elif row=='xy xz yz':
f.write('{0} {1} {2} {3}\n'.format(self.tilt[0],self.tilt[1],self.tilt[2],row))
elif row=='Atoms':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.atoms:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Velocities':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.velocities:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Masses':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.masses:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Shapes':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.shapes:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Dipoles':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.dipoles:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Bonds':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.bonds:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Angles':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.angles:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Dihedrals':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.dihedrals:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
elif row=='Impropers':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.impropers:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Pair Coeffs':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.paircoef:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Bond Coeffs':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.bondcoef:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Angle Coeffs':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.anglecoef:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Dihedral Coeffs':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.dihedralcoef:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
elif row=='Improper Coeffs':
f.write('\n{0}'.format(row)) #new line between header and body portion or two body keywords
f.write('\n') #new line between body keyword and body data
for line in self.impropercoef:
f.write('\n') #creates a new line for adding body data
for item in line:
f.write(' {0}'.format(item)) #adds in each peice of body data with space imbetween
f.write('\n') #allows space to be added between the end of body data and a new keyword
else:
pass
f.close()
def atomorder(self):
"""Takes self.atoms and organizes the atom id from least to greatest.
If the atom ids are allready ordered this algorithm will do nothing."""
current=range(len(self.atoms[0])) # initialize current [assumes self.atoms coloumn #'s does not change]
for i in range(1,len(self.atoms)): #when i=0, self.atoms will not change; therefore its skipped
for k in range(len(self.atoms[i])):
current[k]=self.atoms[i][k]
for j in range(i-1,-1,-1):
for k in range(len(self.atoms[j])):
self.atoms[j+1][k]=self.atoms[j][k]
if int(current[0]) > int(self.atoms[j][0]):
for k in range(len(current)):
self.atoms[j+1][k]=current[k]
break
elif j==0:
for k in range(len(current)):
self.atoms[j][k]=current[k]
#dont need a break here because this is the last j value in the for loop
def addatoms(self, atoms, retlist=False):
"""Appends atoms to self.atoms. Assumes atoms are written in the correct atomtype format.
Change added atom numbers in self.atoms so self.atoms will be in increasing sequential order.
If retlist is True return a list of the modified atom numbers otherwise exit."""
initpos=len(self.atoms) #Store index of first appended atoms
for item in atoms:
self.atoms.append(range(len(item))) #initializing spots for the new atoms to go
numberchange=[] #initiate number change where the changed atom numbers will be stored
count=0
for i in range(initpos,len(self.atoms)):
for j in range(len(self.atoms[i])):
if j==0:
self.atoms[i][0]=str(i+1)
numberchange.append(str(i+1))
else:
self.atoms[i][j]=atoms[count][j]
count+=1
if retlist: return numberchange
else: return #need to redo this algorithm so their is no referencing going on.
def adddata(self,data,keyword):
"""Adds data to a keyword's existing data structure
All body keywords are viable except Atoms which has its own method
All header keywords will do nothing in this algrorithm because they have their own algorithm
changes the body id number so id numbers will be in sequential order"""
if keyword=='Velocities':
pos=len(self.velocities)
for item in data:
self.velocities.append(item) #adds data to existing data structure
self.velocities[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Masses':
pos=len(self.masses)
for item in data:
self.masses.append(item) #adds data to existing data structure
self.masses[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Shapes':
pos=len(self.shapes)
for item in data:
self.shapes.append(item) #adds data to existing data structure
self.shapes[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Dipoles':
pos=len(self.dipoles)
for item in data:
self.dipoles.append(item) #adds data to existing data structure
self.dipoles[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Bonds':
pos=len(self.bonds)
for item in data:
self.bonds.append(item) #adds data to existing data structure
self.bonds[pos][0]=str(pos+1)
pos+=1
elif keyword=='Angles':
pos=len(self.angles)
for item in data:
self.angles.append(item) #adds data to existing data structure
self.angles[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Dihedrals':
pos=len(self.dihedrals)
for item in data:
self.dihedrals.append(item) #adds data to existing data structure
self.dihedrals[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Impropers':
pos=len(self.impropers)
for item in data:
self.impropers.append(item) #adds data to existing data structure
self.impropers[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Pair Coeffs':
pos=len(self.paircoef)
for item in data:
self.paircoef.append(item) #adds data to existing data structure
self.paircoef[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Bond Coeffs':
pos=len(self.bondcoef)
for item in data:
self.bondcoef.append(item) #adds data to existing data structure
self.bondcoef[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Angle Coeffs':
pos=len(self.anglecoef)
for item in data:
self.anglecoef.append(item) #adds data to existing data structure
self.anglecoef[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Dihedral Coeffs':
pos=len(self.dihedralcoef)
for item in data:
self.dihedralcoef.append(item) #adds data to existing data structure
self.dihedralcoef[pos][0]=str(pos+1) #changing body id number
pos+=1
elif keyword=='Improper Coeffs':
pos=len(self.impropercoef)
for item in data:
self.impropercoef.append(item) #adds data to existing data structure
self.impropercoef[pos][0]=str(pos+1) #changing body id number
pos+=1
# def modifiydata(data,keyword): Will not implement
# """for modifying header data"""
def changeatomnum(self,data,originalatoms,atomchanges,vflag=False):
"""Takes data which contains atom ids
and alters those ids from the original atom id system in originalatoms
to the new atom id system in atomchanges."""
#set up boolean array to match the size of data and with all True values
# print atomchanges
array=booleanarray(len(data),len(data[0]),True)
# print originalatoms
if vflag: # for changing atomnumbers for velocities
for i in range(len(originalatoms)): #len of originalatoms should match len of atomchanges
if originalatoms[i][0]==atomchanges[i]: continue
for j in range(len(data)):
for k in range(1):
if data[j][k]==originalatoms[i][0]: #checks if the atom id in data
#matches the atom id in origianalatoms
if array.getelement(j,k): #if the boolean array is true
#set the boolean array to False and
#change the atom id in data to atomchanges
array.setelement(j,k,False)
data[j][k]=atomchanges[i]
break
#the change of boolean array to False insures
#the atom id in data will only be changed once.
else: # for changing atom numbers for everything else
for i in range(len(originalatoms)): #len of originalatoms should match len of atomchanges
if originalatoms[i][0]==atomchanges[i]: continue
for j in range(len(data)):
for k in range(2,len(data[j])):
if data[j][k]==originalatoms[i][0]: #checks if the atom id in data
#matches the atom id in origianalatoms
if array.getelement(j,k): #if the boolean array is true
#set the boolean array to False and
#change the atom id in data to atomchanges
array.setelement(j,k,False)
data[j][k]=atomchanges[i]
break
#the change of boolean array to False insures
#the atom id in data will only be changed once.
# print data
return data
def deletebodydata(self,keyword):
"""Changes a keyword's class data structure to an empty structure"""
if keyword=='Velocities':
self.velocities=[]
elif keyword=='Masses':
self.masses=[]
elif keyword=='Shapes':
self.shapes=[]
elif keyword=='Dipoles':
self.dipoles=[]
elif keyword=='Bonds':
self.bonds=[]
elif keyword=='Angles':
self.angles=[]
elif keyword=='Dihedrals':
self.dihedrals=[]
elif keyword=='Impropers':
self.impropers=[]
elif keyword=='Pair Coeffs':
self.paircoef=[]
elif keyword=='Bond Coeffs':
self.bondcoef=[]
elif keyword=='Angle Coeffs':
self.anglecoef=[]
elif keyword=='Dihedral Coeffs':
self.dihedralcoef=[]
elif keyword=='Improper Coeffs':
self.impropercoef=[]
elif keyword=='Atoms':
self.atoms=[]
def extractmolecules(self,molecule):
"""Takes the variable molecule and
extracts the individual molecules' data back into lmpsdata.
This extraction takes place through a 4 step process.
Step 1: Use a molecule's keywords to alter the lmpsdata data structures to empty.
To acomplish this procedure use the lmpsdata method deletebodydata.
Step 2: Add the molecules' atoms to lmpsdata's atoms using the method addatoms.
Return a list of atom id changes for each molecule.
Step 3: Utilize each molecules list of atom id changes to change their data's atom id numbers.
Uses changeatomnum and returns the altered molecule's data.
Step 4: Add the altered molecules' data to lmpsdata's data using the method adddata"""
#Use molecule index 0's keyword to change the equivalent lmpsdata structures to empty
print 'extracting the molecules back to data'
for keyword in molecule[0].keywords: # step 1
self.deletebodydata(keyword)
atomchanges=[] #initializing step 2
for i in range(len(molecule)): # step 2
atomchanges.append(self.addatoms(molecule[i].atoms,True))
for i in range(len(molecule)): # step 3
for keyword in molecule[i].keywords:
if keyword=='Angles':
molecule[i].angles=self.changeatomnum(molecule[i].angles,molecule[i].atoms,atomchanges[i])
if keyword=='Bonds':
molecule[i].bonds=self.changeatomnum(molecule[i].bonds,molecule[i].atoms,atomchanges[i])
elif keyword=='Dihedrals':
molecule[i].dihedrals=self.changeatomnum(molecule[i].dihedrals,molecule[i].atoms,atomchanges[i])
elif keyword=='Velocities':
molecule[i].velocities=self.changeatomnum(molecule[i].velocities,molecule[i].atoms,atomchanges[i],True)
elif keyword=='Impropers':
molecule[i].impropers=self.changeatomnum(molecule[i].impropers,molecule[i].atoms,atomchanges[i])
for i in range(len(molecule)): # step 4
for keyword in molecule[i].keywords:
if keyword=='Angles':
self.adddata(molecule[i].angles,keyword)
elif keyword=='Bonds':
self.adddata(molecule[i].bonds,keyword)
elif keyword=='Dihedrals':
self.adddata(molecule[i].dihedrals,keyword)
elif keyword=='Velocities':
self.adddata(molecule[i].velocities,keyword)
elif keyword=='Impropers':
self.adddata(molecule[i].impropers,keyword)
def density(self,ringsize,init,final,file=' '):
"""Create spherical shells from the initial radius to the final radius
The spherical shell's thickness is defined by ringsize
Calculate the mass of particles within each spherical shell
Calculate the volume of each spherical shell
Then calculate the density in each spherical shell
Current code is written with the assumption the spheres are centered around the origin
If a file is listed place the results in the variable file otherwise return the results
The distance calculations are written assuming image flags are in the data file"""
rho=[]
r=[init] #initial radius to examine
radius=init+ringsize
while radius<final: #finding the rest of the radii to examine
r.append(radius)
radius+=ringsize
for radius in r:
volume=4.0/3.0*3.1416*((radius+ringsize)**3 - radius**3)
mass=0.0
for row in self.atoms:
if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\
self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\
self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\
self.atomtype=='peri':
l=len(row)-1
dist=float(row[l-3])**2+float(row[l-4])**2+float(row[l-5])**2
if dist<(radius+ringsize)**2 and dist>=radius**2:
if self.atomtype=='atomic' or self.atomtype=='charge' or\
self.atomtype=='colloid' or self.atomtype=='electron' or\
self.atomtype=='granular' or self.atomtype=='peri':
type=int(row[1])-1
else:
type=int(row[2])-1
mass+=float(self.masses[type][1])
#assumes self.masses is written in atomtype order
elif self.atomtype=='dipole':
l=len(row)-1
dist=float(row[l-6])**2+float(row[l-7])**2+float(row[l-8])**2
if dist<(radius+ringsize)**2 and dist>=radius**2:
type=int(row[1])-1
mass+=float(self.masses[type][1])
#assumes self.masses is written in atomtype order
elif self.atomtype=='ellipsoid':
l=len(row)-1
dist=float(row[l-7])**2+float(row[l-8])**2+float(row[l-9])**2
if dist<(radius+ringsize)**2 and dist>=radius**2:
type=int(row[1])
mass+=float(self.masses[type][1])
#assumes self.masses is written in atomtype order
elif self.atomtype=='hybrid':
dist=float(row[2])**2+float(row[3])**2+float(row[4])**2
if dist<(radius+ringsize)**2 and dist>=radius**2:
type=int(row[1])
mass+=float(self.masses[type][1])
#assumes self.masses is written in atomtype order
rho.append(mass/volume)
if file==' ':
return r,rho
else:
f=open(file,'w')
for i in range(len(r)):
f.write('{0} {1}\n'.format(r[i],rho[i]))
f.close()
return
def createxyz(self,file, routine='mass', values=None):
"""Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user.
The mass version is assessed by setting the routine to 'mass' which is the default method.
The other version is assessed by setting the routine to 'atomtype'.
The other version takes values which is a list containing the value the user wants to assign those atomtypes to.
The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1.
This makes it easier to find the atomtype and assign the value from the list
All atom data is assumed to have image flags in the data."""
f=open(file,'w')
f.write('{0}\n'.format(len(self.atoms)))
f.write('atoms\n')
if routine=='mass':
for line in self.atoms:
if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\
self.atomtype=='molecular':
type=int(line[2])-1 #3rd position
else:
type=int(line[1])-1 #2nd position
mass=self.masses[type][1]
if mass=='12.0107': elementnum=6
elif mass=='15.9994': elementnum=8
elif mass=='26.9815':elementnum=13
else:
print 'no matching mass value. The method will exit'
return
l=len(line)-1
if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\
self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\
self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\
self.atomtype=='peri':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3]))
elif self.atomtype=='dipole':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6]))
elif self.atomtype=='ellipsoid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7]))
elif self.atomtype=='hybrid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4]))
f.close()
elif routine=='atomtype':
for line in self.atoms:
if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\
self.atomtype=='molecular':
type=int(line[2])-1 #3rd position
else:
type=int(line[1])-1 #2nd position
elementnum=values[type]
l=len(line)-1
if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\
self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\
self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\
self.atomtype=='peri':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3]))
elif self.atomtype=='dipole':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6]))
elif self.atomtype=='ellipsoid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7]))
elif self.atomtype=='hybrid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4]))
f.close()
class booleanarray:
"""A class that stores boolean values in a list of lists."""
def __init__(self,rownum, colnum, initval):
""" initialise a list of lists (array) with
rownum correspondinig to the number of lists in the list and
colnum corresponding to the number of elements in the list's list.
initval is the value the list of lists will be initialized with.
initval should be a boolean value."""
# initializing self.array
self.array=[]
for i in range(rownum):
self.array.append(range(colnum))
# setting elements of self.array to initval
for i in range(rownum):
for j in range(colnum):
self.setelement(i,j,initval)
def setelement(self,rownum, colnum, value):
"""Assigns value to the list of lists (array) element at rownum and colnum."""
self.array[rownum][colnum]=value
def getelement(self,rownum,colnum):
"""Returns element in the list of lists (array) at rownum and colnum."""
return self.array[rownum][colnum]
def atomdistance(a,b,atomtype):
"""Returns the distance between atom a and b.
Atomtype is the style the atom data is written in.
All atom data is assumed to have image flags in the data.
Atom a and atom b are assumed to be the same atomtype."""
from math import sqrt
if atomtype=='angle' or atomtype=='atomic' or atomtype=='bond' or\
atomtype=='charge' or atomtype=='colloid' or atomtype=='electron' or\
atomtype=='full' or atomtype=='granular' or atomtype=='molecular' or\
atomtype=='peri':
l=len(a)-1
dist=sqrt((float(a[l-3])-float(b[l-3]))**2\
+(float(a[l-4])-float(b[l-4]))**2\
+(float(a[l-5])-float(b[l-5]))**2)
elif atomtype=='dipole':
l=len(a)-1
dist=sqrt((float(a[l-6])-float(b[l-6]))**2\
+(float(a[l-7])-float(b[l-7]))**2\
+(float(a[l-8])-float(b[l-8]))**2)
elif atomtype=='ellipsoid':
l=len(a)-1
dist=sqrt((float(a[l-7])-float(b[l-7]))**2\
+(float(a[l-8])-float(b[l-8]))**2\
+(float(a[l-9])-float(b[l-9]))**2)
elif atomtype=='hybrid':
dist=sqrt((float(a[2])-float(b[2]))**2\
+(float(a[3])-float(b[3]))**2\
+(float(a[4])-float(b[4]))**2)
return dist
def distance(a,coord,atomtype):
"""Returns the distance between atom a and coord.
Atomtype is the style the atom data is written in.
All atom data is assumed to have image flags in the data."""
from math import sqrt #need to alter this slightly
if atomtype=='angle' or atomtype=='atomic' or atomtype=='bond' or\
atomtype=='charge' or atomtype=='colloid' or atomtype=='electron' or\
atomtype=='full' or atomtype=='granular' or atomtype=='molecular' or\
atomtype=='peri':
l=len(a)-1
dist=sqrt((float(a[l-3])-coord[2])**2\
+(float(a[l-4])-coord[1])**2\
+(float(a[l-5])-coord[0])**2)
elif atomtype=='dipole':
l=len(a)-1
dist=sqrt((float(a[l-6])-coord[2])**2\
+(float(a[l-7])-coord[1])**2\
+(float(a[l-8])-coord[0])**2)
elif atomtype=='ellipsoid':
l=len(a)-1
dist=sqrt((float(a[l-7])-coord[2])**2\
+(float(a[l-8])-coord[1])**2\
+(float(a[l-9])-coord[0])**2)
elif atomtype=='hybrid':
dist=sqrt((float(a[2])-coord[0])**2\
+(float(a[3])-coord[1])**2\
+(float(a[4])-coord[0])**2)
return dist
class particlesurface:
def __init__(self,particle,cutoff,atomid,atomtype,shape='sphere'):
"""Builds a particle surface with a specific shape from a particle
The atoms choosen from the surface will have the specific atomid
atomid will be given in terms of an integer and not a string."""
self.particle=particle.atoms
self.atomtype=atomtype
self.cutoff=cutoff
self.surface=[]
self.particlelocator=[] #stores surface particle locations with respect to the particle
self.deletedatoms=[]
if shape=='sphere': self.createspheresurf(atomid)
def createspheresurf(self,atomid):
"""Assumes sphere is centered around 0,0,0.
Finds maximum distance particle with atomid.
Than all atoms within cutoff distance of max distance
will be included into self.surface
Also will build a list of where the surface particles are
in relationship to the particle.
Will create a booleanarray to keep track of any changes made to the surface."""
particledist=[]
for row in self.particle: #Stores all of the distances of the particle in particledist
if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\
self.atomtype=='molecular':
type=int(row[2]) #3rd position
else:
type=int(row[1]) #2nd position
if type==atomid: #only calculates the distance if the atom has the correct atomid number
particledist.append(distance(row,[0,0,0],self.atomtype))
else:
particledist.append(0.0)
self.maxdist=max(particledist)# finds the max dist.
for i in range(len(particledist)):
if particledist[i]>=(self.maxdist-self.cutoff):
self.particlelocator.append(i)
self.surface.append(self.particle[i])
self.bool=booleanarray(len(self.surface),1,True)
def getsurfatom(self,rownum):
return self.surface[rownum]
def getbool(self,rownum):
return self.bool.getelement(rownum,0)
def setbool(self,rownum,value):
self.bool.setelement(rownum,0,value)
def removesurfatom(self,rownum):
"""note this does not actually remove the surface atom.
Instead this adds a value to the list deletedatoms.
This list is later used in extractparticle to actually delete those atoms"""
self.deletedatoms.append(self.particlelocator[rownum])
def extractparticle(self):
""""deletesatoms from particle from the list of deletedatoms
and than returns the particle."""
self.particle=self.deleteatoms(self.particle,self.deletedatoms)
return self.particle
def deleteatoms(self,structure,rows):
"""delete atoms from particle and shifts the structure down"""
new=[]
#mulitple copying of b to the rows being replaced.
if rows==[]:
for line in structure:
new.append(line) #if no rows need replacing copy structure
return new
for i in range(rows[0]):
new.append(structure[i]) #copy structure to new until the first replaced row
count=0
for i in range(rows[0],len(structure)-len(rows)):# to replace rows and shift undeleted rows over
for j in range(i+1+count,len(structure)):
for val in rows:
if val==j:
count+=1
break
if val==j:continue
else:
new.append(structure[j])
break
return new
def addatom(self,atomtype,charge=None,moleculenum=None):
"""Adds an atom to the particle surface between maxdist and the cutoff distance.
This atom is stored in self.particle and self.surface.
In the current coding the particle is centered at (0,0,0).
This method has a required input of atomtype and optional inputs of charge and moleculenum.
This method currently does not support correctly self.atomtype of dipole, electron, ellipsoid, granular, peri or hybrid
Will use the image flags 0,0,0.
Note: The atomtype here is different from the atomtype attribute as part of this class
Note: If the added atom will be required for bonding later than you will need to extract the particle data
and rebuild the surface, because this method doesn't include updates to the required variables due to some coding issues"""
import math, random
# Checks to make sure self.atomtype is not set to an unsupported value
if self.atomtype=='dipole' or self.atomtype=='electron' or\
self.atomtype=='ellipsoid' or self.atomtype=='peri' or self.atomtype=='hybrid':
print self.atomtype, 'is not currently supported. But, you can add support by adding the needed functionality and inputs'
return
print 'adding an atom to the surface'
# Randomly assigns the position of a new atom in a spherical shell between self.cutoff and self.maxdist
foundpoint=False
while (foundpoint==False):
x=random.uniform(-self.maxdist,self.maxdist)
y=random.uniform(-self.maxdist,self.maxdist)
try:
z=random.uniform(math.sqrt(self.cutoff**2-x**2-y**2),math.sqrt(self.maxdist**2-x**2-y**2))
except ValueError:
continue
distance=math.sqrt(x**2+y**2+z**2)
if distance<self.maxdist and distance>=self.maxdist-self.cutoff:
foundpoint=True
# initialize new row in self.particle and self.surface
self.particle.append([])
self.surface.append([])
pposition=len(self.particle)-1 #particle postion
sposition=len(self.surface)-1 #surface position
# Adds the atom id number to the new row
self.particle[pposition].append(str(pposition+1))
self.surface[sposition].append(str(pposition+1))
# Adds the atomtype and the moleculenum if needed to the new row
if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\
self.atomtype=='molecular':
self.particle[pposition].append(str(moleculenum))
self.surface[sposition].append(str(moleculenum))
self.particle[pposition].append(str(atomtype))
self.surface[sposition].append(str(atomtype))
else:
self.particle[pposition].append(str(atomtype))
self.surface[sposition].append(str(atomtype))
# Adds the charge if needed to the new row
if self.atomtype=='charge' or self.atomtype=='full':
self.particle[pposition].append(str(charge))
self.surface[sposition].append(str(charge))
# Adds the atom's position to the new row
self.particle[pposition].append(str(x))
self.surface[sposition].append(str(x))
self.particle[pposition].append(str(y))
self.surface[sposition].append(str(y))
self.particle[pposition].append(str(z))
self.surface[sposition].append(str(z))
# Adds the atom's image flags to the new row
self.particle[pposition].append('0')
self.surface[sposition].append('0')
self.particle[pposition].append('0')
self.surface[sposition].append('0')
self.particle[pposition].append('0')
self.surface[sposition].append('0')
def createxyz(self,file,data,routine='mass', values=None):
"""This shows the particle surface. To show the particle surface after bonding has occured,
you will need to extract the particle than reinsert the particle into the class and use createxyz.
Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user.
The mass version is assessed by setting the routine to 'mass' which is the default method.
The other version is assessed by setting the routine to 'atomtype'.
The other version takes values which is a list containing the value the user wants to assign those atomtypes to.
The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1.
This makes it easier to find the atomtype and assign the value from the list
All atom data is assumed to have image flags in the data."""
f=open(file,'w')
f.write('{0}\n'.format(len(self.surface)))
f.write('atoms\n')
if routine=='mass':
for line in self.surface:
if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\
self.atomtype=='molecular':
type=int(line[2])-1 #3rd position
else:
type=int(line[1])-1 #2nd position
mass=data.masses[type][1]
if mass=='12.0107': elementnum=6
elif mass=='15.9994': elementnum=8
elif mass=='26.981539':elementnum=13
else:
print 'no matching mass value. The method will exit'
return
l=len(line)-1
if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\
self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\
self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\
self.atomtype=='peri':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3]))
elif self.atomtype=='dipole':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6]))
elif self.atomtype=='ellipsoid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7]))
elif self.atomtype=='hybrid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4]))
f.close()
elif routine=='atomtype':
for line in self.surface:
if self.atomtype=='angle' or self.atomtype=='bond' or self.atomtype=='full' or\
self.atomtype=='molecular':
type=int(line[2])-1 #3rd position
else:
type=int(line[1])-1 #2nd position
elementnum=values[type]
l=len(line)-1
if self.atomtype=='angle' or self.atomtype=='atomic' or self.atomtype=='bond' or\
self.atomtype=='charge' or self.atomtype=='colloid' or self.atomtype=='electron' or\
self.atomtype=='full' or self.atomtype=='granular' or self.atomtype=='molecular' or\
self.atomtype=='peri':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3]))
elif self.atomtype=='dipole':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-8],line[l-7],line[l-6]))
elif self.atomtype=='ellipsoid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-9],line[l-8],line[l-7]))
elif self.atomtype=='hybrid':
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[2],line[3],line[4]))
f.close()
def molecules(data,init,final,processors, method='all'):
""" There are two ways to initialize the class Lmpsmolecule.
Method controls which way is chosen.
Acceptable values for Method are 'all', 'atom'.
The default method is 'all'"""
molecule=[]
from multiprocessing import Pool
p=Pool(processors)
for i in range(init,final+1):
molecule.append(p.apply_async(Lmpsmolecule,(i,data,method,)))
for i in range(len(molecule)):
molecule[i]=molecule[i].get()
p.close()
p.join()
return molecule
class Lmpsmolecule: #Technically should be a meta class but written as a seperate class for easier coding.
def __init__(self,moleculenum,data,method):
"""initiates lammps molecule structures
and than extract the appropriate molecular structures from the base class data"""
# ***** Lammps Molecule Structures
self.keywords=['Atoms']#include atoms here since this will be in every instance of this class
self.atoms=[] #extract
self.angles=[] #extract if keyword is there
self.bonds=[] #extract if keyword is there
self.dihedrals=[] #extract if keyword is there
self.impropers=[] #extract if keyword is there
self.velocities=[] #extract if keyword is there
# ***** Molecule number
self.moleculenum=moleculenum
# ***** Extracting the moleculue's atom information from data
self.extract(data,method)
# ***** If methods is 'all' then self.extract will extract all of the molecule's information from data
def extract(self,data,method):
"""uses base class data to extract the molecules atoms
and any other molecule data from molecule moleculenum. This extraction is done in a two step process.
Step 1: extract data.atoms with moleculenum and renumber self.atoms beginning from 1
store extracted atom numbers from data.atoms in a list called temp
Step 2: pass temp and data to method changeatomnum
which extracts the rest of the Lammps Molecule structures from data
and changes the atomnumbers in those structures to match the ones already in self.atoms"""
#checking to make sure atomtype is a valid molecule
#if not a valid molecule print error and exit method
if data.atomtype=='full' or data.atomtype=='molecular':
# ***** Sets the molecule's atomtype
self.atomtype=data.atomtype
else:
print "not a valid molecular structure"
return
#extract the molecule self.moleculenum from data.atom to self.atom
atomnum=1
temp=[]
for i in range(len(data.atoms)):
if int(data.atoms[i][1])==self.moleculenum:
temp.append(data.atoms[i][0]) #store extracted atomnumbers into temp
self.atoms.append(data.atoms[i]) #store extracted atom into self.atom
self.atoms[atomnum-1][0]=str(atomnum) #change extracted atom to the correct atomnumber
atomnum+=1 #update atomnumber
#extract the rest of the Lammps Molecule structures from data using changeatomnum
if method=='atom':
return
else:
self.changeatomnum(temp,data)
def changeatomnum(self,temp,data=' '):
"""changes the atomnumbers in Lmpsmolecule structures angles, bonds, dihedrals, impropers,
and velocities in a three step process.
If data is defined extract the rest of the keywords used for Lmpsmolecule.
Step 1A:If data is defined copy from data one of the above structures.
If data is not defined skip this step.
Step 1B:If data is defined remove the rows of copy which do not contain any values from temp
Step 2: If data is defined extract from copy to an above Lmpsmolecule structure and remove the extracted rows
If data is not defined skip this step.
Step 3: Use temp and self.atoms to convert the atomnumbers in Lmpsmolecule structures
to the correct values.
temp and self.atoms line up, so temp[i] and self.atoms[i][0]
correspond to the current Lmpsmolecule structures atom numbers
and correct values
will use booleanarray class to ensure that lmpsmolecule's structure values are altered only once."""
#Step 1 and Step 2
if data!=' ':
#extract molecular keywords from data.keywords
for keywords in data.keywords:
if keywords=='Angles' or keywords=='Bonds' or keywords=='Dihedrals' or\
keywords=='Impropers' or keywords=='Velocities':
self.keywords.append(keywords)
for keywords in self.keywords:
if keywords!='Atoms': print 'extracting the data from', keywords
if keywords=='Angles':
copy=self.copy(data.angles) #Step 1A
dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep
for j in range(len(copy)): #Step 1B
dnr.append(0) #adds 0 to all list elements of dnr
for item in temp: #finds copied structure that has temp values
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
dnr[j]=1 #changes jth list element of dnr to 1
break
remove=[]
for j in range(len(dnr)): # finds dnr values that are still 0
if dnr[j]==0:
remove.append(j) #and appends their index value to the list remove
copy=self.deleterows(copy,remove) #removes all unneeded rows from copy
structnum=1
print 'the length of data is', len(copy)
for item in temp: #Step 2
found=[]
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
found.append(j)
self.angles.append(copy[j])
l=len(self.angles)-1
self.angles[l][0]=str(structnum)
structnum+=1
break
# print 'the item is', item
# print 'found is', found
copy=self.deleterows(copy,found)
elif keywords=='Bonds':
copy=self.copy(data.bonds) #Step 1A
dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep
for j in range(len(copy)): #Step 1B
dnr.append(0) #adds 0 to all list elements of dnr
for item in temp: #finds copied structure that has temp values
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
dnr[j]=1 #changes jth list element of dnr to 1
break
remove=[]
for j in range(len(dnr)): # finds dnr values that are still 0
if dnr[j]==0:
remove.append(j) #and appends their index value to the list remove
copy=self.deleterows(copy,remove) #removes all unneeded rows from copy
structnum=1
print 'the length of data is', len(copy)
for item in temp: #Step 2
found=[]
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
found.append(j)
self.bonds.append(copy[j])
l=len(self.bonds)-1
self.bonds[l][0]=str(structnum)
structnum+=1
break
# print 'the item is', item
# print 'found is', found
copy=self.deleterows(copy,found)
elif keywords=='Dihedrals':
copy=self.copy(data.dihedrals) #Step 1A
dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep
for j in range(len(copy)): #Step 1B
dnr.append(0) #adds 0 to all list elements of dnr
for item in temp: #finds copied structure that has temp values
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
dnr[j]=1 #changes jth list element of dnr to 1
break
remove=[]
for j in range(len(dnr)): # finds dnr values that are still 0
if dnr[j]==0:
remove.append(j) #and appends their index value to the list remove
copy=self.deleterows(copy,remove) #removes all unneeded rows from copy
structnum=1
print 'the length of data is', len(copy)
structnum=1
for item in temp: #Step 2
found=[]
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
found.append(j)
self.dihedrals.append(copy[j])
l=len(self.dihedrals)-1
self.dihedrals[l][0]=str(structnum)
structnum+=1
break
# print 'the item is', item
# print 'found is', found
copy=self.deleterows(copy,found)
elif keywords=='Impropers':
copy=self.copy(data.impropers) #Step 1B
dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep
for j in range(len(copy)): #Step 1B
dnr.append(0) #adds 0 to all list elements of dnr
for item in temp: #finds copied structure that has temp values
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
dnr[j]=1 #changes jth list element of dnr to 1
break
remove=[]
for j in range(len(dnr)): # finds dnr values that are still 0
if dnr[j]==0:
remove.append(j) #and appends their index value to the list remove
copy=self.deleterows(copy,remove) #removes all unneeded rows from copy
structnum=1
print 'the length of data is', len(copy)
for item in temp: #Step 2
found=[]
for j in range(len(copy)):
for i in range(2,len(copy[j])):
if copy[j][i]==item:
found.append(j)
self.impropers.append(copy[j])
l=len(self.impropers)-1
self.impropers[l][0]=str(structnum)
structnum+=1
break
# print 'the item is', item
# print 'found is', found
copy=self.deleterows(copy,found)
elif keywords=='Velocities':
copy=self.copy(data.velocities) #Step 1
dnr=[] #dnr means do not remove a 0 means remove and a 1 means keep
for j in range(len(copy)): #Step 1B
dnr.append(0) #adds 0 to all list elements of dnr
for item in temp: #finds copied structure that has temp values
for j in range(len(copy)):
for i in range(1):
if copy[j][i]==item:
dnr[j]=1 #changes jth list element of dnr to 1
break
remove=[]
for j in range(len(dnr)): # finds dnr values that are still 0
if dnr[j]==0:
remove.append(j) #and appends their index value to the list remove
copy=self.deleterows(copy,remove) #removes all unneeded rows from copy
structnum=1
print 'the length of data is', len(copy)
for item in temp: #Step 2
found=[]
for j in range(len(copy)):
for i in range(1):
if copy[j][i]==item:
found.append(j)
self.velocities.append(copy[j])
l=len(self.velocities)-1
self.velocities[l][0]=str(structnum)
structnum+=1
break
# print 'the item is', item
# print 'found is', found
copy=self.deleterows(copy,found)
#Step 3
for keywords in self.keywords:
if keywords!='Atoms': print 'altering data structure values for', keywords
if keywords=='Angles':
copy=self.copy(self.angles)
array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values
for i in range(len(temp)):
if temp[i]==self.atoms[i][0]: continue
for j in range(len(copy)):
for k in range(2,len(copy[j])):
if copy[j][k]==temp[i]:
if array.getelement(j,k): #if the boolean array is true
self.angles[j][k]=self.atoms[i][0]
array.setelement(j,k,False)
break
elif keywords=='Bonds':
copy=self.copy(self.bonds)
array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values
for i in range(len(temp)):
if temp[i]==self.atoms[i][0]: continue
for j in range(len(copy)):
for k in range(2,len(copy[j])):
if copy[j][k]==temp[i]:
if array.getelement(j,k): #if the boolean array is true
self.bonds[j][k]=self.atoms[i][0]
array.setelement(j,k,False)
break
elif keywords=='Dihedrals':
copy=self.copy(self.dihedrals)
array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values
for i in range(len(temp)):
if temp[i]==self.atoms[i][0]: continue
for j in range(len(copy)):
for k in range(2,len(copy[j])):
if copy[j][k]==temp[i]:
if array.getelement(j,k): #if the boolean array is true
self.dihedrals[j][k]=self.atoms[i][0]
array.setelement(j,k,False)
break
elif keywords=='Impropers':
copy=self.copy(self.impropers)
array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values
for i in range(len(temp)):
if temp[i]==self.atoms[i][0]: continue
for j in range(len(copy)):
for k in range(2,len(copy[j])):
if copy[j][k]==temp[i]:
if array.getelement(j,k): #if the boolean array is true
self.impropers[j][k]=self.atoms[i][0]
array.setelement(j,k,False)
break
elif keywords=='Velocities':
copy=self.copy(self.velocities)
array=booleanarray(len(copy),len(copy[0]),True) #creating a booleanarray with true values
for i in range(len(temp)):
if temp[i]==self.atoms[i][0]: continue
for j in range(len(copy)):
for k in range(1):
if copy[j][k]==temp[i]:
if array.getelement(j,k): #if the boolean array is true
self.velocities[j][k]=self.atoms[i][0]
array.setelement(j,k,False)
break
def copy(self,structure):
"""copies structure to same and returns same."""
same=[]
for row in structure:
same.append(row)
return same
def deleterows(self,structure,rows): # run through this {structure is list of strings and rows is list
#of numbers}
"""delete rows in a structure and shifts the structure up
rows must be in increasing order for this algorithm to work correctly"""
new=[]
#mulitple copying of b to the rows being replaced.
if rows==[]:
for line in structure:
new.append(line) #if no rows need replacing copy structure
return new
for i in range(rows[0]):
new.append(structure[i]) #copy structure to new until the first replaced row
count=0
for i in range(rows[0],len(structure)-len(rows)):# to replace rows and shift undeleted rows over
for j in range(i+1+count,len(structure)):
for val in rows:
if val==j:
count+=1
break
if val==j:continue
else:
new.append(structure[j])
break
return new
def modifyatom(self,atomnumber,column,newvalue):
"""modifies self.atom[atomnumber-1][column] to newvalue.
*note: newvalue needs to be a string
if column is 0 than changeatomnum method might need runnning for all atoms
that have had their atomnumber(column=0) modified.
changeatomnum method is not ran in this method when column=0"""
self.atoms[atomnumber-1][column]=newvalue
def deleteatoms(self,atomnumbers,atomid):
"""Algorithm to find all atoms bonded to atomnumbers in the direction of atomid.
Atomid can be a list or a single integer. The single integer corresponds to atom's atomtype value.
The list corresponds to the atom's atomid values. The only difference between both cases is in the top portion of code.
When all bonded atoms have been found; ie: (the modified return values from findbonds yields an empty list);
delete those bonded atoms and all molecule structures which contain those atoms.
Convert the list of bonded atoms into rows and delete those rows from molecule.atoms"""
print 'finding atoms to delete'
bondedatoms=[]
# 1st iteration
nextatoms=self.findbonds(atomnumbers)
try: # tests whether atomid is a list or a single integer
atomid[0]
except TypeError:
#need to remove atoms from nextatoms which dont have the proper atom id
testflag=True
i=0
while testflag:
if i>len(nextatoms)-1: break #For the case where i becomes greater than the len of nextatoms break loop
if int(self.atoms[nextatoms[i]-1][2])!=atomid:#uses the atom id for the row in atoms and than checks atom id
del nextatoms[i] #delets the atom at i
i-=1 #and than decreases i by 1 so next atom in the list will line up when i is increased
i+=1 #increase i by 1
else:
#need to remove atoms from nextatoms which dont have the proper atom id
testflag=True
i=0
while testflag:
if i>len(nextatoms)-1: break #For the case where i becomes greater than the len of nextatoms break loop
keep=False
for id in atomid:
if int(self.atoms[nextatoms[i]-1][0])==id: #checking if atomid is in next atom
keep=True #keep this atom
break
if not keep:
del nextatoms[i] #delets the atom at i
i-=1 #and than decreases i by 1 so next atom in the list will line up when i is increased
i+=1 #increase i by 1
#append next atoms into bondedatoms
#copy next atoms into prevatoms
prevatoms=[]
for atom in nextatoms:
bondedatoms.append(atom)
prevatoms.append(atom)
#2nd iteration
if prevatoms==[]:
print 'no bonds were found in first iteration that had atomid criteria'
return
nextatoms=self.findbonds(prevatoms)
#need to remove atoms from nextatoms which are in atomnumbers
for atom in atomnumbers:
for i in range(len(nextatoms)):
if nextatoms[i]==atom: #checking if atom is in next atom
del nextatoms[i] #delete the atom at i
break
if nextatoms==[]: break #all bonds from find bonds have allready been added to bondedatoms
#append next atoms into bondedatoms
#copy next atoms into prevatoms
prevatoms=[]
for atom in nextatoms:
bondedatoms.append(atom)
prevatoms.append(atom)
#iterative proccess for finding the rest of the atoms bonded to atomnumbers in the direction of atomid.
while prevatoms!=[]:
nextatoms=self.findbonds(prevatoms)
#need to remove atoms from nextatoms which are in the prevatoms
for atom in bondedatoms:
for i in range(len(nextatoms)):
if nextatoms[i]==atom: #checking if atom is in next atom
del nextatoms[i] #delete the atom at i
break
if nextatoms==[]: break #all bonds from find bonds have allready been added to bondedatoms
#append next atoms into bondedatoms
#copy next atoms into prevatoms
prevatoms=[]
for atom in nextatoms:
bondedatoms.append(atom)
prevatoms.append(atom)
print 'the atoms to delete are', bondedatoms
print 'deleting atoms from structures'
#delete bonded atoms from the molecule's structures except the atom structure
for i in range(1,len(self.keywords)): #goes through all keywords except the atom keyword
print self.keywords[i]
if self.keywords[i]=='Angles':
rows=self.findatomnumbers(bondedatoms,self.angles,False)
rows=self.listorder(rows) #to order the rows in increasing order
self.angles=self.deleterows(self.angles,rows) #requires that rows be in increasing order
elif self.keywords[i]=='Bonds':
rows=self.findatomnumbers(bondedatoms,self.bonds,False)
rows=self.listorder(rows) #to order the rows in increasing order
self.bonds=self.deleterows(self.bonds,rows)#requires that rows be in increasing order
elif self.keywords[i]=='Dihedrals':
rows=self.findatomnumbers(bondedatoms,self.dihedrals,False)
rows=self.listorder(rows) #to order the rows in increasing order
self.dihedrals=self.deleterows(self.dihedrals,rows) #requires that rows be in increasing order
elif self.keywords[i]=='Impropers':
rows=self.findatomnumbers(bondedatoms,self.impropers,False)
rows=self.listorder(rows) #to order the rows in increasing order
self.impropers=self.deleterows(self.impropers,rows) #requires that rows be in increasing order
elif self.keywords[i]=='Velocities':
rows=self.findatomnumbers(bondedatoms,self.velocities,True)
rows=self.listorder(rows) #to order the rows in increasing order
self.velocities=self.deleterows(self.velocities,rows) #requires that rows be in increasing order
print 'Atoms'
#convert bondedatoms from atom numbers to row numbers
for i in range(len(bondedatoms)):
bondedatoms[i]-=1
bondedatoms=self.listorder(bondedatoms) #to order the row numbers in increasing order
#delete bonded atoms (row numbers) from the atom structure
self.atoms=self.deleterows(self.atoms,bondedatoms) #requires that row numbers be in increasing order
def findatomnumbers(self,atomnumbers,structure,vflag): #need to read through this algorithm
"""Algorithm to find atomnumbers in a molecule structure except.
Atoms structure is not handled in here.
Returns a list of the rows in which the atomnumbers are contained in the molecule structure"""
rows=[]
if vflag: #for handling velocity structure
for atom in atomnumbers:
for i in range(len(structure)):
if int(structure[i][0])==atom:
#duplicate rows for vflag=True are not possible.
rows.append(i)
break
else: #for handling all other structures except atoms
for atom in atomnumbers:
for i in range(len(structure)):
for j in range(2,len(structure[i])):
if int(structure[i][j])==atom:
#need to make sure duplicate value of rows are not being added
if rows==[]:
rows.append(i) #appends the row number
else:
#checking for duplicate values of rows
duplicateflag=0
for k in range(len(rows)):
if rows[k]==i:
duplicateflag=1
break
if duplicateflag==0: #if no duplicates adds bond number
rows.append(i) #appends the row number
break
print 'the finished row is', rows
return rows
def listorder(self,struct):
"""Takes struct and organizes the list from least to greatest.
If the list is allready ordered this algorithm will do nothing."""
if len(struct)==1: return struct #with the length at 1; thier is only one element and theirfore nothing to order
for i in range(1,len(struct)): #when i=0, struct will not change; therefore its skipped
copy=struct[i]
for j in range(i-1,-1,-1):
struct[j+1]=struct[j]
if copy >struct[j]:
struct[j+1]=copy
break
elif j==0:
struct[j]=copy
#dont need a break here because this is the last j value in the for loop
print 'the organized row is', struct
return struct
def findbonds(self,atomnumbers):
"""Algorithm to find all atoms bonded to atomnumbers.
Returns a list of atomnumbers"""
#finds the bonds in which the atomnumbers are located
bondids=[]
for i in range(len(atomnumbers)):
for j in range(len(self.bonds)):
for k in range(2,len(self.bonds[j])):
if int(self.bonds[j][k])==atomnumbers[i]:
if bondids==[]:
bondids.append(int(self.bonds[j][0])) #appends the bond number
else:
#checking for duplicates of bondids
duplicateflag=0
for l in range(len(bondids)):
if bondids[l]==int(self.bonds[j][0]):
duplicateflag=1
break
if duplicateflag==0: #if no duplicates adds bond number
bondids.append(int(self.bonds[j][0]))
break
#Using bondids find the atoms bonded to atomnumbers
bondedatoms=[]
from math import fabs
for id in bondids:
for atoms in atomnumbers:
found=False
for i in range(2,len(self.bonds[id-1])):
if int(self.bonds[id-1][i])==atoms:
j=int(fabs(i-5)) #switches the index from the atomnumber location to the other location
bondedatoms.append(int(self.bonds[id-1][j])) #appends the atomnuber at the other location
found=True
break
if found==True: break
return bondedatoms
def findparticlebondingpoints(self,particle,atomid,cutoffdistance,bondnumber):
"""Particle is a particlesurfaceobject, atomid is an int.
Finds atoms with atomid in the molecule which are less than the cutoffdistance from atoms in the particle.
The found atoms and particles locations are stored in a 3 dimensional list called possiblebonds
The first index corresponds with the molecule's atom
The second index correspons with the molecule/particle combination
The third index corresponds with whether the value is the molecule or the particle
Always bonds the two ends of the molecule that meet cutoff requirement.
All other possible bonds are randomly choosen until the required bondnumbers are met.
After every bond is choosen, the particle object's boolean list is updated, and possiblebonds is updated.
The update to possiblebonds involves removing the row from which the bonded molecule's atom is located
and also removing the particle atom and it's corresponding bonded atom from other rows of possiblebonds.
The final bonds are all stored in self.bonding as a 2 dimensional list"""
possiblebonds=[]
for i in range(len(self.atoms)):
if int(self.atoms[i][2])==atomid:
row=[]
for j in range(len(particle.surface)):
#assumes molecule and particle have same atomtype if not true than atomdistance will give bad value
if atomdistance(self.atoms[i],particle.surface[j],self.atomtype)<=cutoffdistance:
if particle.getbool(j): row.append([i,j]) #if not bonded than can form a possible bond
if row!=[]:possiblebonds.append(row) #need to correct this....
#initiate section which assigns bonds into bonding information
bonds=0
self.bondinginformation=[] #initiate new class member
#Checks to see if no bonds are possible [2 possible cases]
if possiblebonds==[]:
print 'no possible bonds can be formed'
return
if bondnumber==0:
print 'bondnumber is 0; so, no possible bonds can form'
return
#section which assigns a bond to the first molecule atom which can be bonded.
self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,0))
bonds+=1
if bonds==bondnumber: return
del possiblebonds[0] #deletes possible bonds to the first molecule which can be bonded
l=len(self.bondinginformation)-1
possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1]) #updates possiblebonds
#by removing any bonds which contain the newly bonded particle.
if possiblebonds==[]:return
#section which finds the last molecule atom which can be bonded
l=len(possiblebonds)-1
while possiblebonds[l]==[]:# to find the last molecule atom which can be bonded
if possiblebonds==[]:return
del possiblebonds[l] #since there are no more possible bonds in this location delete
l-=1 #go to next possible spot where the last molecule atom could be bonded
#section which assigns a bond to the last molecule atom which can be bonded.
self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,l))
bonds+=1
if bonds==bondnumber: return
del possiblebonds[l] #deletes possible bonds to the last molecule which can be bonded
l=len(self.bondinginformation)-1
possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1])
if bondnumber-bonds>=len(possiblebonds): #the rest of the bonds are assigned in order
while possiblebonds!=[]:
if possiblebonds[0]==[]: #if row of possible bonds is empty then delete the row
del possiblebonds[0]
else: #else find a bond in the row of possible bonds
self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,0))
#dont need to update bonds since possiblebonds will become empty at the same time or before bonds
#is equal to bondnumber
del possiblebonds[0]
l=len(self.bondinginformation)-1
possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1])
return
else: #the rest of the bonds are assigned randomly
from random import randint #use to randomly choose an index to bond.
while bonds<bondnumber:
if possiblebonds==[]:break #exits while loop when possiblebonds has become an empty set.
#this ensures that in the case there are not egnough viable bonds from possiblebonds to
#reach the bondnumber. Than, the while loop will not become infinite.
l=len(possiblebonds)-1
i=randint(0,l)
if possiblebonds[i]==[]: #if row of possible bonds is empty then delete the row
del possiblebonds[i]
else: #else find a bond in the row of possible bonds
self.bondinginformation.append(self.particlebondinginfo(particle,possiblebonds,i))
bonds+=1
del possiblebonds[i]
l=len(self.bondinginformation)-1
possiblebonds=self.updatepossiblebonds(possiblebonds,self.bondinginformation[l][1])
return
def particlebondinginfo(self,particle,possiblebonds,index1):
"""Takes list of possiblebonds and assigns a bond from possiblebonds[index1].
Then assigns false to particle.setbool at the bonding location
to ensure no more bonds can form with that particle.
Returns the assigned bond."""
from random import randint #use to randomly choose 0 and 1 where 1 is keep.
for i in range(len(possiblebonds[index1])):
if i==len(possiblebonds[index1])-1: #automattically bonds under this condition
break
else: #bond has random chance of forming
if randint(0,1)==1: #bond forms
break
else: #no bond forms
continue
particle.setbool(possiblebonds[index1][i][1],False)#makes sure no more bonds can form
return possiblebonds[index1][i]
def updatepossiblebonds(self,possiblebonds,particlenum):
"""finds the particle number in the remaining possible bonds than deletes that bonding information.
This algorithm assumes that particlenum can exist only once in each row of possiblebonds."""
for i in range(len(possiblebonds)):
for j in range(len(possiblebonds[i])):
if possiblebonds[i][j][1]==particlenum:
del possiblebonds[i][j]
break #go to next row of possiblebonds
return possiblebonds
def bondtoparticle(self,particle,atomid,newid,newcharge):
"""Bonds molecule to particle. Particle is a particlesurface object.
Moves the molecule's atoms bonded to the particle to the particle's position.
Alters the atomtype of the molecule's atoms bonded to the particle to newid.
Alters the charge of the molecule's atoms bonded to the particle to newcharge
Because bonds have formed between the molecule's atoms and the particle's atoms,
atoms with atomid on the molecule need to be removed otherwise the molecule's atoms will have to many bonds.
self.deleteatoms will take care of deleting the atoms atomid and
atoms down the polymer chain in the direction away from the molecule/particle bond.
Now remove the bonded particle atoms from the particle surface becaused the bonded molecule atoms
have replaced those particle atoms and their bonds with the rest of the particle.
Newid must be a string. Atomid can now be a list or an integer.
The list is a list of atom's id values and the single integer is atom's atomtype"""
#moves the molecule's atoms bonded to the particle to the particle's position
#need to do this step first before the atom id's and row indexes get out of sync
#which will occur after atoms get deleted.
print 'modifying the bonded molecules atom information'
for i in range(len(self.bondinginformation)):
#patom=particle.getsurfatom(self.bondinginformation[i][1])
#if self.atomtype=='molecular':#3-5->x,y,z[molecular]
# for j in range(3,6):
# self.modifyatom(self.bondinginformation[i][0]+1,j,patom[j])#uses the atomid here
#elif self.atomtype=='full':#4-6->x,y,z[full]
# for j in range(4,7):
# self.modifyatom(self.bondinginformation[i][0]+1,j,patom[j])#uses the atomid here
#alters the atomtype to newid of the molecule's atoms bonded to the particle.
self.modifyatom(self.bondinginformation[i][0]+1,2,newid) #uses the atomid here
if self.atomtype=='full':
# Alters the charge of the molecule's atoms bonded to the particle to newcharge
self.modifyatom(self.bondinginformation[i][0]+1,3,newcharge)
#This is a seperate loop so atom id's and row indexes for previous steps wont get out of sync
#create atomnumbers to begin deletion process.
atomnumbers=[]
for i in range(len(self.bondinginformation)):
atomnumbers.append(self.bondinginformation[i][0]+1)#using actual atomnumbers rather than the row index
#Call algorithm to find all atoms bonded to atomnumbers in the direction of atomid.
#Than delete those atoms and all molecule structures which contain those atoms.
self.deleteatoms(atomnumbers,atomid)
print 'begining deletion process of the surface atoms for which the molecule atoms have replaced'
#Goes through the bondinginformation and superficially removes the surfaceatom
for i in range(len(self.bondinginformation)):
particle.removesurfatom(self.bondinginformation[i][1]) #uses the row number
#Allows the particle extract method used outside this class to remove this atom.
def createxyz(self,file,data,routine='mass', values=None):
"""Two possible routines one to use the masses from data and the other to use the atom type and values supplied by the user.
The mass version is assessed by setting the routine to 'mass' which is the default method.
The other version is assessed by setting the routine to 'atomtype'.
The other version takes values which is a list containing the value the user wants to assign those atomtypes to.
The atomtypes of the values in the list will start at 1 even if no atoms in molecule use 1.
This makes it easier to find the atomtype and assign the value from the list
All atom data is assumed to have image flags in the data."""
f=open(file,'w')
f.write('{0}\n'.format(len(self.atoms)))
f.write('atoms\n')
if routine=='mass':
for line in self.atoms:
type=int(line[2])-1
mass=data.masses[type][1]
if mass=='12.0107': elementnum=6
elif mass=='15.9994': elementnum=8
elif mass=='26.9815':elementnum=13
else:
print 'no matching mass value. The method will exit'
return
l=len(line)-1
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3]))
f.close()
elif routine=='atomtype':
for line in self.atoms:
type=int(line[2])-1
elementnum=values[type]
l=len(line)-1
f.write('{0} {1} {2} {3}\n'.format(elementnum,line[l-5],line[l-4],line[l-3]))
f.close()
c4science · Help