| 1 | import pyMoleCuilder as mol
|
|---|
| 2 | import sys, os, math
|
|---|
| 3 | try:
|
|---|
| 4 | import numpy
|
|---|
| 5 | numpy_present = True
|
|---|
| 6 | except ImportError:
|
|---|
| 7 | numpy_present = False
|
|---|
| 8 |
|
|---|
| 9 | if len(sys.argv) < 5:
|
|---|
| 10 | print('Usage: '+sys.argv[0]+' <input> <path> <steps> <no_atoms> <use bondgraph>')
|
|---|
| 11 | sys.exit(1)
|
|---|
| 12 |
|
|---|
| 13 | steps=int(sys.argv[3])
|
|---|
| 14 | equilibrium_distance=1.6
|
|---|
| 15 | no_atoms=int(sys.argv[4])
|
|---|
| 16 | inputfile=sys.argv[1]
|
|---|
| 17 | forcespath=sys.argv[2]
|
|---|
| 18 | forcesfile="ising.forces"
|
|---|
| 19 | use_bondgraph=sys.argv[5]
|
|---|
| 20 |
|
|---|
| 21 | # creating input file
|
|---|
| 22 | atomstart=7.6-1.6*math.floor(no_atoms/2)
|
|---|
| 23 | print("Creating "+inputfile)
|
|---|
| 24 | with open(inputfile, 'w') as f:
|
|---|
| 25 | f.write("# ATOMDATA\ttype\tId\tx=3\tu=3\tF=3\tneighbors=4\n")
|
|---|
| 26 | f.write("# Box\t20\t0\t0\t0\t20\t0\t0\t0\t20\n")
|
|---|
| 27 | for i in range(1, no_atoms+1):
|
|---|
| 28 | atompos=atomstart+1.6*float(i)
|
|---|
| 29 | if i==math.floor(no_atoms/2+1):
|
|---|
| 30 | atompos=atompos-.5
|
|---|
| 31 | if i==1:
|
|---|
| 32 | f.write("C\t%d\t%lg\t10\t10\t0\t0\t0\t0\t0\t0\t%d\t0\t0\t0\n" % (i, atompos, i+1));
|
|---|
| 33 | elif i==no_atoms:
|
|---|
| 34 | f.write("C\t%d\t%lg\t10\t10\t0\t0\t0\t0\t0\t0\t%d\t0\t0\t0\n" % (i, atompos, i-1));
|
|---|
| 35 | else:
|
|---|
| 36 | f.write("C\t%d\t%lg\t10\t10\t0\t0\t0\t0\t0\t0\t%d\t%d\t0\t0\n" % (i, atompos, i-1, i+1));
|
|---|
| 37 |
|
|---|
| 38 | print("Parsing from "+inputfile)
|
|---|
| 39 | mol.WorldInput(inputfile)
|
|---|
| 40 | mol.SelectionAllAtoms()
|
|---|
| 41 | mol.CommandVerbose("4")
|
|---|
| 42 |
|
|---|
| 43 | # calculate damping factor from finite geometric series
|
|---|
| 44 | # s_n/a = \sum^{n-1}_{k=0} r^k = (1-r^n)/(1-r) -> s_(n+1)/a -1 = \sum^{n}_{k=1} r^k = (1-r^(n+1))/(1-r) - 1
|
|---|
| 45 | # \sum^{n}_{k=1} r^k := 1 and 1 = (1-r^(n+1))/(1-r) - 1 -> 2*(1-r) = 1 - r^(n+1) -> 1 - 2*r + r^(n+1) = 0
|
|---|
| 46 | # find root: p[0] is coefficient of monomial with highest power
|
|---|
| 47 | if numpy_present:
|
|---|
| 48 | p=[0.] * (no_atoms+1)
|
|---|
| 49 | p[0]=1.
|
|---|
| 50 | p[no_atoms-1]=-2.
|
|---|
| 51 | p[no_atoms]=1.
|
|---|
| 52 | zeros=numpy.roots(p)
|
|---|
| 53 | print("Roots of p "+str(p)+" are "+str(zeros))
|
|---|
| 54 | damping=numpy.real(zeros[-1])
|
|---|
| 55 | print("Using damping factor of "+str(damping))
|
|---|
| 56 | else:
|
|---|
| 57 | if no_atoms == 2:
|
|---|
| 58 | damping=0.5
|
|---|
| 59 | elif no_atoms == 5:
|
|---|
| 60 | damping=0.5436890126920764
|
|---|
| 61 |
|
|---|
| 62 | for i in range(0, steps):
|
|---|
| 63 | # TODO: Python interface should have something to iterate over selected atoms
|
|---|
| 64 | # and molecules and get information on their internal status
|
|---|
| 65 |
|
|---|
| 66 | # read current atomic positions
|
|---|
| 67 | outputfile=forcespath+'/'+forcesfile+'.xyz'
|
|---|
| 68 | try:
|
|---|
| 69 | os.remove(outputfile)
|
|---|
| 70 | except: OSError
|
|---|
| 71 | #
|
|---|
| 72 | mol.WorldOutputAs(outputfile)
|
|---|
| 73 | mol.wait()
|
|---|
| 74 | distances=[]
|
|---|
| 75 | coords=[0.,0.,0.]
|
|---|
| 76 | try:
|
|---|
| 77 | skiplines=2+i*(1+1+no_atoms+1) # no_atoms, comment, no_atoms atoms, empty line
|
|---|
| 78 | with open(outputfile) as f:
|
|---|
| 79 | for line in f:
|
|---|
| 80 | if skiplines != 0:
|
|---|
| 81 | skiplines=skiplines-1
|
|---|
| 82 | continue
|
|---|
| 83 | line=line.replace('\t',' ')
|
|---|
| 84 | print("LINE: "+line)
|
|---|
| 85 | [elementtype, X, Y, Z] = line.split(' ', 4)
|
|---|
| 86 | if coords!=[0.,0.,0.]:
|
|---|
| 87 | distances.append(math.sqrt((coords[0]-float(X))**2+(coords[1]-float(Y))**2+(coords[2]-float(Z))**2))
|
|---|
| 88 | coords=[float(X),float(Y),float(Z)]
|
|---|
| 89 | except IOError:
|
|---|
| 90 | print('Warning: '+outputfile+' not readable.')
|
|---|
| 91 | sys.exit(1)
|
|---|
| 92 |
|
|---|
| 93 | assert(len(distances)==no_atoms-1)
|
|---|
| 94 |
|
|---|
| 95 | #
|
|---|
| 96 | # generate Ising model forces and store in file
|
|---|
| 97 | #
|
|---|
| 98 | # i.e. we have spring forces between neighboring atoms depending on their distance
|
|---|
| 99 | forces=[]
|
|---|
| 100 | for d in distances:
|
|---|
| 101 | forces.append( d - equilibrium_distance );
|
|---|
| 102 |
|
|---|
| 103 | # generate new forces file
|
|---|
| 104 |
|
|---|
| 105 | with open(forcespath+'/'+forcesfile, 'w') as f:
|
|---|
| 106 | f.write('# atom\tf_x\tf_y\tf_z\n')
|
|---|
| 107 | for i in range(len(distances)+1):
|
|---|
| 108 | force=0
|
|---|
| 109 | if i!=0:
|
|---|
| 110 | force=force-forces[i-1]
|
|---|
| 111 | if (i != len(distances)):
|
|---|
| 112 | force=force+forces[i]
|
|---|
| 113 | f.write("%d\t%f\t0.\t0.\n" % (i+1, force))
|
|---|
| 114 |
|
|---|
| 115 | mol.WorldStepWorldTime("1")
|
|---|
| 116 | mol.MoleculeForceAnnealing(forcespath+'/'+forcesfile, ".1", "%d" % (steps), "%d" % (no_atoms-1), "%lg" % (damping), use_bondgraph)
|
|---|
| 117 | mol.wait()
|
|---|
| 118 |
|
|---|
| 119 | sys.exit(0)
|
|---|
