| 1 | #!@PYTHON@
|
|---|
| 2 | #
|
|---|
| 3 | # Create the steps from bcc to fcc
|
|---|
| 4 | #
|
|---|
| 5 | # Note: PCP config file should contain two configurations which are linearly interpolated in the given amount of steps
|
|---|
| 6 |
|
|---|
| 7 | import sys, string, re, math, os
|
|---|
| 8 | wrerr = sys.stderr.write
|
|---|
| 9 |
|
|---|
| 10 | if len(sys.argv) < 3:
|
|---|
| 11 | wrerr("Usage: "+sys.argv[0]+" <config> <steps>\n")
|
|---|
| 12 | sys.exit(1)
|
|---|
| 13 |
|
|---|
| 14 | config = open(sys.argv[1], "r")
|
|---|
| 15 | steps = int(sys.argv[2])
|
|---|
| 16 |
|
|---|
| 17 | # find the lines where source and target configuration start
|
|---|
| 18 | nr=0
|
|---|
| 19 | source=0
|
|---|
| 20 | target=0
|
|---|
| 21 | box=0
|
|---|
| 22 | boxlength = []
|
|---|
| 23 | for line in config:
|
|---|
| 24 | nr=nr+1
|
|---|
| 25 | if re.compile("BoxLength").match(line):
|
|---|
| 26 | box=nr
|
|---|
| 27 | if box != 0 and nr > box and nr < box+4:
|
|---|
| 28 | numbers = line.split()
|
|---|
| 29 | for j in range(len(numbers)):
|
|---|
| 30 | boxlength.append(float(numbers[j]))
|
|---|
| 31 | if re.compile("Ion_Type[0-9]+[_]+").match(line) and source == 0:
|
|---|
| 32 | source=nr
|
|---|
| 33 | if "# ====== MD step" in line:
|
|---|
| 34 | target=nr+1
|
|---|
| 35 | if target == 0:
|
|---|
| 36 | wrerr("Could find only a single geometry, need at least two!\n")
|
|---|
| 37 | sys.exit(1)
|
|---|
| 38 | print "Changing from config starting at line %d to the one starting at line %d ..." % (source, target)
|
|---|
| 39 |
|
|---|
| 40 | #for i in range(len(boxlength)):
|
|---|
| 41 | # print "%lg\t" % (boxlength[i])
|
|---|
| 42 | # parse the source configuration
|
|---|
| 43 | config.seek(0)
|
|---|
| 44 | n=0
|
|---|
| 45 | nr=0
|
|---|
| 46 | sourcepos = []
|
|---|
| 47 | sourcemax = [0., 0., 0.]
|
|---|
| 48 | for line in config:
|
|---|
| 49 | nr=nr+1
|
|---|
| 50 | if source > nr:
|
|---|
| 51 | continue
|
|---|
| 52 | else:
|
|---|
| 53 | if re.compile("Ion_Type[0-9]+[_]+").match(line):
|
|---|
| 54 | words = line.split()
|
|---|
| 55 | print words
|
|---|
| 56 | sourcepos.append(words[0]+"="+words[1]+"="+words[2]+"="+words[3])
|
|---|
| 57 | for j in range(3):
|
|---|
| 58 | if sourcemax[j] < float(words[j+1]):
|
|---|
| 59 | sourcemax[j] = float(words[j+1])
|
|---|
| 60 | n=n+1
|
|---|
| 61 | else:
|
|---|
| 62 | if "#" in line: # we break on the next comment line after source
|
|---|
| 63 | break
|
|---|
| 64 | print "I found ",n," atoms in source config."
|
|---|
| 65 |
|
|---|
| 66 | # parse the target configuration
|
|---|
| 67 | config.seek(0)
|
|---|
| 68 | n=0
|
|---|
| 69 | nr=0
|
|---|
| 70 | targetpos = []
|
|---|
| 71 | targetmax= [0., 0., 0.]
|
|---|
| 72 | for line in config:
|
|---|
| 73 | nr=nr+1
|
|---|
| 74 | if target > nr:
|
|---|
| 75 | continue
|
|---|
| 76 | else:
|
|---|
| 77 | if re.compile("Ion_Type[0-9]+[_]+").match(line):
|
|---|
| 78 | words = line.split()
|
|---|
| 79 | print words
|
|---|
| 80 | targetpos.append(words[0]+"="+words[1]+"="+words[2]+"="+words[3])
|
|---|
| 81 | for j in range(3):
|
|---|
| 82 | if targetmax[j] < float(words[j+1]):
|
|---|
| 83 | targetmax[j] = float(words[j+1])
|
|---|
| 84 | n=n+1
|
|---|
| 85 | else:
|
|---|
| 86 | if "#" in line: # we break on the next comment line after source
|
|---|
| 87 | break
|
|---|
| 88 | print "I found ",n," atoms in target config."
|
|---|
| 89 | targetmax[2] = 2.64357
|
|---|
| 90 | # Now as we have the two configs, we need go through the steps desired and write the output files
|
|---|
| 91 | diff = [0., 0., 0.]
|
|---|
| 92 | BoxLength = [0., 0., 0., 0., 0., 0.]
|
|---|
| 93 | for i in range(steps+1):
|
|---|
| 94 | # change the boxlength
|
|---|
| 95 | if i == 0:
|
|---|
| 96 | for j in range(6):
|
|---|
| 97 | BoxLength[j] = boxlength[j]
|
|---|
| 98 | else:
|
|---|
| 99 | BoxLength[0] = boxlength[0] + boxlength[0]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
|
|---|
| 100 | BoxLength[1] = boxlength[1] + boxlength[1]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
|
|---|
| 101 | BoxLength[2] = boxlength[2] + boxlength[2]*(targetmax[1]/sourcemax[1]-1.)*(float(i)/float(steps))
|
|---|
| 102 | BoxLength[3] = boxlength[3] + boxlength[3]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
|
|---|
| 103 | BoxLength[4] = boxlength[4] + boxlength[4]*(targetmax[0]/sourcemax[0]-1.)*(float(i)/float(steps))
|
|---|
| 104 | BoxLength[5] = boxlength[5] + boxlength[5]*(targetmax[2]/sourcemax[2]-1.)*(float(i)/float(steps))
|
|---|
| 105 | # write the initial config file to the new destination
|
|---|
| 106 | newdir="step%02d" % (i)
|
|---|
| 107 | if not os.path.isdir(newdir):
|
|---|
| 108 | if os.path.isfile(newdir):
|
|---|
| 109 | raise OSError("a file with the same name as the desired " \
|
|---|
| 110 | "dir, '%s', already exists." % newdir)
|
|---|
| 111 | else:
|
|---|
| 112 | os.mkdir(newdir)
|
|---|
| 113 | output = file("step%02d/%s" % (i,sys.argv[1]),"w")
|
|---|
| 114 | config.seek(0)
|
|---|
| 115 | nr=0
|
|---|
| 116 | for line in config:
|
|---|
| 117 | nr=nr+1
|
|---|
| 118 | if box < nr and box+4 > nr:
|
|---|
| 119 | j=nr-box-1
|
|---|
| 120 | for m in range(j+1):
|
|---|
| 121 | if j==2:
|
|---|
| 122 | j=j+1
|
|---|
| 123 | output.write("%lg\t" % (BoxLength[j+m]))
|
|---|
| 124 | else:
|
|---|
| 125 | output.write("%lg\t" % (BoxLength[j+m]))
|
|---|
| 126 | output.write("\n")
|
|---|
| 127 | continue
|
|---|
| 128 | output.write(line)
|
|---|
| 129 | if source <= nr+1:
|
|---|
| 130 | break
|
|---|
| 131 | # add the shifted coordinates
|
|---|
| 132 | for j in range(n):
|
|---|
| 133 | sourceatoms = sourcepos[j].split("=")
|
|---|
| 134 | targetatoms = targetpos[j].split("=")
|
|---|
| 135 | for m in range(3):
|
|---|
| 136 | diff[m] = (float(targetatoms[m+1])-float(sourceatoms[m+1]))*float(i)/float(steps)
|
|---|
| 137 | 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]))
|
|---|
| 138 | output.close()
|
|---|
| 139 | config.close()
|
|---|
