import numpy as np
#Dump to lammps positions for file
#Note only works for classic lammps atom dumps

#reading in dump file
name = str(input("What is the name of the dump file?"))
file = open(name, "r")
if(name[len(name)-9:len(name)]!="lammpstrj"):
	print("Error: please abort program")
content = file.read()
data = content.split()

step = str(input("What step would you like to pull? (l for last)"))

#reading in reference file
name2 = str(input("What is the name of the reference system file?"))
file2 = open(name2, "r")
content2 = file2.read()
start = content2.find("Atoms")+5
stop = content2.find("Bonds")
atoms = content2[start:stop]
#print(atoms)
atomArr = atoms.split()

#determining bounds as stored in scaled format --> x stored * range of x + xmin = x real
xmin = -1
xmax = -1
ymin = -1
ymax = -1
zmin = -1
zmax = -1
foundData = False
beginIndex = -1
endIndex = -1

if(step=="l"): #if last easy to find
	endIndex = len(data)
	i=1
	while(i<len(data)):
		if(data[len(data)-i]=="zs"):
			beginIndex=len(data)-i+1
			xmax=float(data[len(data)-i-12])
			xmin=float(data[len(data)-i-11])
			ymax=float(data[len(data)-i-10])
			ymin=float(data[len(data)-i-9])
			zmax=float(data[len(data)-i-8])
			zmin=float(data[len(data)-i-7])
			i=len(data)
		i+=1
		
else: #if not last then find the timestep that matches
	k=int(step)
	i=0
	while(i<len(data)):
		if(data[i]=="TIMESTEP"):
			if(int(data[i+1])==k):
				xmax=float(data[i+13])
				xmin=float(data[i+14])
				ymax=float(data[i+15])
				ymin=float(data[i+16])
				zmax=float(data[i+17])
				zmin=float(data[i+18])
				i+=25
				beginIndex = i
		if(beginIndex != -1 and data[i]=="ITEM:"):
			endIndex = i
			i=len(data)
		i+=1
	
#determine ranges	
xScale=xmax-xmin
yScale=ymax-ymin
zScale=zmax-zmin
#print(xScale)
#print(yScale)
#print(zScale)

newPos = data[beginIndex:endIndex]
#print(newPos)
#print(beginIndex)
#print(endIndex)
#print(len(newPos))
#print(len(atomArr))

#read in, convert, and store positions
#note positions are NOT stored in order
i=0
while(i<int(len(newPos)/5)):
	ind = int(newPos[i*5]) # gives some atom to update
	atomArr[(ind-1)*7+4]=float(newPos[i*5+2])*xScale + xmin
	atomArr[(ind-1)*7+5]=float(newPos[i*5+3])*yScale + ymin
	atomArr[(ind-1)*7+6]=float(newPos[i*5+4])*zScale + zmin
	#print(float(newPos[i*5+2])*xScale + xmin)
	i+=1
	
#Export copiable text to replace atoms section
str1=""

i=0
while(i<int(len(atomArr)/7)):
	str1 += "     " + str(atomArr[i*7]) + "      " + str(atomArr[i*7+1])+ "      " + str(atomArr[i*7+2]) + " " + str(atomArr[i*7+3]) + "     " + str(atomArr[i*7+4]) + "  " + str(atomArr[i*7+5]) + "  " + str(atomArr[i*7+6]) +  "\n"
	#print(i*7)
	i+=1

print(str1)