Inconsistent Units Crash Simulation

[ejmeitz]

When writing code for SHAKE/RATTLE I ran into a weird bug that crashed my simulation. I pasted a MWE below. When I run this code I get the error that "system force units are kcal Å^-1 mol^-1 but encountered force units kcal nm^-1 mol^-1". If I change box_size to be defined in Angstroms or do not define \sigma in the Atom constructor the error goes away. I'm not really sure what's going on here since Unitful should be handling this conversion, might be something wrong my Julia. If someone else could try running this that would be great. Thanks!

Edit: Realized might just be an issue with the check_energy_units& check_force_unitsfunctions as they do not convert units.

e.g. unit(1.0u"nm") == u"Å" will give you false. So the math might be getting done right and this function just throws a random error.

r_cut = 8.5u"Å"
temp = 10.0u"K"

n_atoms = 20
atom_mass = 10.0u"u"
atoms = [Atom(mass=atom_mass, σ=3.4u"Å", ϵ=0.24u"kcal* mol^-1") for i in 1:n_atoms]
box_size = 2.0u"nm" # nm
coords = [box_size .* rand(SVector{3}) for i in 1:n_atoms]

velocities = [random_velocity(atom_mass, temp) for i in 1:length(atoms)]
boundary = CubicBoundary(box_size)

sys = System(
        atoms = atoms,
        coords = coords,
        boundary = boundary,
        velocities=velocities,
        pairwise_inters=(
            LennardJones(
                energy_units = u"kcal * mol^-1",
                force_units = u"kcal * mol^-1 * Å^-1"
                ),
            ),
        energy_units = u"kcal * mol^-1",
        force_units = u"kcal * mol^-1 * Å^-1"
        )

simulator = VelocityVerlet(dt = 2.0u"fs")
simulate!(sys, simulator, 100, n_threads = 1)

[bondrewd]

Hi! I tried your code and got the same error. As you pointed out, it relates to how Unitful treats units when doing basic arithmetic. In your code, the coordinates are defined based on the boundary units. Setting the boundary in u"\AA" makes the coordinates be defined in u"\AA", and the same applies to u"nm". However, after calculating the interaction in the Lennard-Jones potential, you end up with mixed units if you are not careful! For example, if you try Unitful in the Julia REPL, you will find that:

julia> 1.0u"nm" / 10.0u"Å"
0.1 nm Å⁻¹

julia> 1.0u"nm" / 1.0u"nm"
1.0

So the units only cancel out if they are the same. It made me scratch my head for some hours.

Regarding not defining \sigma, the interaction follows a shortcut when either \epsilon or \sigma is equal to zero (and returns 0 with the right units):

https://github.com/JuliaMolSim/Molly.jl/blob/38dec32b1668bbe2275e2c23edf1bf4206b6636f/src/interactions/lennard_jones.jl#L86-L89

And \sigma is set equal to zero by default when you don't specify it in the constructor:

https://github.com/JuliaMolSim/Molly.jl/blob/38dec32b1668bbe2275e2c23edf1bf4206b6636f/src/types.jl#L235-L243

This triggers the shortcut, so you don't see the error anymore. I hope this helps! :)

[jgreener64]

Yes it seems like check_energy_units and check_force_units are too strict and should be re-written to allow a unit conversion. I can take a look later this week.

[jgreener64]

It appears to be non-trivial to allow the unit conversion in a performant way, especially on the GPU.

One option is to document (in both the docs and the error message) that units returned from force/energy functions should be the same as the system force/energy units, not just valid for conversion. In practice this means that the box size needs to be specified in the same units as appear in the system force units, since the box size is used to get the vector between atoms.

Another option is to do some unit conversion during setup - I am less keen on this.

[ejmeitz]

Simplest thing is probably to check that the units match in setup and then throw an error if they dont.

[ejmeitz]

Just adding this cause its relevant to this issue, but changing mass to g/mol instead of u also causes issues.

For the code above it specifically breaks in the velocities of the system. I think this is because the random_velocity function assumes the mass units will be u with the default params.