mosdef-hub / gmso

Flexible storage of chemical topology for molecular simulation
https://gmso.mosdef.org
MIT License
52 stars 33 forks source link

force field units #187

Open ahy3nz opened 4 years ago

ahy3nz commented 4 years ago

When defining parameters in a force field file, specify dimensions (distance, energy, mass, etc.) rather than units (nm, kj/mol, amu), but then infer units from the dimensions given some specification of a unit system {'energy': kJ/mol, 'distance':nm}

umesh-timalsina commented 4 years ago

This can be done like this:

def _parse_units(unit_tag):
    import unyt as u
    if unit_tag is None:
        unit_tag = {}
    units_map = {
        'energy': u.kcal / u.mol,
        'distance': u.nm,
        'mass': u.amu,
        'charge': u.coulomb
    }
    for attrib, val in unit_tag.items():
        try:  # First Pass: Direct Units
            units_map[attrib] = UNITS_MAP[val]
        except KeyError:  # Second Pass: Check for '/' or ^
            units = val.replace('^', '**')
            exponent = 1
            complex_unit = []
            for unit in units.split('/'):
                actual_unit = unit.strip()
                if '**' in actual_unit:
                    exponent = int(actual_unit.split('**')[-1].strip())
                    actual_unit = actual_unit.split('**')[0].strip()
                try:
                    complex_unit.append(UNITS_MAP[actual_unit])
                    if exponent != 1:
                        complex_unit[-1] = complex_unit[-1] ** exponent
                    this_attrib_unit = complex_unit[0]
                    for _unit in complex_unit[1:]:
                        this_attrib_unit = this_attrib_unit / _unit
                    units_map[attrib] = this_attrib_unit
                except KeyError:
                    warnings.warn('No Unit for {0} named {1} exists, going with default {2}'
                                  .format(attrib, val, units_map[attrib]))
    return units_map

Where UNITS_MAP happens to be like this:

def _form_units_dict():
    import unyt as u
    unyts_dict = {}
    for key, value in u.__dict__.items():
        if isinstance(value, (u.Unit, u.unyt_quantity)):
            unyts_dict[key] = value
    del u
    return unyts_dict
UNITS_MAP = _form_units_dict()
umesh-timalsina commented 4 years ago

For a more complicated unit a regex based splitter that converts the infix operations to postfix might be necessary. For example (Miles/sec)**2

umesh-timalsina commented 4 years ago

This might sound so foolish, and I wasted nearly 4 hours for this, but combining is a one liner solution the Unyts module constructor can take in very complicated expressions as well.

In [65]: import unyt                                                                                                                                                                               

In [66]: c = unyt.Unit('kcal/mol*nm**6')                                                                                                                                                           

In [67]: b = unyt.Unit('cal/mol*nm**6')                                                                                                                                                            

In [68]: b*1 +  c * 2                                                                                                                                                                              
Out[68]: unyt_quantity(2001., 'cal*nm**6/mol')