Inconsistent Units Crash Simulation about molly.jl

[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):

Molly.jl/src/interactions/lennard_jones.jl

Lines 86 to 89 in 38dec32

	if !S && (iszero_value(atom_i.ϵ) || iszero_value(atom_j.ϵ) ||
	iszero_value(atom_i.σ) || iszero_value(atom_j.σ))
	return ustrip.(zero(coord_i)) * inter.force_units
	end

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

Molly.jl/src/types.jl

Lines 235 to 243 in 38dec32

	function Atom(;
	index=1,
	charge=0.0,
	mass=1.0u"u",
	σ=0.0u"nm",
	ϵ=0.0u"kJ * mol^-1",
	solute=false)
	return Atom(index, charge, mass, σ, ϵ, solute)
	end

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

[ejmeitz]

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

[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]

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.