source: utils/boxmaker.py@ 32bc47

import re, os, os.path, sys, operator

class c_opt():
    basename = None
    tremofiledir = './'
    potentialsfiledir = './'
    outfilename = 'out'

    source = None
    molarmass = None
    density = None
    temp = None

    number = '1000'
    
    cubicdomain = 'on'
    cubiccell = 'off'
    autorotate = 'off'
    autodim = 'on'
    postprocess = 'on'

    def update(self, name, value):
        if name in dir(self):
            exec('self.%s = "%s"' % (name, value))
        else:
            print 'Warning: Unknown option:', name


def ReadSettings(opt):
    # Obtain basename
    if len(sys.argv) >= 2:
        opt.basename = sys.argv[1]
    else:
        print 'Usage: boxmaker.py <basename> [options]'
        exit()

    # Read settings file
    with open('boxmaker.' + opt.basename) as f:
        for line in f:
            if len(line) > 0 and line[0] != '#':
                L, S, R = line.partition('=')
                opt.update(L.strip(), R.strip())

    # Parse parameters
    i = 2
    while i < len(sys.argv):
        L = sys.argv[i]

        if L[0] in '+-':
            LN = L[1:]

            if L[0] == '+':
                R = 'on'
            else:
                R = 'off'
        else:
            LN = L
            i += 1
            R = sys.argv[i]
        
        opt.update(LN, R)
        i += 1


def ReadUnits(opt):
    lines = [] # The file needs to be processed twice, so we save the lines in the first run

    with open(opt.tremofiledir + opt.basename + '.tremolo') as f:
        for line in f:
            if len(line) > 0 and line[0] != '#':
                line = line.strip()
                lines.append(line)

                if 'systemofunits' in line:
                    L, S, SOU = line.partition('=')
                    SOU = SOU.strip()[:-1] # Remove semicolon

    if SOU == 'custom':
        units = {}
        quantities = ['length', 'mass']

        for quantity in quantities:
            units[quantity] = [None, None] # Init with scaling factor and unit 'None'.

        for line in lines:
            for quantity in quantities:
                if quantity in line:
                    L, S, R = line.partition('=')
                    R = R.strip()[:-1] # Remove semicolon

                    if 'scalingfactor' in line:
                        units[quantity][0] = R
                    else:
                        units[quantity][1] = R

    elif SOU == 'kcalpermole':
        units = {'length': ['1', 'angstrom'], 'mass': ['1', 'u']}

    elif SOU == 'evolt':
        units = {'length': ['1', 'angstrom'], 'mass': ['1', 'u']}

    else: # SI
        units = {'length': ['1', 'm'], 'mass': ['1', 'kg']}

    return units


def ConvertUnits(have, want):
    # Redo with pipes?
    ret = os.system("units '%s' '%s' > temp_units_output" % (have, want))

    if ret == 0:
        with open('temp_units_output') as f:
            line = f.readline()

        os.system('rm temp_units_output')

        return float(line[3:-1])
    else:
        raise NameError('UnitError')


def UpdateSettings(opt):
    # Map boolean values
    for name in ['cubicdomain', 'cubiccell', 'autorotate', 'autodim', 'postprocess']:
        value = eval('opt.' + name)
        
        if value == 'on':
            value = True
        elif value == 'off':
            value = False
        else:
            print 'Not a boolean value:', value
            exit()
            
        exec('opt.' + name + '= value')

    # Convert dimensions
    if opt.autodim:
        units = ReadUnits(opt)

        have = opt.molarmass
        want = '%s*%s / mol' % tuple(units['mass'])
        opt.molarmass = ConvertUnits(have, want)

        have = opt.density
        want = '(%s*%s) ' % tuple(units['mass']) + '/ (%s*%s)**3' % tuple(units['length'])
        opt.density = ConvertUnits(have, want)
    else:
        opt.molarmass = float(opt.molarmass)
        opt.density = float(opt.density)

    # Number might be an integer or a 3-vector
    nvec = opt.number.split()
    if len(nvec) == 3:
        opt.number = [0]*3

        for i in range(3):
            opt.number[i] = int(nvec[i])
    else:
        opt.number = int(opt.number)
        
    # Temp should be float or None
    if opt.temp:
        opt.temp = float(opt.temp)


def FindBestCube(opt):
    newroot = int( round(opt.number**(1./3)) )
    newnumber = newroot**3

    if newnumber != opt.number:
        print 'Warning: Number changed to %d.' % newnumber

    return [newroot] * 3


def FindBestCuboid(opt):
    n = opt.number

    # Prime factors of n
    factors = []

    for i in [2, 3]:
        while n % i == 0:
            factors.append(i)
            n /= 2

    t = 5
    diff = 2

    while t*t <= n:
        while n % t == 0:
            factors.append(t)
            n /= t

        t = t + diff
        diff = 6 - diff

    if n > 1:
        factors.append(n)

    # Even distribution of current biggest prime to each vector -> similar sizes
    if len(factors) < 3:
        print 'Warning: Not enough prime factors - falling back to cubic placement'
        return FindBestCube(opt)

    factors.sort()
    distri = [[],[],[]]
    current = 0

    for factor in factors:
        distri[current].append(factor)
        current += 1
        if current == 3:
            current = 0

    result = [0]*3

    print '======== CUBOID USED:',

    for i in range(3):
        result[i] = int( reduce(operator.mul, distri[i]) )

    print result
    return result


def GetSourceBBabs(opt):
    bbmax = [0.0]*3
    bbmin = [float('inf')]*3

    s_name_ext = os.path.basename(opt.source).rsplit('.', 1)
    s_namepart = s_name_ext[0]

    if len(s_name_ext) > 1:
        s_ext = s_name_ext[1]
    else:
        s_ext = ''

    # Convert from any format to xyz
    os.system('molecuilder -o xyz --parse-tremolo-potentials %s -i temp_source.xyz -l %s' % (opt.potentialsfiledir+opt.basename+'.potentials', opt.source))

    # Calculate bounding box from xyz-file
    with open('temp_source.xyz') as f:
        N = int(f.readline())
        comment = f.readline()

        for i in xrange(N):
            buf = f.readline()
            xyz = buf.split()[1:]

            for i in range(3):
                bbmax[i] = max(bbmax[i], float(xyz[i]))
                bbmin[i] = min(bbmin[i], float(xyz[i]))

    bb = [0.0]*3

    for i in range(3):
        bb[i] = abs(bbmax[i] - bbmin[i])

    os.system('rm temp_source.*')
    return bb

# Global options with sensible default parameters
opt = c_opt()

ReadSettings(opt)
UpdateSettings(opt)

if type(opt.number) == type([]):
    # Number is a vector - use it without any modification
    nbox = opt.number
else:
    if opt.cubicdomain:
        nbox = FindBestCube(opt)
    else:
        nbox = FindBestCuboid(opt)
        
# Autorotate
if opt.autorotate:
    os.system('molecuilder --parse-tremolo-potentials %s -i rotated_temp_source.data -l %s --rotate-to-principal-axis-system "1, 0, 0"' % (opt.potentialsfiledir+opt.basename+'.potentials', opt.source))
    opt.source = 'rotated_temp_source.data'

avogadro =  6.022143e23
VolumePerMolecule = opt.molarmass / (avogadro * opt.density)
cell = [VolumePerMolecule**(1./3)] * 3

if not opt.cubiccell:
    try:
        bb = GetSourceBBabs(opt)
        print '======== BBOX:', bb
        # Scaling factor - the molecules bounding box is scaled to fit the volume suiting the density
        s = (VolumePerMolecule / (bb[0]*bb[1]*bb[2])) ** (1./3)

        if s < 1:
            print 'Warning: Molecular cells will overlap.'

        for i in range(3):
            cell[i] = bb[i]*s
    except ZeroDivisionError:
        print 'Warning:  Singularity in bounding box, falling back to cubic cell.'
        

print '======== CELL: ', cell

import pyMoleCuilder as mol
mol.CommandVerbose('0')
mol.ParserParseTremoloPotentials(opt.potentialsfiledir + opt.basename + '.potentials')
mol.WorldInput(opt.source)
mol.WorldCenterInBox('%f 0 0 %f 0 %f' % tuple(cell))
mol.WorldRepeatBox('%d %d %d' % tuple(nbox))
mol.WorldOutput(opt.outfilename + '.data')
mol.WorldOutput(opt.outfilename + '.xyz')

domain = [0.0]*3

for i in range(3):
    domain[i] = cell[i]*nbox[i]

print  '======== DOMAIN: ', domain

# Postprocessing

if opt.postprocess:
    with open(opt.outfilename + '.data') as f:
        ofile = f.read()
    
    with open(opt.outfilename + '.data', 'w') as f:
        f.write('# INPUTCONV shift center\n')
        
        if opt.temp:
            f.write('# INPUTCONV temp %.4f\n' % opt.temp)
            
        f.write(ofile)
        
if opt.autorotate:
    os.system('rm ' + opt.source)