Units and sanity checking for box potential with non-periodic b.c.

[clonker]

Units

Introduced a runtime dependency to pint. Units can now be used in all relevant api functions:

---

Construct a system with a simulation box size of (10, 10, 10) and units of:
- length in nanometers
- time in nanoseconds
- energy in kJ/mol
The temperature is defaulted to 2.437 kJ/mol, all boundaries are set to be periodic.

import readdy
system = readdy.ReactionDiffusionSystem([10., 10., 10.])

The simulation box size is the only required constructor argument of a reaction diffusion system.

---

Construct a system with a simulation box size of (10, 10, 10) and units of:
- length in kilometers
- time in hours
- energy in kcal/mol

import readdy
system = readdy.ReactionDiffusionSystem([10, 10, 10] * readdy.units.meters,
            length_unit='kilometer', time_unit='hour', energy_unit='kilocal/mol')
print(system.box_size)
>>> [ 0.01  0.01  0.01] kilometer
print(system.kbt)
>>> 0.5824569789674953 kilocalorie / mole

The room temperature as well as any other units of time/length/energy will be automatically converted to this particular set of system units.

---

Construct a unitless system with a simulation box size of (10, 10, 10):

import readdy
system = readdy.ReactionDiffusionSystem([10, 10, 10], length_unit=None, time_unit=None,
            energy_unit=None)
print(system.box_size)
>>> [ 10.  10.  10.]
print(system.kbt)
>>> 2.437

If time, length, and energy units are set to None or empty strings, the system will be treated as completely unitless system.

---

If properties are set without providing a particular unit, e.g.,

system.potentials.add_harmonic_repulsion("A", "B", force_constant=10., 
    interaction_distance=5.)

the system's units are assumed without further conversion, i.e., energy/length**2 for the force constant, and length for the interaction distance. Providing a unit will trigger a conversion:

from readdy import units as units
sys = readdy.ReactionDiffusionSystem(box_size=[1., 1., 1.])
sys.add_species("A", 1.)
sys.add_species("B", 1.)
sys.potentials.add_harmonic_repulsion("A", "B", 
    force_constant=10. * units.kilocal / (units.mol * units.mile**2), 
    interaction_distance=5.)

or also trigger an error if the dimensionality does not match:

sys.potentials.add_harmonic_repulsion("A", "B", 
    force_constant=10. * units.kilocal / (units.mol * units.mile**3), 
    interaction_distance=5.)

raises an

DimensionalityError: Cannot convert from 'kilocalorie / mile ** 3 / mole' ([mass] / [length] / [substance] / [time] ** 2) to 'kilojoule / mole / nanometer ** 2' ([mass] / [substance] / [time] ** 2)

PBC and box potentials

Upon simulation start there will be a check on the particle types together with periodicity of the simulation box and box potentials, in particular:

If the simulation box is not completely periodic and if the particle types in that box have no box potential to keep them contained in the non-periodic directions, i.e.,

box_lower_left[dim] >= potential_lower_left[dim] 
    or box_upper_right[dim] <= potential_upper_right[dim] 

is true for any non-periodic dim, an exception is thrown.

[clonker]

Added a bond breaking preset for structural topology reactions, example:

sys = readdy.ReactionDiffusionSystem(box_size=[150,150,150])
sys.topologies.add_type("Toptype")
sys.add_topology_species("T", diffusion_constant=1.)
sys.topologies.configure_harmonic_bond("T", "T", force_constant=20., length=2.)
# breaks bonds with a topology-global rate of .01 * n_edges
sys.topologies.add_topology_dissociation("Toptype", .01)

[franknoe]

Temperature is not a Energy. I think you mean that kT is 2.437 kJ/mol, and the temperature would probably be 300 K or so. If you are using units, I suggest that the temperature is the unit to be set, and kT is computed. Without units the convention would be that energies are measured in kT.

[clonker]

I have introduced a temperature unit (default Kelvin) and removed the setter property for kbt in favor of propertiy getter and setter temperature in units of temperature, internally expressed in kT. Also the constructor now takes a unit_system (which is defaulted to 'default') that can be overridden by providing a dictionary. In particular:

# default units
system = readdy.ReactionDiffusionSystem(box_size=[1., 1., 1.])
# temperature in kelvin, internally expressed in kT/mol
system.temperature = 293
print(system.temperature)
>>> 293
print(system.kbt)
>>> 2.4361374086224026 kJ/mol

# custom units
system = ReactionDiffusionSystem(box_size=[1., 1., 1.], unit_system={
    'energy_unit': 'kcal/mol',
    'temperature_unit': 'rankine',
    'time_unit': 'hour'
})
print(system.temperature.to(system.units.kelvin))
>>> 293 kelvin

[franknoe]

OK, but the sentence "internally expressed in kT/mol" does not make sense.

[clonker]

Yeah right, sorry, I meant kJ/mol.

[franknoe]

No, the temperature cannot be expressed in kJ/mol. kT could be expressed in kJ/mol. T is the temperature, kT is the thermal energy.

[franknoe]

BTW, I don't know what "rankine" is

[clonker]

Ok thanks! I have updated the documentation and will merge when the tests are passing if that is alright.

0° Rankine are, same as 0° Kelvin, absolute zero, however Rankine operates on the Fahrenheit scale.