search

JuliaMolSim / Molly.jl

Molecular simulation in Julia
Other
393 stars 53 forks source link

Crystals v2 #125

Closed ejmeitz closed 1 year ago

ejmeitz commented 1 year ago

Don't even know what happened in the last branch. Unbroke whatever merge issues I had from pulling the CUDA branch.

Things to remember:

codecov[bot] commented 1 year ago

Codecov Report

Patch coverage: 88.23% and no project coverage change.

Comparison is base (5cf4d29) 72.82% compared to head (660a66c) 72.82%.

Additional details and impacted files ```diff @@ Coverage Diff @@ ## master #125 +/- ## ======================================= Coverage 72.82% 72.82% ======================================= Files 34 34 Lines 4798 4821 +23 ======================================= + Hits 3494 3511 +17 - Misses 1304 1310 +6 ``` | [Impacted Files](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/125?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim) | Coverage Δ | | |---|---|---| | [src/Molly.jl](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/125?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL01vbGx5Lmps) | `100.00% <ø> (ø)` | | | [src/types.jl](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/125?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim#diff-c3JjL3R5cGVzLmps) | `72.36% <88.23%> (+0.93%)` | :arrow_up: | ... and [3 files with indirect coverage changes](https://app.codecov.io/gh/JuliaMolSim/Molly.jl/pull/125/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=JuliaMolSim)

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

ejmeitz commented 1 year ago

@jgreener64 if you get the chance could you look at the crystal constructor I added to types.jl and help me fix the type instability? I think the instability is because of the Boundary if-else block but I'm not really sure how to fix that as the Boundary is something that is different depending on what crystal gets passed (e.g. cubic for FCC triclinic for others)

Probably can just dispatch on the type somehow but then Ill have like 4 crystal constructors

jgreener64 commented 1 year ago

What do you mean by type instability in this context? Can you give example code that shows it?

ejmeitz commented 1 year ago

Yeah for sure, below is an example of constructing a system from a 4x4x4 unit-cell FCC system. The construction of the System struct is type-unstable and when I call simulate!() with this struct the run-times are much slower than I expect. The results of the simulation are what I expect but the output from @code_warntype on the construction of the System indicates that the compiler cannot assign concrete types to everything. This will break auto-diff break and slow down execution.

using Molly
using SimpleCrystals

fcc_crystal = FCC(5.4u"Å", :Ar, SVector(4,4,4))

n_atoms = length(fcc_crystal)
atom_mass = atomic_mass(fcc_crystal,1)
velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms]

sys = @code_warntype System(
    fcc_crystal,
    velocities=velocities,
    pairwise_inters=(LennardJones(),),
    loggers=(ke=KineticEnergyLogger(10),
            pe=PotentialEnergyLogger(10)),
    energy_units = u"kJ * mol^-1",
    force_units = u"kJ * mol^-1 * Å^-1")

If you actually want this code to run with LJ you have to use that extra convenience constructor and modify all the atoms in the system to have the proper sigma and epsilon values.

σ = 3.4u"Å" 
ϵ = (4.184*0.24037)u"kJ * mol^-1"
updated_atoms = []
for i in range(1,length(sys.atoms))
    push!(updated_atoms, Molly.Atom(index = sys.atoms[i].index, charge = sys.atoms[i].charge,
     mass = sys.atoms[i].mass, σ = σ, ϵ = ϵ, solute = sys.atoms[i].solute))
end

sys = System(sys, atoms = updated_atoms)
jgreener64 commented 1 year ago

The constructor itself can be type unstable since it is only called once at the start. The concern is what is returned, i.e. typeof(sys). In the first case this is concrete so all is well - here this can be checked by seeing that all the type parameters in typeof(sys) are concrete.

The problem is that updated_atoms = [] creates a Vector{Any} which is not concrete. A quick fix is sys = System(sys, atoms = [updated_atoms...]).

ejmeitz commented 1 year ago

Cool yeah

[updated_atoms...]

fixed the call to simulate! which is what matters. I guess I'm unhappy now that the user has to avoid that pitfall or will suffer 50% performance loss.

The interface works from the few systems I've ran but I still am not satisfied with the extra work/knowledge put onto the user with the current implementation. The code 2 comments above is basically what you will have to do every time. I'm open to suggestions of how I could change SimpleCrystals or the code I've added to Molly. For now my only idea requires breaking changes to Molly (moving $\sigma$ and $\epsilon$ and similar params to potential definitions).

jgreener64 commented 1 year ago

σ and ϵ only need to be set in cases where they will be used, i.e. there is a Lennard Jones or similar interaction. I am less familiar with these systems - how often would a Lennard Jones-type interaction be used when simulating? Always, sometimes or never?

Moving the parameters to the interactions is probably the right way to go long term but requires a lot of changes and careful thought about how to not make setup harder in other contexts.

ejmeitz commented 1 year ago

I use LJ to test that the system set up is correct. As far as giving accurate properties LJ is really only good for Argon, but I know lots of people use it to study phenomenological things because it is the cheapest potential to evaluate.

However useful LJ is, other potentials like Buckingham require atom properties (A_ij, B_ij) as well which are not even a part of the base Atom struct. I'd be willing to undertake this change as I think long term it is the correct thing, but it might require changing how interactions are defined (I have to look more into this). Right now mixing rules are automatically calculated I know we have InteractionList2Atoms etc. which I've never used.

I think there should be limited breaking changes but lots of internal changes.

jgreener64 commented 1 year ago

I'm not going to stop you giving it a go, but I would warn you that it is really a lot of work and probably not worth it just to make the above setup process easier. There are considerations about differentiable simulations and GPU performance, for example. It might be better to just document the setup procedure clearly.

The lack of properties for Buckingham potentials etc. is an issue but it is possible to use a custom atom type without a performance penalty, one is used in the tests for Buckingham:https://github.com/JuliaMolSim/Molly.jl/blob/2c0631401a2b51d5215308fc1d877456b8f60411/test/interactions.jl#L80-L85

ejmeitz commented 1 year ago

Sounds good, I guess for now I'll just write documentation and add some tests.

ejmeitz commented 1 year ago

Tests should pass, in conclusion this PR adds:

Writing the tests for this made me realize that my NoseHoover implementation doesn't work for solids, made another PR to add that to the docs. Might need help to fix that.