Commit fabd4009 authored by Matthias Schroeder's avatar Matthias Schroeder

initial release

parents
#!/usr/bin/python3
#
import numpy as np
import os
MAX_MASS = 2000
LINES_MOL = 2000
FILEPAT = "mol_%04d"
OUTDIR = 'masses/'
# description of elements
elements = [
{ "formel":'C', "weight":1.01, "min":1, "max":20 },
{ "formel":'N', "weight":1.02, "min":0, "max":50 },
{ "formel":'H', "weight":1.01, "min":1, "max":20 },
]
# create outdir if not exists
if not os.path.isdir(OUTDIR):
os.mkdir(OUTDIR)
# from formel to idx
element_idx= {}
for i in range(len(elements)):
element_idx[elements[i]['formel']] = i
# prepare some arrays
moltype = np.dtype([('strct', np.int16, len(elements)), ('weight', np.float64)])
mols= np.ndarray(shape=(MAX_MASS,LINES_MOL), dtype=moltype)
posi= [0] * MAX_MASS
def check_strct(s,w):
""" Check for relations between element, total mass, bla bla,
return True for seemingly valid element,
False otherwise
"""
if w > 3000:
return False
if s[element_idx['C']] > s[element_idx['H']]:
return False
return True
def format_strct(s):
""" Formats molecule in a nice way """
m=''
for i in range(len(elements)):
if s[0][i]>0:
m += elements[i]['formel'] + '_' + str(s[0][i]) + ' '
return "%16.10f\t%s" % (s[1], m)
def calc_weight(s):
""" return total weight of molekule descripte by s """
w=0
for i in range(len(elements)):
w= w + s[i] * elements[i]['weight']
return w
def flush_mass(w):
""" writes out masses for group w, w is integer for masses m
with w<= m < w+1
in case there is something """
if posi[w]>0:
filename= FILEPAT % w
print("Filename %s, cnt=%d" % (filename,posi[w]))
f = open(os.path.join(OUTDIR, filename), 'a')
for l in range(posi[w]):
f.write(format_strct(mols[w][l]) + "\n")
f.close()
posi[w]=0
if __name__ == "__main__":
# set start values
strct= [ e['min'] for e in elements ]
# start the work
while True:
w= calc_weight(strct)
if check_strct(strct, w):
# we have a valid molecule
widx= int(w)
mols[widx][posi[widx]]= (strct[:], w)
posi[widx]+= 1
if posi[widx]>=LINES_MOL:
flush_mass(widx)
# build next molecule
inc_done = False
for i in range(len(elements)):
if strct[i]>=elements[i]['max']:
# hit top, -> fall down and continue to increment next position
strct[i]=elements[i]['min']
else:
# increment next position and done
strct[i]= strct[i]+1
inc_done =True
break
# no increment - so we had had reached the last element
if not inc_done:
break
# flush the rest ....
for l in range(MAX_MASS):
flush_mass(l)
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment