#
#
#	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()