search

JuliaMolSim / Molly.jl

Molecular simulation in Julia
Other
393 stars 53 forks source link

Boltzmann Constant in LAMMPS 'real' units does not work #127

Closed ejmeitz closed 1 year ago

ejmeitz commented 1 year ago

I've been trying to make a system using the units below. In LAMMPS these are called "real" units and are one of the most commonly used unit systems there. I'm not the biggest fan of it but a lot existing data and constants are already in these units. When I try to make a system from these units the convert_k() function does not like the units on the k I've passed even though they are the correct units for this system. The only difference between this and the default is kJ --> kcal but I think the problem is that the k I passed in is per mol.

If I leave it as the Unitful.k then it converts to the incorrect units. Why do we convert the units on k after the user passes it in? Is it in case they change the energy units but do not change k?

energy_units = u"kcal * mol^-1",
force_units = u"kcal * mol^-1 * Å^-1",
k = 1.987e-3u"kcal * mol^-1 * K^-1"
jgreener64 commented 1 year ago

Can you paste code giving the error? As I recall we convert k to make kT the same units as energy_units, with any molar term removed.

ejmeitz commented 1 year ago

Will do later, if you construct a System with those arguments it should give the error though.

jgreener64 commented 1 year ago

I get an error with all 3 arguments, but it works without specifying k. In this case sys.k ends up as 3.299830305927343e-27 kcal K^-1 which is equal to k / Unitful.Na.

Alternatively passing in k=3.299491e-27u"kcal * K^-1" also works.

ejmeitz commented 1 year ago

Yeah I got the same results, I was not expecting k to be converted to away from being molar when the energy I specify is molar. At least for my use case kB needed to be molar for the math to work out.

Internally where is kB being used that the unit conversion is important? I didn't look that closely.

Also is there a reason to not support passing kB as molar?

jgreener64 commented 1 year ago

Internally where is kB being used that the unit conversion is important? I didn't look that closely.

Off the top of my head it is things like velocity generation and temperature. There might be a more principled way to do it.

More generally k is only supposed to be supplied to the constructor if it is numerically different to the real Boltzmann constant, which is apparently a useful feature in some cases. In most cases, including here, you can just omit k and it gets converted appropriately.