| [e1a46d] | 1 | #!@PYTHON@
|
|---|
| 2 | #
|
|---|
| 3 | # This code calculates the pair correlation between the water molecules and the non-convex surface of the cluster (i.e. the boundary defined by the Si and Ca atoms)
|
|---|
| 4 |
|
|---|
| 5 | import sys, string, re, math, os
|
|---|
| 6 | wrerr = sys.stderr.write
|
|---|
| 7 | wrout = sys.stdout.write
|
|---|
| 8 |
|
|---|
| 9 | if (len(sys.argv) < 4):
|
|---|
| 10 | print "Usage: "+sys.argv[0]+" <xyz> <dbond> <dat>"
|
|---|
| 11 | sys.exit(1)
|
|---|
| 12 |
|
|---|
| 13 | # parse the pdb file and obtain maxid
|
|---|
| 14 | print "Parsing the XYZ file "+sys.argv[1]+"."
|
|---|
| 15 | xyzfile = open(sys.argv[1], "r")
|
|---|
| 16 | atoms = {}
|
|---|
| 17 | id=1
|
|---|
| 18 | xyzfile.readline()
|
|---|
| 19 | xyzfile.readline()
|
|---|
| 20 | emptyline=re.compile('[ \t]*\n')
|
|---|
| 21 | for line in xyzfile:
|
|---|
| 22 | #if "END" in line or "CONECT" in line or "MASTER" in line:
|
|---|
| 23 | if emptyline.match(line):
|
|---|
| 24 | continue
|
|---|
| 25 | fields = line.split()
|
|---|
| 26 | # id, element, x, y, z and ten neighbour ids
|
|---|
| 27 | atoms[str(id)] = [id, fields[0], float(fields[1]), float(fields[2]), float(fields[3]), -1, -1, -1, -1, -1, -1, -1, -1, -1, -1]
|
|---|
| 28 | id+=1
|
|---|
| 29 | print "The biggest Id is %d" % (id)
|
|---|
| 30 | xyzfile.close()
|
|---|
| 31 |
|
|---|
| 32 | # Parsing the dbond file
|
|---|
| 33 | print "Parsing the DBOND file "+sys.argv[2]+"."
|
|---|
| 34 | dbondfile = open(sys.argv[2], "r")
|
|---|
| 35 | test = dbondfile.readline() # skip the first two lines
|
|---|
| 36 | test = dbondfile.readline()
|
|---|
| 37 | for line in dbondfile:
|
|---|
| 38 | fields = line.split()
|
|---|
| 39 | ids = [int(fields[0]), int(fields[1])]
|
|---|
| 40 | set = [0,0]
|
|---|
| 41 | if str(ids[0]) in atoms and str(ids[1]) in atoms:
|
|---|
| 42 | for i in range(5,15):
|
|---|
| 43 | for j in range(2):
|
|---|
| 44 | if set[j] == 0 and (atoms[ str(ids[j]) ])[i] == -1:
|
|---|
| 45 | (atoms[ str(ids[j]) ])[i] = ids [ (j+1) %2 ]
|
|---|
| 46 | set[j] = 1
|
|---|
| 47 | if set[0] == 1 and set[1] == 1:
|
|---|
| 48 | break
|
|---|
| 49 | if set[0] == 0 or set[1] == 0:
|
|---|
| 50 | print "ERROR: Either id already had ten neighbours, increase number in code!"
|
|---|
| 51 | sys.exit(1)
|
|---|
| 52 | else:
|
|---|
| 53 | print "ERROR: Either %d or %d are not listed in the PDB file!" % (ids[0], ids[1])
|
|---|
| 54 | sys.exit(1)
|
|---|
| 55 |
|
|---|
| 56 | # parse the tecplot dat file (the first three lines are tecplot INFO)
|
|---|
| 57 | datfile=open(sys.argv[3], "r")
|
|---|
| 58 | line=datfile.readline()
|
|---|
| 59 | line=datfile.readline()
|
|---|
| 60 | line=datfile.readline()
|
|---|
| 61 | points=[]
|
|---|
| 62 | pcounter=0
|
|---|
| 63 | for line in datfile:
|
|---|
| 64 | if not line.strip(): # end on first empty line, indicating end of points and beginning of triangles
|
|---|
| 65 | break
|
|---|
| 66 | else:
|
|---|
| 67 | fields=line.split()
|
|---|
| 68 | points.append([pcounter, float(fields[0]), float(fields[1]), float(fields[2])])
|
|---|
| 69 | pcounter=pcounter+1
|
|---|
| 70 | print "There are %d endpoints of triangles." % (pcounter)
|
|---|
| 71 |
|
|---|
| 72 | tcounter=0
|
|---|
| 73 | triangles=[]
|
|---|
| 74 | for line in datfile:
|
|---|
| 75 | fields=line.split()
|
|---|
| 76 | triangles.append([tcounter, int(fields[0]), int(fields[1]), int(fields[2])])
|
|---|
| 77 | corners=[points[ int(fields[0])-1 ], points[ int(fields[1])-1 ], points[ int(fields[2])-1 ]]
|
|---|
| 78 | # calculate middle point and add to points
|
|---|
| 79 | centroid=[0.,0.,0.]
|
|---|
| 80 | for i in range(3):
|
|---|
| 81 | for j in range(3):
|
|---|
| 82 | centroid[j] += corners[i][1+j]
|
|---|
| 83 | for i in range(3):
|
|---|
| 84 | centroid[i] = centroid[i]/3.
|
|---|
| 85 | points.append([pcounter+tcounter, centroid[0], centroid[1], centroid[2]])
|
|---|
| 86 | #print ("Centroid is at (%g,%g,%g) for triangle at (%g,%g,%g), (%g,%g,%g), (%g,%g,%g)." % (centroid[0], centroid[1], centroid[2], corners[0][1],corners[0][2],corners[0][3], corners[1][1],corners[1][2],corners[1][3], corners[2][1],corners[2][2],corners[2][3]))
|
|---|
| 87 |
|
|---|
| 88 | tcounter=tcounter+1
|
|---|
| 89 | print "There are %d triangles." % (tcounter)
|
|---|
| 90 |
|
|---|
| 91 | # go through all Si and Ca and calculate CoG
|
|---|
| 92 | CoG = [ 0, 0, 0]
|
|---|
| 93 | key=atoms.keys()[0]
|
|---|
| 94 | max=[atoms[key][2],atoms[key][3],atoms[key][4]]
|
|---|
| 95 | min=[atoms[key][2],atoms[key][3],atoms[key][4]]
|
|---|
| 96 | nr=0
|
|---|
| 97 | for atom in atoms:
|
|---|
| 98 | #print atoms[atom]
|
|---|
| 99 | if "Si" in atoms[atom][1] or "Ca" in atoms[atom][1]:
|
|---|
| 100 | for i in range(3):
|
|---|
| 101 | CoG[i] += atoms[atom][2+i]
|
|---|
| 102 | nr = nr + 1
|
|---|
| 103 | for i in range(3):
|
|---|
| 104 | if min[i] > atoms[atom][2+i]:
|
|---|
| 105 | min[i] = atoms[atom][2+i]
|
|---|
| 106 | if max[i] < atoms[atom][2+i]:
|
|---|
| 107 | max[i] = atoms[atom][2+i]
|
|---|
| 108 | for i in range(3):
|
|---|
| 109 | CoG[i] /= nr
|
|---|
| 110 | print "The Center of Gravity is at %g,%g,%g." % (CoG[0], CoG[1], CoG[2])
|
|---|
| 111 | print "Bounds of the cluster's atoms are [%f,%f,], [%f,%f], [%f,%f]." % (min[0],max[0], min[1],max[1], min[2],max[2])
|
|---|
| 112 |
|
|---|
| 113 | # check bounds against bounds of cluster
|
|---|
| 114 | max=[points[0][1],points[0][2],points[0][3]]
|
|---|
| 115 | min=[points[0][1],points[0][2],points[0][3]]
|
|---|
| 116 | for j in range(len(points)):
|
|---|
| 117 | for i in range(3):
|
|---|
| 118 | if min[i] > points[j][1+i]:
|
|---|
| 119 | min[i] = points[j][1+i]
|
|---|
| 120 | if max[i] < points[j][1+i]:
|
|---|
| 121 | max[i] = points[j][1+i]
|
|---|
| 122 | print "Bounds of the cluster's tesselated surface are [%f,%f,], [%f,%f], [%f,%f]." % (min[0],max[0], min[1],max[1], min[2],max[2])
|
|---|
| 123 |
|
|---|
| 124 | # find radius of cluster
|
|---|
| 125 | bounds=[0.,0.,0.]
|
|---|
| 126 | maximum=0
|
|---|
| 127 | maxid=-1
|
|---|
| 128 | for i in range(3):
|
|---|
| 129 | bounds[i] = max[i] - min[i]
|
|---|
| 130 | if maximum < bounds[i]:
|
|---|
| 131 | maxid=i
|
|---|
| 132 | maximum = bounds[i]
|
|---|
| 133 | print "Cylinder direction is probably "+str(maxid)+"."
|
|---|
| 134 | R=(bounds[(maxid+1)%3]+bounds[(maxid+2)%3])/4.
|
|---|
| 135 |
|
|---|
| 136 | # search atom closest to CoG to allow visual check
|
|---|
| 137 | mindist=1e+16
|
|---|
| 138 | minid=-1
|
|---|
| 139 | for atom in atoms:
|
|---|
| 140 | dist=0
|
|---|
| 141 | for i in range(3):
|
|---|
| 142 | temp = CoG[i] - atoms[atom][2+i]
|
|---|
| 143 | dist += temp*temp
|
|---|
| 144 | if dist < mindist:
|
|---|
| 145 | mindist = dist
|
|---|
| 146 | minid = atom
|
|---|
| 147 | print "The atom closest to CoG is %d, element %s (%f,%f,%f) with distance %f." % (atoms[minid][0], atoms[minid][1], atoms[minid][2], atoms[minid][3], atoms[minid][4], math.sqrt(mindist))
|
|---|
| 148 |
|
|---|
| 149 | # initialize the bin array
|
|---|
| 150 | NrBins=100
|
|---|
| 151 | BinWidth=0.5
|
|---|
| 152 | bins = []
|
|---|
| 153 | for i in range(2*NrBins+1):
|
|---|
| 154 | bins.append(0)
|
|---|
| 155 |
|
|---|
| 156 | # get path
|
|---|
| 157 | line=sys.argv[1].rpartition('/')
|
|---|
| 158 | path=line[0]+line[1]
|
|---|
| 159 | print "Path is "+path+"."
|
|---|
| 160 |
|
|---|
| 161 | # open a number of xyz files
|
|---|
| 162 | innerwaterfile=open(path+"innerwater","w")
|
|---|
| 163 | outerwaterfile=open(path+"outerwater","w")
|
|---|
| 164 | hyperdryclusterfile=open(path+"hyperdrycluster","w")
|
|---|
| 165 | dryclusterfile=open(path+"drycluster","w")
|
|---|
| 166 | waterdryclusterfile=open(path+"waterdrycluster","w")
|
|---|
| 167 | layerclusterfile=open(path+"layercluster","w")
|
|---|
| 168 | outerwaterfile=open(path+"outerwater","w")
|
|---|
| 169 | nrcluster=0
|
|---|
| 170 | for atom in atoms:
|
|---|
| 171 | if "Si" in atoms[atom][1] or "Ca" in atoms[atom][1]:
|
|---|
| 172 | #innerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 173 | #outerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 174 | hyperdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 175 | dryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 176 | waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 177 | layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 178 | nrcluster+=1
|
|---|
| 179 | nohyperdry=nrcluster
|
|---|
| 180 | nodry=nrcluster
|
|---|
| 181 | nowaterdry=nrcluster
|
|---|
| 182 | nolayer=nrcluster
|
|---|
| 183 |
|
|---|
| 184 | # go through every atom
|
|---|
| 185 | elements=["H", "O"]
|
|---|
| 186 | nrwaters = 0
|
|---|
| 187 | nrothers = 0
|
|---|
| 188 | inner=0
|
|---|
| 189 | outer=0
|
|---|
| 190 | for element in elements:
|
|---|
| 191 | print "Calculating for element "+element+"."
|
|---|
| 192 | for i in range(2*NrBins+1):
|
|---|
| 193 | bins[i] = 0.
|
|---|
| 194 | for atom in atoms:
|
|---|
| 195 | if element in atoms[atom][1]:
|
|---|
| 196 | water = 1
|
|---|
| 197 | # check whether its really water (i.e. not bound to Si)
|
|---|
| 198 | for i in range(5,15):
|
|---|
| 199 | if atoms[atom][i] == -1:
|
|---|
| 200 | break
|
|---|
| 201 | if str(atoms[atom][i]) in atoms:
|
|---|
| 202 | if "Si" in atoms[ str(atoms[atom][i]) ][1] or "Ca" in atoms[ str(atoms[atom][i]) ][1]:
|
|---|
| 203 | water = 0
|
|---|
| 204 | break
|
|---|
| 205 | else:
|
|---|
| 206 | print "ERROR: key %d is not in list!" % (atoms[atom][i])
|
|---|
| 207 | sys.exit(1)
|
|---|
| 208 | # if it's water find closest boundary point and check whether atom is in- or outside of cluster
|
|---|
| 209 | if water == 1:
|
|---|
| 210 | mindist=1e+16
|
|---|
| 211 | vector=[0.,0.,0.]
|
|---|
| 212 | vector2=[0.,0.,0.]
|
|---|
| 213 | minptr=points[0]
|
|---|
| 214 | for ptr in points:
|
|---|
| 215 | dist = 0
|
|---|
| 216 | for i in range(3):
|
|---|
| 217 | vector[i] = ptr[1+i] - atoms[atom][2+i]
|
|---|
| 218 | dist += (vector[i])*(vector[i])
|
|---|
| 219 | if dist < mindist:
|
|---|
| 220 | minptr=ptr
|
|---|
| 221 | mindist = dist
|
|---|
| 222 | sign=1
|
|---|
| 223 | skp=0
|
|---|
| 224 | diff=0
|
|---|
| 225 | for i in range(3):
|
|---|
| 226 | vector2[i] = CoG[i] - atoms[atom][2+i]
|
|---|
| 227 | skp += vector[i]*vector2[i]
|
|---|
| 228 | diff += vector2[i]*vector2[i]
|
|---|
| 229 | if (skp < 0) or (diff < dist):
|
|---|
| 230 | sign=-1
|
|---|
| 231 | mindist = int(math.sqrt(mindist)/BinWidth)
|
|---|
| 232 | if mindist < NrBins:
|
|---|
| 233 | #if mindist < 8 and abs(sign) == 1:
|
|---|
| 234 | #for i in range(3):
|
|---|
| 235 | #vector[i] = minptr[1+i] - atoms[atom][2+i]
|
|---|
| 236 | #vector2[i] = CoG[i] - atoms[atom][2+i]
|
|---|
| 237 | #print "For atom %d vector to boundary is (%f,%f,%f) and vector to CoG (%f,%f,%f)." % (atoms[atom][0], vector[0], vector[1], vector[2], vector2[0], vector2[1], vector2[2])
|
|---|
| 238 | #print "atom %d vector is (%f,%f,%f) and boundary point vector is (%f,%f,%f)." % (atoms[atom][0], atoms[atom][2], atoms[atom][3], atoms[atom][4],minptr[1], minptr[2], minptr[3])
|
|---|
| 239 | #print "Atom %d goes into bin %d." % (atoms[atom][0], mindist)
|
|---|
| 240 | if sign == -1:
|
|---|
| 241 | bins[NrBins-mindist] = bins[NrBins-mindist] + 1
|
|---|
| 242 | innerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 243 | dryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 244 | waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 245 | layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 246 | nodry += 1
|
|---|
| 247 | nowaterdry += 1
|
|---|
| 248 | nolayer += 1
|
|---|
| 249 | inner += 1
|
|---|
| 250 | else:
|
|---|
| 251 | bins[NrBins+mindist] = bins[NrBins+mindist] + 1
|
|---|
| 252 | outerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 253 | if mindist < 8:
|
|---|
| 254 | waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 255 | nowaterdry += 1
|
|---|
| 256 | if mindist < 11:
|
|---|
| 257 | layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 258 | nolayer += 1
|
|---|
| 259 | outer += 1
|
|---|
| 260 | nrwaters = nrwaters+1
|
|---|
| 261 | else:
|
|---|
| 262 | print "NOTE: Have to throw "+str(atoms[atom][0])+" away, as it would be put into Bin "+str(mindist)+" which is out of range 0,"+str(NrBins-1)+"."
|
|---|
| 263 | else:
|
|---|
| 264 | if element == elements[1]:
|
|---|
| 265 | nrothers = nrothers + 1
|
|---|
| 266 | nrcluster = nrcluster + 1
|
|---|
| 267 | #innerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 268 | #outerwaterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 269 | dryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 270 | hyperdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 271 | waterdryclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 272 | layerclusterfile.write("%s\t%f\t%f\t%f\n" % (atoms[atom][1], atoms[atom][2], atoms[atom][3], atoms[atom][4]))
|
|---|
| 273 | nodry += 1
|
|---|
| 274 | nohyperdry += 1
|
|---|
| 275 | nowaterdry += 1
|
|---|
| 276 | nolayer += 1
|
|---|
| 277 |
|
|---|
| 278 | #print "The Bins are as follows:"
|
|---|
| 279 | output = open(path+"histo_pair_correlation_"+element+".dat", "w")
|
|---|
| 280 | output.write("#Distance [A]\tnumber in bin\n")
|
|---|
| 281 | for i in range(NrBins):
|
|---|
| 282 | #print "%g\t%d" % (i*BinWidth, bins[i])
|
|---|
| 283 | output.write("%g\t%d\n" % (-(NrBins-i)*BinWidth, bins[i]))
|
|---|
| 284 | for i in range(NrBins,2*NrBins):
|
|---|
| 285 | #print "%g\t%d" % (i*BinWidth, bins[i])
|
|---|
| 286 | output.write("%g\t%d\n" % ((i-NrBins)*BinWidth, bins[i]))
|
|---|
| 287 | output.close()
|
|---|
| 288 |
|
|---|
| 289 |
|
|---|
| 290 | # calculate numerically the bins for homogeneous water distribution where we assume the cluster to be a cylinder (hence, just 2d and circle)
|
|---|
| 291 | # centered in the box of length 44 with a certain radius R
|
|---|
| 292 | a=22.0 # 44/2
|
|---|
| 293 | N=2200
|
|---|
| 294 | lastarea=0
|
|---|
| 295 | waterdensity=0.033456344 # (1g/18amu)/cm^3 -> 1/angstrom^3, i.e. number of hydrogen particles in Angstroem^3 volume
|
|---|
| 296 | if "H" in element:
|
|---|
| 297 | waterdensity = waterdensity*2
|
|---|
| 298 | output=open(path+"homogeneous_water_histo_"+element+".dat","w")
|
|---|
| 299 | for i in range(1,2*NrBins):
|
|---|
| 300 | area=0
|
|---|
| 301 | radius = float(i)*BinWidth
|
|---|
| 302 | h=(a)/float(N)
|
|---|
| 303 | for j in range(N):
|
|---|
| 304 | x = j*h
|
|---|
| 305 | if x < a and x < radius:
|
|---|
| 306 | y = (radius)*(radius) - x*x
|
|---|
| 307 | if y > a*a:
|
|---|
| 308 | y=a*a
|
|---|
| 309 | dx = h
|
|---|
| 310 | if (x+dx > a):
|
|---|
| 311 | dx -= a-x
|
|---|
| 312 | if (x+dx > radius):
|
|---|
| 313 | dx -= radius-x
|
|---|
| 314 | if (dx < 0):
|
|---|
| 315 | print "ERROR: dx is "+str(dx)+"."
|
|---|
| 316 | area += dx * math.sqrt(y)
|
|---|
| 317 | #print ("Area for radius %g is %g." % (radius, area-lastarea))
|
|---|
| 318 | volume = int((2.*a) * 4.* (area - lastarea) * waterdensity) # 2a is length of cylinder, 4* due to only looking at a quadrant
|
|---|
| 319 | lastarea = area
|
|---|
| 320 | if radius <= R:
|
|---|
| 321 | volume = 0
|
|---|
| 322 | output.write("%g\t%d\n" % ((i*BinWidth)-R, volume))
|
|---|
| 323 | output.close()
|
|---|
| 324 |
|
|---|
| 325 | print "There were "+str(nrwaters)+" hydrogen and oxygen atoms and "+str(nrothers)+" connected to silica."
|
|---|
| 326 | print "There were "+str(inner)+" inner and "+str(outer)+" outer hydrogen and oxygen atoms."
|
|---|
| 327 |
|
|---|
| 328 | innerwaterfile.close()
|
|---|
| 329 | outerwaterfile.close()
|
|---|
| 330 | hyperdryclusterfile.close()
|
|---|
| 331 | dryclusterfile.close()
|
|---|
| 332 | waterdryclusterfile.close()
|
|---|
| 333 | layerclusterfile.close()
|
|---|
| 334 |
|
|---|
| 335 | # put number of atoms in front
|
|---|
| 336 | innerwaterfile=open(path+"innerwater.xyz","w")
|
|---|
| 337 | innerwaterfileold=open(path+"innerwater","r")
|
|---|
| 338 | innerwaterfile.write(str(inner)+"\n\tCreated by "+sys.argv[0]+"\n")
|
|---|
| 339 | for line in innerwaterfileold:
|
|---|
| 340 | innerwaterfile.write(line)
|
|---|
| 341 | innerwaterfile.close()
|
|---|
| 342 | innerwaterfileold.close()
|
|---|
| 343 | os.remove(path+"innerwater")
|
|---|
| 344 |
|
|---|
| 345 | outerwaterfile=open(path+"outerwater.xyz","w")
|
|---|
| 346 | outerwaterfileold=open(path+"outerwater","r")
|
|---|
| 347 | outerwaterfile.write(str(outer)+"\n\tCreated by "+sys.argv[0]+"\n")
|
|---|
| 348 | for line in outerwaterfileold:
|
|---|
| 349 | outerwaterfile.write(line)
|
|---|
| 350 | outerwaterfile.close()
|
|---|
| 351 | outerwaterfileold.close()
|
|---|
| 352 | os.remove(path+"outerwater")
|
|---|
| 353 |
|
|---|
| 354 | hyperdryclusterfile=open(path+"hyperdrycluster.xyz","w")
|
|---|
| 355 | hyperdryclusterfileold=open(path+"hyperdrycluster","r")
|
|---|
| 356 | hyperdryclusterfile.write(str(nohyperdry)+"\n\tCreated by "+sys.argv[0]+"\n")
|
|---|
| 357 | for line in hyperdryclusterfileold:
|
|---|
| 358 | hyperdryclusterfile.write(line)
|
|---|
| 359 | hyperdryclusterfile.close()
|
|---|
| 360 | hyperdryclusterfileold.close()
|
|---|
| 361 | os.remove(path+"hyperdrycluster")
|
|---|
| 362 |
|
|---|
| 363 | dryclusterfile=open(path+"drycluster.xyz","w")
|
|---|
| 364 | dryclusterfileold=open(path+"drycluster","r")
|
|---|
| 365 | dryclusterfile.write(str(nodry)+"\n\tCreated by "+sys.argv[0]+"\n")
|
|---|
| 366 | for line in dryclusterfileold:
|
|---|
| 367 | dryclusterfile.write(line)
|
|---|
| 368 | dryclusterfile.close()
|
|---|
| 369 | dryclusterfileold.close()
|
|---|
| 370 | os.remove(path+"drycluster")
|
|---|
| 371 |
|
|---|
| 372 | waterdryclusterfile=open(path+"waterdrycluster.xyz","w")
|
|---|
| 373 | waterdryclusterfileold=open(path+"waterdrycluster","r")
|
|---|
| 374 | waterdryclusterfile.write(str(nowaterdry)+"\n\tCreated by "+sys.argv[0]+"\n")
|
|---|
| 375 | for line in waterdryclusterfileold:
|
|---|
| 376 | waterdryclusterfile.write(line)
|
|---|
| 377 | waterdryclusterfile.close()
|
|---|
| 378 | waterdryclusterfileold.close()
|
|---|
| 379 | os.remove(path+"waterdrycluster")
|
|---|
| 380 |
|
|---|
| 381 | layerclusterfile=open(path+"layercluster.xyz","w")
|
|---|
| 382 | layerclusterfileold=open(path+"layercluster","r")
|
|---|
| 383 | layerclusterfile.write(str(nolayer)+"\n\tCreated by "+sys.argv[0]+"\n")
|
|---|
| 384 | for line in layerclusterfileold:
|
|---|
| 385 | layerclusterfile.write(line)
|
|---|
| 386 | layerclusterfile.close()
|
|---|
| 387 | layerclusterfileold.close()
|
|---|
| 388 | os.remove(path+"layercluster")
|
|---|
| 389 |
|
|---|
| 390 | print "finished."
|
|---|
