import numpy as np
import math

#types1 = np.array( (12.011,12.011,15.999, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008))
#types2 = np.array((12.011, 12.011, 15.999, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 12.011, 12.011, 12.011, 15.999, 15.999, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008, 1.008))

file = open("PLA3.lmp", "r")
#This converter is intended for raw LigParGen files to be utilized to create
#further LAMMPS data input files
#Only error surfaces in the pair coeffs section
# LigParGen gives ID value1 value2
# LAMMPS expects  ID ID value1 value2 to define LJ well

# Read the entire content of the file
content = file.read()

#Maybe Unnecessary Tools
#Masses and types per entry (maybe unnecessary)
start = content.find("Masses")+6
stop = content.find("Pair Coeffs")
masses = content[start:stop]
intermediate = masses.split()
types1 = np.array(())

i=1
while(i<len(intermediate)):
	types1= np.append(types1,intermediate[i])

	i+=2

conv = {'12.011':"C",'15.999':"O",'1.008':"H"}
#conv2 = {12.011:"C",15.999:"O",1.008:"H"}


i=0
str1 = ""
while(i<len(types1)):
	str1 += conv[types1[i]]
	i+=1

#print("Atoms: ",str1)

#Pair Coeff + finding duplicates (maybe unnecessary)
start = content.find("Pair Coeffs")+11
stop = content.find("Bond Coeffs")
pairs = content[start:stop]
intermediate = pairs.split()
types2 = np.array([[None]*2]*int(len(intermediate)/3)) #4 if LAMMPS Ready
lrPair = ""

#print(intermediate)
i=1 #2 if LAMMPS ready
while(i<len(intermediate)):
	types2[int((i-1)/3)] = [intermediate[i],intermediate[i+1]] #4 if LAMMPS ready
	lrPair += "     " + str(int((i-1)/3+1)) + "  "+ str(int((i-1)/3+1))+ "   "+ str(intermediate[i])+"   "+str(intermediate[i+1]) +"\n"
	i+=3 #4 if LAMMPS ready
	
#print(types2)

pairTypes = {(-1, -1): 'Null'}
pairList = np.array(())

i=0
while(i<len(types2)):
	try: #if pair type has already been defined
		pType = pairTypes[(types2[i][0],types2[i][1])]
		pairList = np.append(pairList,pType)
	except: #if it hasn't, add it
		pairTypes[(types2[i][0],types2[i][1])] = len(pairTypes)
		pType = pairTypes[(types2[i][0],types2[i][1])] #should be len(pairTypes)
		pairList = np.append(pairList,pType)
	i+=1


#print(pairTypes)
#pairTypes = {('0.066', '3.5000000'): 1, ('0.170', '3.1200000'): 2, ('0.070', '3.5500000'): 3, ('0.210', '2.9600000'): 4, ('0.140', '2.9000000'): 5, ('0.030', '2.5000000'): 6, ('0.000', '0.0000000'): 7}
#print(pairList)
#print(bool(len(pairList)==len(types2)))
#print("Atoms: ",str1)

#All 30's in this section are referring to the final hydrogen and only hold for file PLA3.lmp


#Bond Coeff + finding duplicates (maybe unnecessary)
start = content.find("Bond Coeffs")+11
stop = content.find("Angle Coeffs")
bondC = content[start:stop]
intermediate = bondC.split()
types3 = np.array([[None]*2]*int(len(intermediate)/3)) 

#print(intermediate)
i=1 
while(i<len(intermediate)):
	types3[int((i-1)/3)] = [intermediate[i],intermediate[i+1]]
	i+=3 
	
#print(types2)

bondTypes = {(-1,-1):'Null'}
bondList = np.array(())

i=0
while(i<len(types3)):
	try: #if pair type has already been defined
		bType = bondTypes[(types3[i][0],types3[i][1])]
		bondList = np.append(bondList,bType)
	except: #if it hasn't, add it
		bondTypes[(types3[i][0],types3[i][1])] = len(bondTypes)
		bType = bondTypes[(types3[i][0],types3[i][1])] #should be len(pairTypes)
		bondList = np.append(bondList,bType)
	i+=1



#Analyzing file to determine necessary additions for extended monomers
#Bonds + finding relevant bonds to each monomer
start = content.find("Bonds")+5
stop = content.find("Angles")
bond1 = content[start:stop]
intermediate = bond1.split()
bonds = np.array([[None]*3]*int(len(intermediate)/4))

#addedAtoms = np.array([7,8,9,10,11,17,18,19,20]) #reference numbers for atom types to be added for PLA2.lmp
addedAtoms = np.array([7,8,9,10,11,22,23,24,25]) #reference numbers for atom types to be added for PLA3.lmp (for safety)

i=0
while(i<len(bonds)):
	bonds[int(i)][0] = intermediate[4*i+1]
	bonds[int(i)][1] = intermediate[4*i+2]
	bonds[int(i)][2] = intermediate[4*i+3]
	i+=1
	#print(i)

bondIndex = i #to be used later
bondsRel = np.array(())
bondsFin = np.array(()) #finding bonds and other relationships with final hydrogen
bondsSt = np.array(())

i=0
while(i<len(bonds)):
	j=0
	#print(bonds[i])
	while(j<len(addedAtoms)):
		if((int(bonds[i][1])==addedAtoms[j] or int(bonds[i][2])==addedAtoms[j]) and not(int(bonds[i][1])==12 or int(bonds[i][2])==12)):
			#exclude 12 as each monomer should not link to the next before it has been created
			bondsRel = np.append(bondsRel, i)
			bondsSt = np.append(bondsSt, i)
			j=len(addedAtoms)
			#print("here")
		j+=1
	if(int(bonds[i][1])==30 or int(bonds[i][2])==30):
		bondsFin = np.append(bondsFin, i)
	i+=1
bondIndex -= len(bondsRel)
#print(bonds)
#print(bondsRel)
#print(bondsFin)
#print(bondsSt)

#Angles + finding relevant angles to each monomer
start = content.find("Angles")+6
stop = content.find("Dihedrals")
angles1 = content[start:stop]
intermediate = angles1.split()
angles = np.array([[None]*4]*int(len(intermediate)/5))

#addedAtoms = np.array([7,8,9,10,11,17,18,19,20]) #reference numbers for atom types to be added for PLA2.lmp
addedAtoms = np.array([7,8,9,10,11,22,23,24,25]) #reference numbers for atom types to be added for PLA3.lmp (for safety)

i=0
while(i<len(angles)):
	angles[int(i)][0] = intermediate[5*i+1]
	angles[int(i)][1] = intermediate[5*i+2]
	angles[int(i)][2] = intermediate[5*i+3]
	angles[int(i)][3] = intermediate[5*i+4]
	i+=1
	#print(i)

anglesIndex = i #to be used later
anglesRel = np.array(())
anglesFin = np.array(())
anglesSt = np.array(())

i=0
while(i<len(angles)):
	j=0
	#print(bonds[i])
	while(j<len(addedAtoms)):
		if(((int(angles[i][1])==addedAtoms[j] or int(angles[i][2])==addedAtoms[j]) or int(angles[i][3])==addedAtoms[j]) and not((int(angles[i][1])==12 or int(angles[i][2])==12) or int(angles[i][3])==12)):
			#exclude 12 as each monomer should not link to the next before it has been created
			anglesRel = np.append(anglesRel, i)
			j=len(addedAtoms)
			#print("here")
		j+=1
	if((int(angles[i][1])==30 or int(angles[i][2])==30) or int(angles[i][3])==30):
		anglesFin = np.append(anglesFin, i)
	if((int(angles[i][1])<26 and int(angles[i][2])<26) and int(angles[i][3])<26):
		anglesSt = np.append(anglesSt, int(i))
	i+=1
anglesIndex -= len(anglesRel)

#print(angles)
#print(anglesRel)
#print(anglesFin)


#Repeating for Dihedrals
start = content.find("Dihedrals")+9
stop = content.find("Impropers")
dihe1 = content[start:stop]
intermediate = dihe1.split()
dihe = np.array([[None]*5]*int(len(intermediate)/6))

#addedAtoms = np.array([7,8,9,10,11,17,18,19,20]) #reference numbers for atom types to be added for PLA2.lmp
addedAtoms = np.array([7,8,9,10,11,22,23,24,25]) #reference numbers for atom types to be added for PLA3.lmp (for safety)

i=0
while(i<len(dihe)):
	dihe[int(i)][0] = intermediate[6*i+1]
	dihe[int(i)][1] = intermediate[6*i+2]
	dihe[int(i)][2] = intermediate[6*i+3]
	dihe[int(i)][3] = intermediate[6*i+4]
	dihe[int(i)][4] = intermediate[6*i+5]
	i+=1
	#print(i)

diheIndex = i #to be used later
diheRel = np.array(())
diheFin = np.array(())
diheSt = np.array(())

i=0
while(i<len(dihe)):
	j=0
	#print(bonds[i])
	while(j<len(addedAtoms)):
		if(((int(dihe[i][1])==addedAtoms[j] or int(dihe[i][2])==addedAtoms[j]) or (int(dihe[i][3])==addedAtoms[j] or int(dihe[i][4])==addedAtoms[j])) and not((int(dihe[i][1])==12 or int(dihe[i][2])==12) or (int(dihe[i][3])==12 or int(dihe[i][4])==12))):
			#exclude 12 as each monomer should not link to the next before it has been created also kinda like fencepost problem
			diheRel = np.append(diheRel, i)
			j=len(addedAtoms)
			#print("here")
		j+=1
	if((int(dihe[i][1])==30 or int(dihe[i][2])==30) or (int(dihe[i][3])==30 or int(dihe[i][4])==30)):
		diheFin = np.append(diheFin, i)
	if((int(dihe[i][1])<26 and int(dihe[i][2])<26) and (int(dihe[i][3])<26 and int(dihe[i][4])<26)):
		diheSt = np.append(diheSt, int(i))
	i+=1
diheIndex -= len(diheRel)

#print(dihe)
#print(diheRel)
#print(diheFin)


#Repeating for Impropers
start = content.find("Impropers")+9
imp1 = content[start:len(content)-1]
intermediate = imp1.split()
imp = np.array([[None]*5]*int(len(intermediate)/6))

#addedAtoms = np.array([7,8,9,10,11,17,18,19,20]) #reference numbers for atom types to be added for PLA2.lmp
addedAtoms = np.array([7,8,9,10,11,22,23,24,25]) #reference numbers for atom types to be added for PLA3.lmp (for safety)

i=0
while(i<len(imp)):
	imp[int(i)][0] = intermediate[6*i+1]
	imp[int(i)][1] = intermediate[6*i+2]
	imp[int(i)][2] = intermediate[6*i+3]
	imp[int(i)][3] = intermediate[6*i+4]
	imp[int(i)][4] = intermediate[6*i+5]
	i+=1
	#print(i)

impIndex = i #to be used later
impRel = np.array(())
impFin = np.array(())
impSt = np.array(())

i=0
while(i<len(imp)):
	j=0
	#print(bonds[i])
	while(j<len(addedAtoms)):
		if(((int(imp[i][1])==addedAtoms[j] or int(imp[i][2])==addedAtoms[j]) or (int(imp[i][3])==addedAtoms[j] or int(imp[i][4])==addedAtoms[j])) and not((int(imp[i][1])==12 or int(imp[i][2])==12) or (int(imp[i][3])==12 or int(imp[i][4])==12))):
			#exclude 12 as each monomer should not link to the next before it has been created also kinda like fencepost problem
			impRel = np.append(impRel, i)
			j=len(addedAtoms)
			#print("here")
		j+=1
	if((int(imp[i][1])==30 or int(imp[i][2])==30) or (int(imp[i][3])==30 or int(imp[i][4])==30)):
		impFin = np.append(impFin, i)
	if((int(imp[i][1])<26 and int(imp[i][2])<26) and (int(imp[i][3])<26 and int(imp[i][4])<26)):
		impSt = np.append(impSt, int(i))
	i+=1
impIndex -= len(impRel)
#print(imp)
#print(impRel)
#print(impFin)

#Atoms in System Definition
start = content.find("Atoms")+5
stop = content.find("Bonds")
atoms = content[start:stop]
intermediate = atoms.split()
finalIndex =int(len(intermediate)/7) #of the .lmp file which is size (since +1)
atomsInfo = np.array([[None]*6]*finalIndex)


i=2
#print(intermediate)
while(i<len(intermediate)):
	#print(i)
	atomsInfo[int((i-2)/7)][0]=str1[(int(intermediate[i])-1)] #atom type from above
	atomsInfo[int((i-2)/7)][1]=intermediate[i] #number form type from above
	atomsInfo[int((i-2)/7)][2]=intermediate[i+1] #charge
	atomsInfo[int((i-2)/7)][3]=intermediate[i+2] #x
	atomsInfo[int((i-2)/7)][4]=intermediate[i+3] #y
	atomsInfo[int((i-2)/7)][5]=intermediate[i+4] #z
	i+=7


#Properties
PLA6Coords = np.array([[1.000, 1.00000, 0.00000],[-0.520, 1.00000, 0.00000],[-0.961, 1.00000,1.35602],[-1.069, 2.18618, -0.80654],[-1.758, 3.08034, -0.32996],[-0.681, 2.09957, -2.10694],[-1.181, 3.12465, -2.98529],[-0.332, 4.38381, -2.82296],[-1.128, 2.64282, -4.44133],[-1.025, 3.38413, -5.41166],[-1.210, 1.28250, -4.51827],[-1.171, 0.72548, -5.84361],[0.287, 0.64455, -6.29951],[-1.806, -0.67431, -5.87685],[-1.496, -1.55031, -6.67656],[-2.797, -0.80953, -4.94792],[-3.454, -2.09600, -4.90159],[-2.582, -3.07106, -4.11687],[-4.823, -1.95941, -4.21498],[-5.287, -0.91590, -3.77299],[-5.468, -3.16478, -4.15971],[-6.769, -3.16991, -3.53343],[-7.816, -2.70719, -4.54088],[-7.115, -4.59023, -3.05758],[-6.395, -5.57033, -3.18978],[-8.350, -4.61781, -2.45651],[-8.775, -5.91439, -1.97147],[-8.164, -6.19189, -0.60493],[-10.300, -5.93898, -1.84061],[-10.946, -6.90055, -1.45563],[-10.889, -4.76936, -2.19682],[1.382, 0.17265, 0.60663],[1.387, 1.92577, 0.44035],[1.404, 0.90725, -1.01250],[-0.903, 0.08105, -0.45570],[-0.813, 1.89827, 1.70574],[-2.236, 3.33081, -2.76680],[0.720, 4.17496, -3.04952],[-0.365, 4.75570, -1.79393],[-0.673, 5.18509, -3.48664],[-1.757, 1.34345, -6.53537],[0.872, 0.01579, -5.61887],[0.760, 1.63207, -6.29973],[0.370, 0.22610, -7.30768],[-3.643, -2.44954, -5.92298],[-2.416, -2.70152, -3.09876],[-1.596, -3.18322, -4.57884],[-3.043, -4.06113, -4.04940],[-6.751, -2.53137, -2.64187],[-7.836, -3.37455, -5.40988],[-7.594, -1.70346, -4.91436],[-8.819, -2.69389, -4.10244],[-8.515, -6.68762, -2.70421],[-8.442, -5.41082, 0.11165],[-7.071, -6.20479, -0.65429],[-8.504, -7.15339, -0.20524],[-10.214, -4.11904, -2.47532]]) #all coords given from PLA6
xOffset = 0
yOffset = 0
zOffset = 0
forOvito = False

#str2 = "\n When copying this into a file be sure to remove the last hydrogen and its bond, angles, dihedrals, impropers!\n\n"
i=0
k = int(input("How many monomers? (n>2)")) -2#(number of additional monomers added to 2)
str2 = ""

#21 atoms in PLA2
while(i<21-1): #minus 1 prevents last hydrogen from getting added
	if(forOvito):
		if(i==10 and k>0):
			str2 += "    " + str(i+1) + "      " + str(6) + " " #sets oxygen type to non-hydroxy for polymerization
		else:
			if(i<11):
				str2 += "    " + str(i+1) + "      " + str(i+1) + " " #gives ID, mol #, and atom type
			else: #12 should be a hydrogen bound to carbon +5 to 17,18,19
				str2 += "    " + str(i+1) + "      " + str(i+1+5) + " " #gives ID, mol #, and atom type
	
		if(i < 11):
			str2 += str(PLA6Coords[i][0]+xOffset) + " " + str(PLA6Coords[i][1]+yOffset) + " " + str(PLA6Coords[i][2]+zOffset) +"  "
		else:
			str2 += str(PLA6Coords[i+20][0]+xOffset) + " " + str(PLA6Coords[i+20][1]+yOffset) + " " + str(PLA6Coords[i+20][2]+zOffset) +"  "
		
		if(i==10 and k>0):
			str2 += str(atomsInfo[6][2]) + "    1 " #charge
		else:
			if(i<10):
				str2 += str(atomsInfo[i][2]) + "    1 " #charge
			else:
				str2 += str(atomsInfo[i+5][2]) + "    1 " #charge plus 5 to account for misallignment PLA2 to PLA3
	else:
		if(i==10 and k>0):
			str2 += "    " + str(i+1) + "      " + "1" + "      " + str(6) + " " #sets oxygen type to non-hydroxy for polymerization
			str2 += str(atomsInfo[6][2]) + "    " #charge
		else:
			if(i<10):
				str2 += "    " + str(i+1) + "      " + "1" + "      " + str(i+1) + " " #gives ID, mol #, and atom type
				str2 += str(atomsInfo[i][2]) + "    " #charge
			else: 
				str2 += "    " + str(i+1) + "      " + "1" + "      " + str(i+1+5) + " " #gives ID, mol #, and atom type
				str2 += str(atomsInfo[i+5][2]) + "    " #charge
	
		if(i < 11):
			str2 += str(PLA6Coords[i][0]+xOffset) + " " + str(PLA6Coords[i][1]+yOffset) + " " + str(PLA6Coords[i][2]+zOffset)
		else:
			str2 += str(PLA6Coords[i+20][0]+xOffset) + " " + str(PLA6Coords[i+20][1]+yOffset) + " " + str(PLA6Coords[i+20][2]+zOffset)
	str2 += "\n"
	i+=1
#print(str2)
addedAtoms = np.array([7,8,9,10,11,17,18,19,20]) #reference numbers for atom types to be added for PLA2.lmp
#bonds and other interactions detailed in source file

#longer polymers
changeX = -5.478 #the change cross two monomers in linear section (measured from ovito)
changeY = -3.73577
changeZ = 2.70408

j=0
while(j<k):
	
	#atoms addition (for use on PLA3 only) used to create monomers to replace last too.
	finalIndex = 21 # so as to make it remake the last polymer
	i=0
	while(i<9):
		addOn="    "
		if(forOvito):
			if(k-1==j and i==4): #if last monomer and oxygen
				addOn += str(finalIndex+i+j*9) + "      "
				addOn += str(addedAtoms[i]) +" "#atom type
			elif(i==4):	#if joining oxygen
				addOn += str(finalIndex+i+j*9) + "      " 
				addOn += str(6) +" "#atom type
			else:
				addOn += str(finalIndex+i+j*9) + "      " 
				addOn += str(addedAtoms[i]) +" "#atom type
			if(j<4): # if in well defined location domain
				if(j==3 and (i==4 or i==3)): #if 4th monomer and oxygen
					#accesses oxygen from two monomers back and adds standard change for each
					addOn += str(PLA6Coords[5*j+1+i][0]+xOffset+changeX) + " " + str(PLA6Coords[5*j+1+i][1]+yOffset+changeY) + " " + str(PLA6Coords[5*j+1+i][2]+zOffset+changeZ) + "  "
				elif(i<5): #atoms 7-11 mapped to locations of 12-16 then 17-21 etc.
					addOn += str(PLA6Coords[5*j+11+i][0]+xOffset) + " " + str(PLA6Coords[5*j+11+i][1]+yOffset) + " " + str(PLA6Coords[5*j+11+i][2]+zOffset) + "  "
					#print(i, "  ", addedAtoms[i], "  ", finalIndex+i+j*9, " ", PLA6Coords[5*j+11+i][0] +xOffset, " ", PLA6Coords[5*j+11+i][1]+yOffset, " ",PLA6Coords[5*j+11+i][2]+zOffset)
				else: #atoms 17-20 mapped to locations of 37-40 then 41-44 etc.
					addOn += str(PLA6Coords[4*j+41+i-5-1][0]+xOffset) + " " + str(PLA6Coords[4*j+41+i-5-1][1]+yOffset) + " " + str(PLA6Coords[4*j+41+i-5-1][2]+zOffset) + "  "
			else:
				if(j%2==0): #if even addition
					if(i<5):#accesses from two monomers back and adds standard change for each
						addOn += str(PLA6Coords[5*4+1+i][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+i][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+i][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
					else: #atoms 17-20 mapped to locations of 37-40 then 41-44 etc.
						addOn += str(PLA6Coords[4*4+33+i-5-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*4+33+i-5-1][1]+yOffset+ changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*4+33+i-5-1][2]+zOffset+ changeZ*(math.floor(j/2)-1)) + "  "
				else: #if even addition
					if(i<5):#accesses from two monomers back and adds standard change for each
						if(i==3): #switch oxygen 10 and 11 since they are flipped for reference last
							addOn += str(PLA6Coords[5*5+1+4][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+4][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+4][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
						elif(i==4):
							addOn += str(PLA6Coords[5*5+1+3][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+3][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+3][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
						else:
							addOn += str(PLA6Coords[5*5+1+i][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+i][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+i][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
					else: #atoms 17-20 mapped to locations of 37-40 then 41-44 etc.
						addOn += str(PLA6Coords[4*5+33+i-5-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*5+33+i-5-1][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*5+33+i-5-1][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
			if(k-1==j and i==4): #if last monomer and oxygen
				addOn += str(atomsInfo[addedAtoms[i]][2]) + "    1 " #charge
			elif(i==4):	#if joining oxygen
				addOn += str(atomsInfo[6][2]) + "     1 " #charge
			else:
				addOn += str(atomsInfo[addedAtoms[i]][2]) + "     1 " #charge
		
		else:
			if(k-1==j and i==4): #if last monomer and oxygen
				addOn += str(finalIndex+i+j*9) + "      " + "1" + "      " #1 is for molecule number
				addOn += str(addedAtoms[i]) +" "#atom type
				addOn += str(atomsInfo[addedAtoms[i]][2]) + "    " #charge
			elif(i==4):	#if joining oxygen
				addOn += str(finalIndex+i+j*9) + "      " + "1" + "      " #1 is for molecule number
				addOn += str(6) +" "#atom type
				addOn += str(atomsInfo[6][2]) + "    " #charge
			else:
				addOn += str(finalIndex+i+j*9) + "      " + "1" + "      " #1 is for molecule number
				addOn += str(addedAtoms[i]) +" "#atom type
				addOn += str(atomsInfo[addedAtoms[i]][2]) + "    " #charge
			if(j<4): # if in well defined location domain
				if(j==3 and (i==4 or i==3)): #if 4th monomer and oxygen 10 and 11
					#accesses oxygen from two monomers back and adds standard change for each
					addOn += str(PLA6Coords[5*j+1+i][0]+xOffset+changeX) + " " + str(PLA6Coords[5*j+1+i][1]+yOffset+changeY) + " " + str(PLA6Coords[5*j+1+i][2]+zOffset+changeZ) + "  "
				elif(i<5): #atoms 7-11 mapped to locations of 12-16 then 17-21 etc.
					addOn += str(PLA6Coords[5*j+11+i][0]+xOffset) + " " + str(PLA6Coords[5*j+11+i][1]+yOffset) + " " + str(PLA6Coords[5*j+11+i][2]+zOffset)
				else: #atoms 17-20 mapped to locations of 37-40 then 41-44 etc.
					addOn += str(PLA6Coords[4*j+41+i-5-1][0]+xOffset) + " " + str(PLA6Coords[4*j+41+i-5-1][1]+yOffset) + " " + str(PLA6Coords[4*j+41+i-5-1][2]+zOffset)
			else:
				if(j%2==0): #if even addition
					if(i<5):#accesses from two monomers back and adds standard change for each
						if(i==3):
							addOn += str(PLA6Coords[5*4+1+4][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+4][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+4][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
						elif(i==4):
							addOn += str(PLA6Coords[5*4+1+3][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+3][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+3][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
						else:
							addOn += str(PLA6Coords[5*4+1+i][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+i][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*4+1+i][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
					else: #atoms 17-20 mapped to locations of 37-40 then 41-44 etc.
						addOn += str(PLA6Coords[4*4+33+i-5-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*4+33+i-5-1][1]+yOffset+ changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*4+33+i-5-1][2]+zOffset+ changeZ*(math.floor(j/2)-1)) + "  "
				else: #if even addition
					if(i<5):#accesses from two monomers back and adds standard change for each
						if(i==3):
							addOn += str(PLA6Coords[5*5+1+4][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+4][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+4][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
						elif(i==4):
							addOn += str(PLA6Coords[5*5+1+3][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+3][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+3][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
						else:
							addOn += str(PLA6Coords[5*5+1+i][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+i][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[5*5+1+i][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
					else: #atoms 17-20 mapped to locations of 37-40 then 41-44 etc.
						addOn += str(PLA6Coords[4*5+33+i-5-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*5+33+i-5-1][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[4*5+33+i-5-1][2]+zOffset+changeZ*(math.floor(j/2)-1)) + "  "
						

		addOn +="\n"
		str2 += addOn
		i+=1
	j+=1

#adding the final hydrogen
if(forOvito):
	str2+= "    " + str(finalIndex+9+(k-1)*9)+ "      " + str(30) + "    " 
	if(k<4):
		str2 += str(PLA6Coords[5*(k+1)+7][0]+xOffset) + " " + str(PLA6Coords[5*(k-1)+7][1]+yOffset) + " " + str(PLA6Coords[5*(k-1)+7][2]+zOffset)
	elif(k==4):
		str2 += str(PLA6Coords[56][0]+xOffset) + " " + str(PLA6Coords[56][1]+yOffset) + " " + str(PLA6Coords[56][2]+zOffset)
	else:
		if(j%2==1): #odd addition
			str2 += str(PLA6Coords[27-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[27-1][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[27-1][2]+zOffset+changeZ*(math.floor(j/2)-1))
		else:
			str2 += str(PLA6Coords[22-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[22-1][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[22-1][2]+zOffset+changeZ*(math.floor(j/2)-1))
	str2+= " " + str(atomsInfo[29][2]) + "      " + "1" 
else:
	str2+= "    " + str(finalIndex+9+(k-1)*9)+ "      " + "1" + "      " + str(30) +" " + str(atomsInfo[29][2]) + "    " 
	if(k<4):
		str2 += str(PLA6Coords[5*(k+1)+7][0]+xOffset) + " " + str(PLA6Coords[5*(k-1)+7][1]+yOffset) + " " + str(PLA6Coords[5*(k-1)+7][2]+zOffset)
	elif(k==4):
		str2 += str(PLA6Coords[56][0]+xOffset) + " " + str(PLA6Coords[56][1]+yOffset) + " " + str(PLA6Coords[56][2]+zOffset)
	else:
		if(j%2==1): #odd addition
			str2 += str(PLA6Coords[27-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[27-1][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[27-1][2]+zOffset+changeZ*(math.floor(j/2)-1))
		else:
			str2 += str(PLA6Coords[22-1][0]+xOffset+changeX*(math.floor(j/2)-1)) + " " + str(PLA6Coords[22-1][1]+yOffset+changeY*(math.floor(j/2)-1)) + " " + str(PLA6Coords[22-1][2]+zOffset+changeZ*(math.floor(j/2)-1))
str2+="\n\n"
str2+="Bonds\n"
str2+="\n"

#print(bondIndex)
#print(anglesIndex)
#print(diheIndex)
#print(impIndex)

#known bonds for starting
str2+= "     1      1      2      1\n     2      2      3      2\n     3      3      4      2\n     4      4      5      4\n     5      5      6      4\n     6      6      7      6\n     7      7      8      7\n     8      8      9      7\n     9      9     10      9\n    10     10     11      9\n    11     16     12      1\n    12     17     13      1\n    13     18     14      1\n    14     19     15      2\n    15     20     16      3\n    16     21     17      7\n    17     22     18      8\n    18     23     19      8\n    19     24     20      8"


#for ease of converting between the general forms of bonds and other interactions as constructed 
converterA = {
	7:(finalIndex),
	8:(finalIndex+1),
	9:(finalIndex+2),
	10:(finalIndex+3),
	11:(finalIndex+4),
	22:(finalIndex+5),
	23:(finalIndex+6),
	24:(finalIndex+7),
	25:(finalIndex+8),
	6:(finalIndex+4-9), #6 is the oxygen defined by 11 but in the previous monomer]
	16:(finalIndex+4+9), 
	4:(finalIndex+2-9), #corresponds to C 9
	14:(finalIndex+2+9),
	2:(finalIndex-9), #corresponds to C 7
	12:(finalIndex+9),
	5:(finalIndex+3-9), #corresponds to O 10
	15:(finalIndex+3+9),
	30:(finalIndex+9*2)
} #only good if last addition was in fact a monomer in this format so good only j!=0

converterB = {
	7:(finalIndex),
	8:(finalIndex+1),
	9:(finalIndex+2),
	10:(finalIndex+3),
	11:(finalIndex+4),
	22:(finalIndex+5),
	23:(finalIndex+6),
	24:(finalIndex+7),
	25:(finalIndex+8),
	6:11, #constructing monomers from 2
	4:9,
	2:7,
	5:10
}


#bonds addition
str2 += "\n" #to separate atoms addition to bonds addition
#bondIndex was defined earlier to be the index used next for bonds
#bondIndex + 1 is where should start but bondIndex alone accounts for removing last hydrogen

j=0
while(j<k):
	i=0
	while(i<len(bondsRel)):
		addOn="    "
		addOn += str(bondIndex) + "   "
		addOn += str(int(bondsRel[i]+1)) + "   " #gives the bondtype number referenced in the coeff section
		if(j!=0):
			addOn += str(converterA[int(bonds[int(bondsRel[i])][1])]+9*j) + "   "
			addOn += str(converterA[int(bonds[int(bondsRel[i])][2])]+9*j)
		else:
			addOn += str(converterB[int(bonds[int(bondsRel[i])][1])]+9*j) + "   "
			addOn += str(converterB[int(bonds[int(bondsRel[i])][2])]+9*j)
		
		addOn +="\n"
		str2 += addOn

		bondIndex +=1
		i+=1
	j+=1

#Adding bond(s) for final hydrogen
str2+="    "+str(bondIndex) + "   " + str(int(bondsFin[0]+1)) + "   " 
str2 += str(converterA[int(bonds[int(bondsFin[0])][1])]+9*(k-2)) + "   "
str2 += str(converterA[int(bonds[int(bondsFin[0])][2])]+9*(k-2)) +"\n"
	

#angles addition
str2 += "\n" #to separate bonds addition to angles addition
str2 += "Angles\n"
str2 +="\n"
#known angles start
str2 += "     1      1      1      2      3\n     2      2      1      2      4\n     3      3      2      4      5\n     4      4      2      4      6\n     5      5      4      6      7\n     6      6      6      7      8\n     7      7      6      7      9\n     8      8      7      9     10\n     9      9      7      9     11\n    10     15      2      1     12\n    11     16      2      1     13\n    12     17      2      1     14\n    13     18      1      2     15\n    14     19      2      3     16\n    15     20      6      7     17\n    16     21      7      8     18\n    17     22      7      8     19\n    18     23      7      8     20\n    19     29      5      4      6\n    20     30      3      2     15\n    21     31     18      8     19\n    22     34      8      7     17\n    23     36      4      2     15\n    24     37     10      9     11\n    25     38      8      7      9\n    26     39     13      1     14\n    27     40     18      8     20\n    28     41     12      1     14\n    29     42      9      7     17\n    30     43      3      2      4\n    31     46     12      1     13\n    32     47     19      8     20\n"
#angleIndex was defined earlier to be the index used next for angles 
#anglesIndex + 1 is where should start but anglesIndex alone accounts for removing last hydrogen

#print(len(anglesRel))
#print(len(anglesFin))

j=0
while(j<k):	
	i=0
	while(i<len(anglesRel)):
		addOn="    "
		addOn += str(anglesIndex) + "   "
		addOn += str(int(anglesRel[i]+1)) + "   " #gives the bondtype number referenced in the coeff section
		if(j!=0):
			addOn += str(converterA[int(angles[int(anglesRel[i])][1])]+9*j) + "   "
			addOn += str(converterA[int(angles[int(anglesRel[i])][2])]+9*j) + "   "
			addOn += str(converterA[int(angles[int(anglesRel[i])][3])]+9*j)
		else:
			addOn += str(converterB[int(angles[int(anglesRel[i])][1])]+9*j) + "   "
			addOn += str(converterB[int(angles[int(anglesRel[i])][2])]+9*j) + "   "
			addOn += str(converterB[int(angles[int(anglesRel[i])][3])]+9*j)
	
		addOn +="\n"
		str2 += addOn

		anglesIndex +=1
		i+=1
	j+=1

#Adding angles for final hydrogen
n=0
while(n<len(anglesFin)):
	str2 += "    " + str(anglesIndex) + "   "+ str(int(anglesFin[n]+1)) + "   "
	str2 += str(converterA[int(angles[int(anglesFin[n])][1])]+9*(k-2)) + "   "
	str2 += str(converterA[int(angles[int(anglesFin[n])][2])]+9*(k-2)) + "   "
	str2 += str(converterA[int(angles[int(anglesFin[n])][3])]+9*(k-2)) + "\n"
	anglesIndex +=1
	n+=1


#dihedrals addition
str2 += "\n" #to separate angles addition to dihedrals addition
str2 += "Dihedrals\n"
str2 +="\n"
#known dihedrals starting
str2 +="     1      1      5      4      2      1\n     2      2      7      6      4      2\n     3      3      8      7      6      4\n     4      4     10      9      7      6 \n     5      8     12      1      2      3 \n     6      9     16      3      2      1 \n     7     10     18      8      7      6 \n     8     13     16      3      2      4 \n     9     15     11      9      7      6 \n    10     16     16      3      2     15 \n    11     18     18      8      7      9\n    12     20     14      1      2      4 \n    13     21      5      4      2      3 \n    14     22     17      7      9     10 \n    15     23     19      8      7     17 \n    16     25     15      2      4      5 \n    17     26     15      2      1     13 \n    18     27      9      7      6      4 \n    19     28     13      1      2      4 \n    20     29     14      1      2      3 \n    21     30     15      2      1     12 \n    22     32     17      7      6      4 \n    23     34     17      7      9     11 \n    24     36     19      8      7      9 \n    25     37     20      8      7     17 \n    26     38     15      2      1     14 \n    27     39      6      4      2      1 \n    28     44     19      8      7      6 \n    29     46      6      4      2      3 \n    30     47     11      9      7      8 \n    31     48     20      8      7      6 \n    32     49     18      8      7     17 \n    33     50     15      2      4      6 \n    34     51     20      8      7      9 \n    35     53     12      1      2      4 \n    36     54     13      1      2      3\n    37     55     10      9      7      8\n    38     58      7      6      4      5\n"
#diheIndex was defined earlier to be the index used next for dihedral 
#again diheIndex+1 where should start but diheIndex +1 -2 accounts for last hydrogen
diheIndex -=1

j=0
while(j<k):	
	i=0
	while(i<len(diheRel)):
		addOn="    "
		addOn += str(diheIndex) + "   "
		addOn += str(int(diheRel[i]+1)) + "   " #gives the bondtype number referenced in the coeff section
		if(j!=0):
			addOn += str(converterA[int(dihe[int(diheRel[i])][1])]+9*j) + "   "
			addOn += str(converterA[int(dihe[int(diheRel[i])][2])]+9*j) + "   "
			addOn += str(converterA[int(dihe[int(diheRel[i])][3])]+9*j) + "   "
			addOn += str(converterA[int(dihe[int(diheRel[i])][4])]+9*j)
		else:
			addOn += str(converterB[int(dihe[int(diheRel[i])][1])]+9*j) + "   "
			addOn += str(converterB[int(dihe[int(diheRel[i])][2])]+9*j) + "   "
			addOn += str(converterB[int(dihe[int(diheRel[i])][3])]+9*j) + "   "
			addOn += str(converterB[int(dihe[int(diheRel[i])][4])]+9*j)
			
		addOn +="\n"
		str2 += addOn

		diheIndex +=1
		i+=1
	j+=1

#Adding dihedrals for final hydrogen
n=0
while(n<len(diheFin)):
	str2 += "    " + str(diheIndex) + "   "+ str(int(diheFin[n]+1)) + "   "
	str2 += str(converterA[int(dihe[int(diheFin[n])][1])]+9*(k-2)) + "   "
	str2 += str(converterA[int(dihe[int(diheFin[n])][2])]+9*(k-2)) + "   "
	str2 += str(converterA[int(dihe[int(diheFin[n])][3])]+9*(k-2)) + "   "
	str2 += str(converterA[int(dihe[int(diheFin[n])][4])]+9*(k-2)) + "\n"
	diheIndex +=1
	n+=1


#impropers addition
str2 += "\n" #to separate dihedrals addition to impropers addition
str2 += "Impropers\n"
str2 += "\n"
#known impropers for starting
str2 += "     1      1      2      1      3      4 \n     2      2      4      2      5      6 \n     3      3      7      9      6      8 \n     4      4      9     10     11      7 \n     5      7      1     13     12      2\n     6      8      1      2     14     12\n     7      9      2      1      3     15 \n     8     10      7      8     17      6\n     9     11      8     18      7     19 \n    10     12      8     20      18     7 \n"
#impIndex was defined earlier to be the index used next for impropers
#No need to account for last hydrogen as doesn't impact impropers
impIndex +=1 

j=0
while(j<k):	
	i=0
	while(i<len(impRel)):
		addOn="    "
		addOn += str(impIndex) + "   "
		addOn += str(int(impRel[i]+1)) + "   " #gives the bondtype number referenced in the coeff section
		if(j!=0):
			addOn += str(converterA[int(imp[int(impRel[i])][1])]+9*j) + "   "
			addOn += str(converterA[int(imp[int(impRel[i])][2])]+9*j) + "   "
			addOn += str(converterA[int(imp[int(impRel[i])][3])]+9*j) + "   "
			addOn += str(converterA[int(imp[int(impRel[i])][4])]+9*j)
		else:
			addOn += str(converterB[int(imp[int(impRel[i])][1])]+9*j) + "   "
			addOn += str(converterB[int(imp[int(impRel[i])][2])]+9*j) + "   "
			addOn += str(converterB[int(imp[int(impRel[i])][3])]+9*j) + "   "
			addOn += str(converterB[int(imp[int(impRel[i])][4])]+9*j)
		
		addOn +="\n"
		str2 += addOn

		impIndex +=1
		i+=1
	j+=1

#Adding impropers for final hydrogen
n=0
while(n<len(impFin)):
	str2 += "    " + str(impIndex) + "   "+ str(int(impFin[n]+1)) + "   "
	str2 += str(converterA[int(imp[int(impFin[n])][1])]+9*(k-2)) + "   "
	str2 += str(converterA[int(imp[int(impFin[n])][2])]+9*(k-2)) + "   "
	str2 += str(converterA[int(imp[int(impFin[n])][3])]+9*(k-2)) + "   "
	str2 += str(converterA[int(imp[int(impFin[n])][4])]+9*(k-2)) + "\n"
	impIndex +=1
	n+=1		


#str2 += "  " + str(len(atomsInfo)-1) + "    1  " + str(len(atomsInfo)-1)	


#print(str2)
#print(bondsFin)

#New start section
#print("To replace openning section:")
startS = "LAMMPS data file Created by LigParGen Modified and Extended by Harrison S.\n\n"
startS +=  "      " + str(finalIndex+9+(k-1)*9) + " atoms\n"
startS += "      " + str(bondIndex) + " bonds\n"
startS += "      " + str(anglesIndex-1) + " angles\n"
startS += "      " + str(diheIndex-1) + " dihedrals\n"
startS += "      " + str(impIndex-1) + " impropers\n"
i1 = content.find("30 atom types")
i2 = content.find("Pair Coeffs")+11
i3 = content.find("Bond Coeffs")
i4 = content.find("Atoms")+5

copiedData = "      " + content[i1:i2] +"\n\n"+lrPair+"\n" +content[i3:i4]+"\n" #All simply copied data (to start of atoms section)


str2 = startS +"\n\n\n" + copiedData + "\n"+ str2

print("Start file:")
print(str2)