#!@PYTHON@
#
# Create the steps from bcc to fcc
#
# Note: PCP config file should contain two configurations which are linearly interpolated in the given amount of steps

import sys, string, re, math, os
wrerr = sys.stderr.write

if len(sys.argv) < 3:
  wrerr("Usage: "+sys.argv[0]+" <config> <steps>\n")
  sys.exit(1)

config = open(sys.argv[1], "r")
steps = int(sys.argv[2])

# find the lines where source and target configuration start
nr=0
source=0
target=0
box=0
boxlength = []
for line in config:
    nr=nr+1
    if re.compile("BoxLength").match(line):
        box=nr
    if box != 0 and nr > box and nr < box+4:
        numbers = line.split()
        for j in range(len(numbers)):
            boxlength.append(float(numbers[j]))
    if re.compile("Ion_Type[0-9]+[_]+").match(line) and source == 0:
            source=nr
    if "# ====== MD step" in line:
        target=nr+1
if target == 0:
    wrerr("Could find only a single geometry, need at least two!\n")
    sys.exit(1)
print "Changing from config starting at line %d to the one starting at line %d ..." % (source, target)

#for i in range(len(boxlength)):
 #   print "%lg\t" % (boxlength[i])
# parse the source configuration
config.seek(0)
n=0
nr=0
sourcepos = []
sourcemax = [0., 0., 0.]
for line in config:
    nr=nr+1
    if source > nr:
        continue
    else:
        if re.compile("Ion_Type[0-9]+[_]+").match(line):
            words = line.split()
            print words
            sourcepos.append(words[0]+"="+words[1]+"="+words[2]+"="+words[3])
            for j in range(3):
                 if sourcemax[j] < float(words[j+1]):
                     sourcemax[j] = float(words[j+1])  
            n=n+1
        else:
            if "#" in line:	# we break on the next comment line after source
                break
print "I found ",n," atoms in source config."

# parse the target configuration
config.seek(0)
n=0
nr=0
targetpos = []
targetmax= [0., 0., 0.]
for line in config:
    nr=nr+1
    if target > nr:
        continue
    else:
        if re.compile("Ion_Type[0-9]+[_]+").match(line):
            words = line.split()
            print words
            targetpos.append(words[0]+"="+words[1]+"="+words[2]+"="+words[3])
            for j in range(3):
                 if targetmax[j] < float(words[j+1]):
                     targetmax[j] = float(words[j+1])  
            n=n+1
        else:
            if "#" in line:	# we break on the next comment line after source
                break
print "I found ",n," atoms in target config."
targetmax[2] = 2.64357
# Now as we have the two configs, we need go through the steps desired and write the output files
diff = [0., 0., 0.]
BoxLength = [0., 0., 0., 0., 0., 0.]
for i in range(steps+1):
    # change the boxlength
    if i == 0:
        for j in range(6):
            BoxLength[j] = boxlength[j]
    else:
        BoxLength[0] = boxlength[0] + boxlength[0]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
        BoxLength[1] = boxlength[1] + boxlength[1]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
        BoxLength[2] = boxlength[2] + boxlength[2]*(targetmax[1]/sourcemax[1]-1.)*(float(i)/float(steps))
        BoxLength[3] = boxlength[3] + boxlength[3]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
        BoxLength[4] = boxlength[4] + boxlength[4]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
        BoxLength[5] = boxlength[5] + boxlength[5]*(targetmax[2]/sourcemax[2]-1.)*(float(i)/float(steps))
    # write the initial config file to the new destination
    newdir="step%02d" % (i)
    if not os.path.isdir(newdir):
        if os.path.isfile(newdir):
            raise OSError("a file with the same name as the desired " \
                      "dir, '%s', already exists." % newdir)    
        else:
            os.mkdir(newdir)
    output = file("step%02d/%s" % (i,sys.argv[1]),"w")
    config.seek(0)
    nr=0
    for line in config:
        nr=nr+1
        if box < nr and box+4 > nr:
            j=nr-box-1
            for m in range(j+1):
                if j==2:
                    j=j+1
                    output.write("%lg\t" % (BoxLength[j+m]))
                else:
                    output.write("%lg\t" % (BoxLength[j+m]))
            output.write("\n")
            continue
        output.write(line)
        if source <= nr+1:
            break
    # add the shifted coordinates
    for j in range(n):
        sourceatoms = sourcepos[j].split("=")
        targetatoms = targetpos[j].split("=")
        for m in range(3):
            diff[m] = (float(targetatoms[m+1])-float(sourceatoms[m+1]))*float(i)/float(steps)
        output.write("%s\t%2.10f\t%2.10f\t%2.10f\t#\n" % (sourceatoms[0], float(sourceatoms[1])+diff[0], float(sourceatoms[2])+diff[1], float(sourceatoms[3])+diff[2]))
    output.close()
config.close()