juliamolsim / molly.jl Goto Github PK

Hey,
first and foremost, thank you for the great package!
As disclaimer, I am super new to the whole MD thing, so this might be totally on my end :>

I was trying to load a .pdb file of an Alanine-Dipeptide but got the error message that no "H" key was found.
Digging deeper I then found out that this was due to the default rename_terminal_res=true in System.
This replaced the ALA with NALA, but my .pdb only supplied an H atom, not the H1, H2, H3 required in the force field for NALA.

If I understand correctly this is a feature, not a bug. However with only the error message that no key was found I was pretty lost to start with. I cannot judge how sane the default rename_terminal_res setting is, but maybe there should be some kind of check whether ALA -> NALA is actually possible, or kind of a better error message?

It's fine for me now (just setting rename_terminal_res=false), but I thought I'd leave this here at least for others stepping into the same trap :)

I updated CellListMap v0.7.21 ⇒ v0.8.3 and received the following errors on testing the master branch.

Neighbor lists: Error During Test at /home/leios/projects/CESMIX/Molly.jl/test/basic.jl:117
  Got exception outside of a @test
  ArgumentError: zero(Quantity{Float64}) not defined.
  Stacktrace:
    [1] zero(x::Type{Quantity{Float64}})
      @ Unitful ~/.julia/packages/Unitful/ApCuY/src/quantities.jl:389
    [2] (::CellListMap.var"#95#96"{SMatrix{3, 3, Quantity{Float64}, 9}})(el::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}})
      @ CellListMap ~/.julia/packages/CellListMap/rqPY5/src/CellOperations.jl:547
    [3] (::StaticArrays.var"#194#195"{CellListMap.var"#95#96"{SMatrix{3, 3, Quantity{Float64}, 9}}})(x::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}})
      @ StaticArrays ~/.julia/packages/StaticArrays/atiPw/src/mapreduce.jl:259
    [4] macro expansion
      @ ~/.julia/packages/StaticArrays/atiPw/src/mapreduce.jl:140 [inlined]
    [5] _mapfoldl
      @ ~/.julia/packages/StaticArrays/atiPw/src/mapreduce.jl:115 [inlined]
    [6] _mapreduce
      @ ~/.julia/packages/StaticArrays/atiPw/src/mapreduce.jl:113 [inlined]
    [7] #count#193
      @ ~/.julia/packages/StaticArrays/atiPw/src/mapreduce.jl:259 [inlined]
    [8] count
      @ ~/.julia/packages/StaticArrays/atiPw/src/mapreduce.jl:259 [inlined]
    [9] check_unit_cell(unit_cell_matrix::SMatrix{3, 3, Quantity{Float64}, 9}, cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}; printerr::Bool)
      @ CellListMap ~/.julia/packages/CellListMap/rqPY5/src/CellOperations.jl:547
   [10] check_unit_cell(unit_cell_matrix::SMatrix{3, 3, Quantity{Float64}, 9}, cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}})
      @ CellListMap ~/.julia/packages/CellListMap/rqPY5/src/CellOperations.jl:537
   [11] CellListMap.Box(unit_cell_matrix::SMatrix{3, 3, Quantity{Float64}, 9}, cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, lcell::Int64, #unused#::Type{CellListMap.OrthorhombicCell})
      @ CellListMap ~/.julia/packages/CellListMap/rqPY5/src/Box.jl:167
   [12] Box
      @ ~/.julia/packages/CellListMap/rqPY5/src/Box.jl:268 [inlined]
   [13] #Box#14
      @ ~/.julia/packages/CellListMap/rqPY5/src/Box.jl:270 [inlined]
   [14] CellListMap.Box(sides::SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}, cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}})
      @ CellListMap ~/.julia/packages/CellListMap/rqPY5/src/Box.jl:270
   [15] CellListMapNeighborFinder(; nb_matrix::BitMatrix, matrix_14::BitMatrix, n_steps::Int64, x0::Nothing, unit_cell::Nothing, number_of_batches::Tuple{Int64, Int64}, dist_cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}})
      @ Molly ~/projects/CESMIX/Molly.jl/src/neighbors.jl:316
   [16] macro expansion
      @ ~/projects/CESMIX/Molly.jl/test/basic.jl:119 [inlined]
   [17] macro expansion
      @ ~/builds/julia-1.7.1/share/julia/stdlib/v1.7/Test/src/Test.jl:1283 [inlined]
   [18] top-level scope
      @ ~/projects/CESMIX/Molly.jl/test/basic.jl:118
   [19] include(fname::String)
      @ Base.MainInclude ./client.jl:451
   [20] top-level scope
      @ ~/projects/CESMIX/Molly.jl/test/runtests.jl:70
   [21] include(fname::String)
      @ Base.MainInclude ./client.jl:451
   [22] top-level scope
      @ none:6
   [23] eval
      @ ./boot.jl:373 [inlined]
   [24] exec_options(opts::Base.JLOptions)
      @ Base ./client.jl:268
   [25] _start()
      @ Base ./client.jl:495
Test Summary:  | Pass  Error  Total
Neighbor lists |    4      1      5
ERROR: LoadError: Some tests did not pass: 4 passed, 0 failed, 1 errored, 0 broken.
in expression starting at /home/leios/projects/CESMIX/Molly.jl/test/basic.jl:117
in expression starting at /home/leios/projects/CESMIX/Molly.jl/test/runtests.jl:61
ERROR: Package Molly errored during testing

It seems like the zero function in Unitful will purposefully throw an error to a Quantity Type without a unit?

zero(x::AbstractQuantity) = Quantity(zero(x.val), unit(x))
zero(x::AffineQuantity) = Quantity(zero(x.val), absoluteunit(x))
zero(x::Type{<:AbstractQuantity{T}}) where {T} = throw(ArgumentError("zero($x) not defined."))
zero(x::Type{<:AbstractQuantity{T,D}}) where {T,D} = zero(T) * upreferred(D)
zero(x::Type{<:AbstractQuantity{T,D,U}}) where {T,D,U<:ScalarUnits} = zero(T)*U()
zero(x::Type{<:AbstractQuantity{T,D,U}}) where {T,D,U<:AffineUnits} = zero(T)*absoluteunit(U())

The zero function is used in the following block for CellListMaps:

function check_unit_cell(unit_cell_matrix::SMatrix{3},cutoff;printerr=true)
    ...
    if count(el -> el < zero(eltype(unit_cell_matrix)), unit_cell_matrix) != 0
         printerr && println("UNIT CELL CHECK FAILED: unit cell matrix components be strictly positive.")
         check = false
    end
    ...
end

I've gotten this error on two different machines, but I might still be doing something wrong locally. Here is my ] st:

     Project Molly v0.13.0
      Status `~/projects/CESMIX/Molly.jl/Project.toml`
  [a963bdd2] AtomsBase v0.2.3
  [de9282ab] BioStructures v1.2.1
  [052768ef] CUDA v3.12.0
  [69e1c6dd] CellListMap v0.8.3
  [d360d2e6] ChainRulesCore v1.15.5
  [46823bd8] Chemfiles v0.10.3
  [5ae59095] Colors v0.12.8
  [861a8166] Combinatorics v1.0.2
  [864edb3b] DataStructures v0.18.13
  [b4f34e82] Distances v0.10.7
  [31c24e10] Distributions v0.25.73
  [8f5d6c58] EzXML v1.1.0
  [cc61a311] FLoops v0.2.0
  [f6369f11] ForwardDiff v0.10.32
  [5ab0869b] KernelDensity v0.6.5
  [b8a86587] NearestNeighbors v0.4.11
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.0
  [90137ffa] StaticArrays v1.5.7
  [1986cc42] Unitful v1.12.0
  [f31437dd] UnitfulChainRules v0.1.2
  [e88e6eb3] Zygote v0.6.44
  [37e2e46d] LinearAlgebra
  [9a3f8284] Random
  [2f01184e] SparseArrays
  [10745b16] Statistics

note that the tests do pass on [email protected]

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"

Hello,

I'm currently working on improving the interoperability between InteratomicPotentials.jl and Molly. I have a mostly working prototype (based off of some past work) using Lennard-Jones as an example that can be seen here. The near-term goal is to get various machine-learned interatomic potentials working as well.

The primary issue is that the forces outputted from InteratomicPotentials.jl are in atomic units (i.e. Hartree/Bohr), but I typically use units of eV/Å. However, if I set my ForceUnits=u"eV/Å" in my Molly system, it appears that no unit checking occurs (for the general interactions), the units are stripped, and the user-specified units are incorrectly appended, with no unit conversion. While I can manage the units on my end, I think it would be safer/more reliable if these units were converted appropriately within Molly.

I've discussed this issue a bit with @leios already, but any additional help is appreciated! Thank you!

There seems to be some changes in Zygote.jl in the latest release in which the function Zygote._backmean was removed. Molly provides a redefinition of this function in

Due to this the package fails to precompile.

See the commit: FluxML/Zygote.jl@af434d6#diff-a9e025ac90a30d27e7512546971c5d92ea7c3496ba759336ae6bf1cace6db4b2

If we use some value for m and n such that (n / (n - m)) * (n / m) ^ (m / (n - m)) gives a floating point number then,

m, n, mn_fac = promote(m, n, (n / (n - m)) * (n / m) ^ (m / (n - m)))

This might be because of the re-direct thing I merged last week.

There is a special case when the calculation of angle force breaks down, namely when the cross product of the bond vectors vanishes. This may occur in initial configurations with completely straight polymer chains. The behavior is somewhat random due to floating point errors, making the identification of the cause for the NaN values a bit tricky for the user.

Perhaps the algorithm could be modified to work in all scenarios?

Line #1 of force.jl:

# See https://udel.edu/~arthij/MD.pdf for information on forces

contains a dead link.

As the title says, I think Molly needs a cool logo. Most Julia packages have made their logos using Luxor.jl maybe one can be made for Molly as well (although I have no experience with Luxor).

in expression starting at /home/ian/.julia/packages/Molly/7GX09/test/runtests.jl:29
ERROR: Package Molly errored during testing

_ _ ()_ | Documentation: https://docs.julialang.org
() | () () |
_ _ | | __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ ` | |
| | || | | | (| | | Version 1.7.3 (2022-05-06)
/ |_'|||_'_| | Official https://julialang.org/ release
|__/ |

Installed ChainRules ────────── v1.39.0
Installed MLStyle ───────────── v0.4.13
Installed BufferedStreams ───── v1.1.0
Installed PeriodicTable ─────── v1.1.2
Installed UnitfulChainRules ─── v0.1.0
Installed Twiddle ───────────── v1.1.2
Installed BioAlignments ─────── v2.0.0
Installed PooledArrays ──────── v1.4.2
Installed Format ────────────── v1.3.2
Installed KernelDensity ─────── v0.6.4
Installed FLoops ────────────── v0.2.0
Installed IntervalTrees ─────── v1.0.0
Installed PrettyTables ──────── v1.3.1
Installed JuliaVariables ────── v0.2.4
Installed YAML ──────────────── v0.4.7
Installed Chemfiles_jll ─────── v0.10.2+0
Installed StableRNGs ────────── v0.1.2
Installed BioStructures ─────── v1.2.1
Installed UnitfulAtomic ─────── v1.0.0
Installed CellListMap ───────── v0.7.17
Installed BioSequences ──────── v2.0.5
Installed Molly ─────────────── v0.12.0
Downloaded artifact: Chemfiles
Updating ~/.julia/environments/v1.7/Project.toml
[aa0f7f06] + Molly v0.12.0
Updating ~/.julia/environments/v1.7/Manifest.toml
[a963bdd2] + AtomsBase v0.2.2
[67c07d97] + Automa v0.8.2
[b99e7846] + BinaryProvider v0.5.10
[00701ae9] + BioAlignments v2.0.0
[37cfa864] + BioCore v2.0.5
[47718e42] + BioGenerics v0.1.1
[7e6ae17a] + BioSequences v2.0.5
[de9282ab] + BioStructures v1.2.1
[3c28c6f8] + BioSymbols v4.0.4
[e1450e63] + BufferedStreams v1.1.0
[69e1c6dd] + CellListMap v0.7.17
[082447d4] + ChainRules v1.39.0
[46823bd8] + Chemfiles v0.10.2
[6add18c4] + ContextVariablesX v0.1.2
[a93c6f00] + DataFrames v1.3.4
[8f5d6c58] + EzXML v1.1.0
[cc61a311] + FLoops v0.2.0
[b9860ae5] + FLoopsBase v0.1.1
[1fa38f19] + Format v1.3.2
[7869d1d1] + IRTools v0.4.6
[1cb3b9ac] + IndexableBitVectors v1.0.0
[524e6230] + IntervalTrees v1.0.0
[41ab1584] + InvertedIndices v1.1.0
[b14d175d] + JuliaVariables v0.2.4
[5ab0869b] + KernelDensity v0.6.4
[259c3a9c] + MMTF v1.0.0
[626554b9] + MetaGraphs v0.7.1
[aa0f7f06] + Molly v0.12.0
[99f44e22] + MsgPack v1.1.0
[71a1bf82] + NameResolution v0.1.5
[7b2266bf] + PeriodicTable v1.1.2
[2dfb63ee] + PooledArrays v1.4.2
[8162dcfd] + PrettyPrint v0.2.0
[08abe8d2] + PrettyTables v1.3.1
[92933f4c] + ProgressMeter v1.7.2
[c1ae055f] + RealDot v0.1.0
[fdea26ae] + SIMD v3.4.1
[7b38b023] + ScanByte v0.3.3
[860ef19b] + StableRNGs v0.1.2
[69024149] + StringEncodings v0.3.5
[7200193e] + Twiddle v1.1.2
[a7773ee8] + UnitfulAtomic v1.0.0
[f31437dd] + UnitfulChainRules v0.1.0
[ddb6d928] + YAML v0.4.7
[e88e6eb3] + Zygote v0.6.40
[78a364fa] + Chemfiles_jll v0.10.2+0
[9abbd945] + Profile
Building Chemfiles → ~/.julia/scratchspaces/44cfe95a-1eb2-52ea-b672-e2afdf69b78f/a1d1f4b68126ebead888d4d9119a58a80516df63/build.log
Precompiling project...
47 dependencies successfully precompiled in 342 seconds (328 already precompiled)

(@v1.7) pkg> test Molly
Testing Molly
Status /tmp/jl_mVMQ6u/Project.toml
[4c88cf16] Aqua v0.5.5
[a963bdd2] AtomsBase v0.2.2
[de9282ab] BioStructures v1.2.1
[052768ef] CUDA v3.11.0
[69e1c6dd] CellListMap v0.7.17
[d360d2e6] ChainRulesCore v1.15.0
[46823bd8] Chemfiles v0.10.2
[5ae59095] Colors v0.12.8
[861a8166] Combinatorics v1.0.2
[864edb3b] DataStructures v0.18.13
[b4f34e82] Distances v0.10.7
[31c24e10] Distributions v0.25.62
[8f5d6c58] EzXML v1.1.0
[cc61a311] FLoops v0.2.0
[26cc04aa] FiniteDifferences v0.12.24
[f6369f11] ForwardDiff v0.10.30
[e9467ef8] GLMakie v0.6.12
[5ab0869b] KernelDensity v0.6.4
[aa0f7f06] Molly v0.12.0
[b8a86587] NearestNeighbors v0.4.11
[189a3867] Reexport v1.2.2
[ae029012] Requires v1.3.0
[90137ffa] StaticArrays v1.4.7
[1986cc42] Unitful v1.11.0
[f31437dd] UnitfulChainRules v0.1.0
[e88e6eb3] Zygote v0.6.40
[8bb1440f] DelimitedFiles @stdlib/DelimitedFiles
[37e2e46d] LinearAlgebra @stdlib/LinearAlgebra
[9a3f8284] Random @stdlib/Random
[2f01184e] SparseArrays @stdlib/SparseArrays
[10745b16] Statistics @stdlib/Statistics
[8dfed614] Test @stdlib/Test
Status /tmp/jl_mVMQ6u/Manifest.toml
[621f4979] AbstractFFTs v1.1.0
[1520ce14] AbstractTrees v0.4.2
[79e6a3ab] Adapt v3.3.3
[27a7e980] Animations v0.4.1
[4c88cf16] Aqua v0.5.5
[dce04be8] ArgCheck v2.3.0
[ec485272] ArnoldiMethod v0.2.0
[a963bdd2] AtomsBase v0.2.2
[67c07d97] Automa v0.8.2
[13072b0f] AxisAlgorithms v1.0.1
[ab4f0b2a] BFloat16s v0.2.0
[198e06fe] BangBang v0.3.36
[9718e550] Baselet v0.1.1
[b99e7846] BinaryProvider v0.5.10
[00701ae9] BioAlignments v2.0.0
[37cfa864] BioCore v2.0.5
[47718e42] BioGenerics v0.1.1
[7e6ae17a] BioSequences v2.0.5
[de9282ab] BioStructures v1.2.1
[3c28c6f8] BioSymbols v4.0.4
[e1450e63] BufferedStreams v1.1.0
[fa961155] CEnum v0.4.2
[052768ef] CUDA v3.11.0
[49dc2e85] Calculus v0.5.1
[69e1c6dd] CellListMap v0.7.17
[082447d4] ChainRules v1.39.0
[d360d2e6] ChainRulesCore v1.15.0
[9e997f8a] ChangesOfVariables v0.1.3
[46823bd8] Chemfiles v0.10.2
[944b1d66] CodecZlib v0.7.0
[a2cac450] ColorBrewer v0.4.0
[35d6a980] ColorSchemes v3.19.0
[3da002f7] ColorTypes v0.11.4
[c3611d14] ColorVectorSpace v0.9.9
[5ae59095] Colors v0.12.8
[861a8166] Combinatorics v1.0.2
[bbf7d656] CommonSubexpressions v0.3.0
[34da2185] Compat v3.45.0
[a33af91c] CompositionsBase v0.1.1
[187b0558] ConstructionBase v1.3.0
[6add18c4] ContextVariablesX v0.1.2
[d38c429a] Contour v0.6.2
[a8cc5b0e] Crayons v4.1.1
[9a962f9c] DataAPI v1.10.0
[a93c6f00] DataFrames v1.3.4
[864edb3b] DataStructures v0.18.13
[e2d170a0] DataValueInterfaces v1.0.0
[244e2a9f] DefineSingletons v0.1.2
[b429d917] DensityInterface v0.4.0
[163ba53b] DiffResults v1.0.3
[b552c78f] DiffRules v1.11.0
[b4f34e82] Distances v0.10.7
[31c24e10] Distributions v0.25.62
[ffbed154] DocStringExtensions v0.8.6
[fa6b7ba4] DualNumbers v0.6.8
[e2ba6199] ExprTools v0.1.8
[8f5d6c58] EzXML v1.1.0
[c87230d0] FFMPEG v0.4.1
[7a1cc6ca] FFTW v1.5.0
[cc61a311] FLoops v0.2.0
[b9860ae5] FLoopsBase v0.1.1
[5789e2e9] FileIO v1.14.0
[1a297f60] FillArrays v0.12.8
[26cc04aa] FiniteDifferences v0.12.24
[53c48c17] FixedPointNumbers v0.8.4
[1fa38f19] Format v1.3.2
[59287772] Formatting v0.4.2
[f6369f11] ForwardDiff v0.10.30
[b38be410] FreeType v4.0.0
[663a7486] FreeTypeAbstraction v0.9.9
[f7f18e0c] GLFW v3.4.1
[e9467ef8] GLMakie v0.6.12
[0c68f7d7] GPUArrays v8.4.0
[46192b85] GPUArraysCore v0.1.0
[61eb1bfa] GPUCompiler v0.16.1
[5c1252a2] GeometryBasics v0.4.2
[a2bd30eb] Graphics v1.1.2
[86223c79] Graphs v1.7.1
[3955a311] GridLayoutBase v0.9.0
[42e2da0e] Grisu v1.0.2
[34004b35] HypergeometricFunctions v0.3.10
[7869d1d1] IRTools v0.4.6
[a09fc81d] ImageCore v0.9.4
[82e4d734] ImageIO v0.6.6
[1cb3b9ac] IndexableBitVectors v1.0.0
[9b13fd28] IndirectArrays v1.0.0
[d25df0c9] Inflate v0.1.2
[22cec73e] InitialValues v0.3.1
[a98d9a8b] Interpolations v0.13.6
[8197267c] IntervalSets v0.7.1
[524e6230] IntervalTrees v1.0.0
[3587e190] InverseFunctions v0.1.7
[41ab1584] InvertedIndices v1.1.0
[92d709cd] IrrationalConstants v0.1.1
[f1662d9f] Isoband v0.1.1
[c8e1da08] IterTools v1.4.0
[82899510] IteratorInterfaceExtensions v1.0.0
[033835bb] JLD2 v0.4.22
[692b3bcd] JLLWrappers v1.4.1
[682c06a0] JSON v0.21.3
[b835a17e] JpegTurbo v0.1.1
[b14d175d] JuliaVariables v0.2.4
[5ab0869b] KernelDensity v0.6.4
[929cbde3] LLVM v4.14.0
[b964fa9f] LaTeXStrings v1.3.0
[8cdb02fc] LazyModules v0.3.1
[2ab3a3ac] LogExpFunctions v0.3.15
[d8e11817] MLStyle v0.4.13
[259c3a9c] MMTF v1.0.0
[1914dd2f] MacroTools v0.5.9
[ee78f7c6] Makie v0.17.12
[20f20a25] MakieCore v0.3.6
[dbb5928d] MappedArrays v0.4.1
[7eb4fadd] Match v1.2.0
[0a4f8689] MathTeXEngine v0.4.3
[7269a6da] MeshIO v0.4.10
[626554b9] MetaGraphs v0.7.1
[128add7d] MicroCollections v0.1.2
[e1d29d7a] Missings v1.0.2
[66fc600b] ModernGL v1.1.4
[aa0f7f06] Molly v0.12.0
[e94cdb99] MosaicViews v0.3.3
[99f44e22] MsgPack v1.1.0
[77ba4419] NaNMath v0.3.7
[71a1bf82] NameResolution v0.1.5
[b8a86587] NearestNeighbors v0.4.11
[f09324ee] Netpbm v1.0.2
[510215fc] Observables v0.5.1
[6fe1bfb0] OffsetArrays v1.12.6
[52e1d378] OpenEXR v0.3.2
[bac558e1] OrderedCollections v1.4.1
[90014a1f] PDMats v0.11.13
[f57f5aa1] PNGFiles v0.3.16
[19eb6ba3] Packing v0.4.2
[5432bcbf] PaddedViews v0.5.11
[d96e819e] Parameters v0.12.3
[69de0a69] Parsers v2.3.2
[7b2266bf] PeriodicTable v1.1.2
[eebad327] PkgVersion v0.1.1
[995b91a9] PlotUtils v1.3.0
[647866c9] PolygonOps v0.1.2
[2dfb63ee] PooledArrays v1.4.2
[21216c6a] Preferences v1.3.0
[8162dcfd] PrettyPrint v0.2.0
[08abe8d2] PrettyTables v1.3.1
[92933f4c] ProgressMeter v1.7.2
[4b34888f] QOI v1.0.0
[1fd47b50] QuadGK v2.4.2
[74087812] Random123 v1.5.0
[e6cf234a] RandomNumbers v1.5.3
[c84ed2f1] Ratios v0.4.3
[c1ae055f] RealDot v0.1.0
[3cdcf5f2] RecipesBase v1.2.1
[189a3867] Reexport v1.2.2
[05181044] RelocatableFolders v0.3.0
[ae029012] Requires v1.3.0
[708f8203] Richardson v1.4.0
[79098fc4] Rmath v0.7.0
[fdea26ae] SIMD v3.4.1
[7b38b023] ScanByte v0.3.3
[6c6a2e73] Scratch v1.1.1
[efcf1570] Setfield v0.8.2
[65257c39] ShaderAbstractions v0.2.9
[992d4aef] Showoff v1.0.3
[73760f76] SignedDistanceFields v0.4.0
[699a6c99] SimpleTraits v0.9.4
[45858cf5] Sixel v0.1.2
[a2af1166] SortingAlgorithms v1.0.1
[276daf66] SpecialFunctions v2.1.6
[171d559e] SplittablesBase v0.1.14
[860ef19b] StableRNGs v0.1.2
[cae243ae] StackViews v0.1.1
[90137ffa] StaticArrays v1.4.7
[82ae8749] StatsAPI v1.4.0
[2913bbd2] StatsBase v0.33.17
[4c63d2b9] StatsFuns v1.0.1
[69024149] StringEncodings v0.3.5
[09ab397b] StructArrays v0.6.11
[3783bdb8] TableTraits v1.0.1
[bd369af6] Tables v1.7.0
[62fd8b95] TensorCore v0.1.1
[731e570b] TiffImages v0.6.0
[a759f4b9] TimerOutputs v0.5.20
[3bb67fe8] TranscodingStreams v0.9.6
[28d57a85] Transducers v0.4.73
[7200193e] Twiddle v1.1.2
[3a884ed6] UnPack v1.0.2
[1cfade01] UnicodeFun v0.4.1
[1986cc42] Unitful v1.11.0
[a7773ee8] UnitfulAtomic v1.0.0
[f31437dd] UnitfulChainRules v0.1.0
[efce3f68] WoodburyMatrices v0.5.5
[ddb6d928] YAML v0.4.7
[e88e6eb3] Zygote v0.6.40
[700de1a5] ZygoteRules v0.2.2
[6e34b625] Bzip2_jll v1.0.8+0
[83423d85] Cairo_jll v1.16.1+1
[78a364fa] Chemfiles_jll v0.10.2+0
[5ae413db] EarCut_jll v2.2.3+0
[2e619515] Expat_jll v2.4.8+0
[b22a6f82] FFMPEG_jll v4.4.2+0
[f5851436] FFTW_jll v3.3.10+0
[a3f928ae] Fontconfig_jll v2.13.93+0
[d7e528f0] FreeType2_jll v2.10.4+0
[559328eb] FriBidi_jll v1.0.10+0
[0656b61e] GLFW_jll v3.3.6+0
[78b55507] Gettext_jll v0.21.0+0
[7746bdde] Glib_jll v2.68.3+2
[3b182d85] Graphite2_jll v1.3.14+0
[2e76f6c2] HarfBuzz_jll v2.8.1+1
[905a6f67] Imath_jll v3.1.2+0
[1d5cc7b8] IntelOpenMP_jll v2018.0.3+2
[aacddb02] JpegTurbo_jll v2.1.2+0
[c1c5ebd0] LAME_jll v3.100.1+0
[dad2f222] LLVMExtra_jll v0.0.16+0
[dd4b983a] LZO_jll v2.10.1+0
[e9f186c6] Libffi_jll v3.2.2+1
[d4300ac3] Libgcrypt_jll v1.8.7+0
[7e76a0d4] Libglvnd_jll v1.3.0+3
[7add5ba3] Libgpg_error_jll v1.42.0+0
[94ce4f54] Libiconv_jll v1.16.1+1
[4b2f31a3] Libmount_jll v2.35.0+0
[38a345b3] Libuuid_jll v2.36.0+0
[856f044c] MKL_jll v2022.0.0+0
[e7412a2a] Ogg_jll v1.3.5+1
[18a262bb] OpenEXR_jll v3.1.1+0
[458c3c95] OpenSSL_jll v1.1.17+0
[efe28fd5] OpenSpecFun_jll v0.5.5+0
[91d4177d] Opus_jll v1.3.2+0
[2f80f16e] PCRE_jll v8.44.0+0
[30392449] Pixman_jll v0.40.1+0
[f50d1b31] Rmath_jll v0.3.0+0
[02c8fc9c] XML2_jll v2.9.14+0
[aed1982a] XSLT_jll v1.1.34+0
[4f6342f7] Xorg_libX11_jll v1.6.9+4
[0c0b7dd1] Xorg_libXau_jll v1.0.9+4
[935fb764] Xorg_libXcursor_jll v1.2.0+4
[a3789734] Xorg_libXdmcp_jll v1.1.3+4
[1082639a] Xorg_libXext_jll v1.3.4+4
[d091e8ba] Xorg_libXfixes_jll v5.0.3+4
[a51aa0fd] Xorg_libXi_jll v1.7.10+4
[d1454406] Xorg_libXinerama_jll v1.1.4+4
[ec84b674] Xorg_libXrandr_jll v1.5.2+4
[ea2f1a96] Xorg_libXrender_jll v0.9.10+4
[14d82f49] Xorg_libpthread_stubs_jll v0.1.0+3
[c7cfdc94] Xorg_libxcb_jll v1.13.0+3
[c5fb5394] Xorg_xtrans_jll v1.4.0+3
[9a68df92] isoband_jll v0.2.3+0
[a4ae2306] libaom_jll v3.4.0+0
[0ac62f75] libass_jll v0.15.1+0
[f638f0a6] libfdk_aac_jll v2.0.2+0
[b53b4c65] libpng_jll v1.6.38+0
[075b6546] libsixel_jll v1.8.6+1
[f27f6e37] libvorbis_jll v1.3.7+1
[1270edf5] x264_jll v2021.5.5+0
[dfaa095f] x265_jll v3.5.0+0
[0dad84c5] ArgTools @stdlib/ArgTools
[56f22d72] Artifacts @stdlib/Artifacts
[2a0f44e3] Base64 @stdlib/Base64
[ade2ca70] Dates @stdlib/Dates
[8bb1440f] DelimitedFiles @stdlib/DelimitedFiles
[8ba89e20] Distributed @stdlib/Distributed
[f43a241f] Downloads @stdlib/Downloads
[7b1f6079] FileWatching @stdlib/FileWatching
[9fa8497b] Future @stdlib/Future
[b77e0a4c] InteractiveUtils @stdlib/InteractiveUtils
[4af54fe1] LazyArtifacts @stdlib/LazyArtifacts
[b27032c2] LibCURL @stdlib/LibCURL
[76f85450] LibGit2 @stdlib/LibGit2
[8f399da3] Libdl @stdlib/Libdl
[37e2e46d] LinearAlgebra @stdlib/LinearAlgebra
[56ddb016] Logging @stdlib/Logging
[d6f4376e] Markdown @stdlib/Markdown
[a63ad114] Mmap @stdlib/Mmap
[ca575930] NetworkOptions @stdlib/NetworkOptions
[44cfe95a] Pkg @stdlib/Pkg
[de0858da] Printf @stdlib/Printf
[9abbd945] Profile @stdlib/Profile
[3fa0cd96] REPL @stdlib/REPL
[9a3f8284] Random @stdlib/Random
[ea8e919c] SHA @stdlib/SHA
[9e88b42a] Serialization @stdlib/Serialization
[1a1011a3] SharedArrays @stdlib/SharedArrays
[6462fe0b] Sockets @stdlib/Sockets
[2f01184e] SparseArrays @stdlib/SparseArrays
[10745b16] Statistics @stdlib/Statistics
[4607b0f0] SuiteSparse @stdlib/SuiteSparse
[fa267f1f] TOML @stdlib/TOML
[a4e569a6] Tar @stdlib/Tar
[8dfed614] Test @stdlib/Test
[cf7118a7] UUIDs @stdlib/UUIDs
[4ec0a83e] Unicode @stdlib/Unicode
[e66e0078] CompilerSupportLibraries_jll @stdlib/CompilerSupportLibraries_jll
[deac9b47] LibCURL_jll @stdlib/LibCURL_jll
[29816b5a] LibSSH2_jll @stdlib/LibSSH2_jll
[c8ffd9c3] MbedTLS_jll @stdlib/MbedTLS_jll
[14a3606d] MozillaCACerts_jll @stdlib/MozillaCACerts_jll
[4536629a] OpenBLAS_jll @stdlib/OpenBLAS_jll
[05823500] OpenLibm_jll @stdlib/OpenLibm_jll
[83775a58] Zlib_jll @stdlib/Zlib_jll
[8e850b90] libblastrampoline_jll @stdlib/libblastrampoline_jll
[8e850ede] nghttp2_jll @stdlib/nghttp2_jll
[3f19e933] p7zip_jll @stdlib/p7zip_jll
Precompiling project...
✗ GLMakie
39 dependencies successfully precompiled in 927 seconds (223 already precompiled)
1 dependency errored. To see a full report either run import Pkg; Pkg.precompile() or load the package
Testing Running tests...
┌ Warning: This file does not include all the tests for Molly.jl due to CI time limits, see the test directory for more
└ @ Main ~/.julia/packages/Molly/7GX09/test/runtests.jl:15
WARNING: Makie.MakieLayout is deprecatedThe module MakieLayout has been removed and integrated into Makie, so simply replace all usage of MakieLayout with Makie.
likely near none:1
WARNING: importing deprecated binding Makie.MakieLayout into GLMakie.
WARNING: Wrapping Vararg directly in UnionAll is deprecated (wrap the tuple instead).
libGL error: MESA-LOADER: failed to open crocus: /usr/lib/dri/crocus_dri.so: cannot open shared object file: No such file or directory (search paths /usr/lib/x86_64-linux-gnu/dri:$${ORIGIN}/dri:/usr/lib/dri, suffix _dri)
libGL error: failed to load driver: crocus
libGL error: MESA-LOADER: failed to open swrast: /usr/lib/dri/swrast_dri.so: cannot open shared object file: No such file or directory (search paths /usr/lib/x86_64-linux-gnu/dri:$${ORIGIN}/dri:/usr/lib/dri, suffix _dri)
libGL error: failed to load driver: swrast
┌ Warning: GLFW couldn't create an OpenGL window.
│ This likely means, you don't have an OpenGL capable Graphic Card,
│ or you don't have an OpenGL 3.3 capable video driver installed.
│ Have a look at the troubleshooting section in the GLMakie readme:
│ https://github.com/JuliaPlots/Makie.jl/tree/master/GLMakie#troubleshooting-opengl.
└ @ GLMakie ~/.julia/packages/GLMakie/IhyZ5/src/screen.jl:391
ERROR: LoadError: GLFWError (VERSION_UNAVAILABLE): GLX: Failed to create context: GLXBadFBConfig
Stacktrace:
[1] _ErrorCallbackWrapper(code::Int32, description::Cstring)
@ GLFW ~/.julia/packages/GLFW/BWxfF/src/callback.jl:43
[2] CreateWindow(width::Int64, height::Int64, title::String, monitor::GLFW.Monitor, share::GLFW.Window)
@ GLFW ~/.julia/packages/GLFW/BWxfF/src/glfw3.jl:499
[3] GLFW.Window(; name::String, resolution::Tuple{Int64, Int64}, debugging::Bool, major::Int64, minor::Int64, windowhints::Vector{Tuple{UInt32, Integer}}, contexthints::Vector{Tuple{UInt32, Integer}}, visible::Bool, focus::Bool, fullscreen::Bool, monitor::Nothing, share::GLFW.Window)
@ GLFW ~/.julia/packages/GLFW/BWxfF/src/glfw3.jl:344
[4] GLMakie.Screen(; resolution::Tuple{Int64, Int64}, visible::Bool, title::String, start_renderloop::Bool, kw_args::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ GLMakie ~/.julia/packages/GLMakie/IhyZ5/src/screen.jl:383
[5] top-level scope
@ ~/.julia/packages/GLMakie/IhyZ5/src/precompiles.jl:5
[6] include(mod::Module, _path::String)
@ Base ./Base.jl:418
[7] include
@ ~/.julia/packages/GLMakie/IhyZ5/src/GLMakie.jl:1 [inlined]
[8] macro expansion
@ ~/.julia/packages/GLMakie/IhyZ5/src/precompiles.jl:19 [inlined]
[9] macro expansion
@ ~/.julia/packages/Makie/h8AEO/src/precompiles.jl:96 [inlined]
[10] top-level scope
@ ~/.julia/packages/GLMakie/IhyZ5/src/precompiles.jl:14
[11] include(mod::Module, _path::String)
@ Base ./Base.jl:418
[12] include(x::String)
@ GLMakie ~/.julia/packages/GLMakie/IhyZ5/src/GLMakie.jl:1
[13] top-level scope
@ ~/.julia/packages/GLMakie/IhyZ5/src/GLMakie.jl:65
[14] include
@ ./Base.jl:418 [inlined]
[15] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
@ Base ./loading.jl:1318
[16] top-level scope
@ none:1
[17] eval
@ ./boot.jl:373 [inlined]
[18] eval(x::Expr)
@ Base.MainInclude ./client.jl:453
[19] top-level scope
@ none:1
in expression starting at /home/ian/.julia/packages/Makie/h8AEO/precompile/shared-precompile.jl:4
in expression starting at /home/ian/.julia/packages/GLMakie/IhyZ5/src/precompiles.jl:13
in expression starting at /home/ian/.julia/packages/GLMakie/IhyZ5/src/GLMakie.jl:1
ERROR: LoadError: Failed to precompile GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a] to /home/ian/.julia/compiled/v1.7/GLMakie/jl_BU3qQc.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:33
[2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, ignore_loaded_modules::Bool)
@ Base ./loading.jl:1466
[3] compilecache(pkg::Base.PkgId, path::String)
@ Base ./loading.jl:1410
[4] _require(pkg::Base.PkgId)
@ Base ./loading.jl:1120
[5] require(uuidkey::Base.PkgId)
@ Base ./loading.jl:1013
[6] require(into::Module, mod::Symbol)
@ Base ./loading.jl:997
[7] top-level scope
@ ~/.julia/packages/Molly/7GX09/test/runtests.jl:30
[8] include(fname::String)
@ Base.MainInclude ./client.jl:451
[9] top-level scope
@ none:6
in expression starting at /home/ian/.julia/packages/Molly/7GX09/test/runtests.jl:29
ERROR: Package Molly errored during testing

(@v1.7) pkg>

I encountered the below error when attempting to compile Molly on the above system configurations.
Let me know if there is additional information that may be of use.

julia> import Pkg; Pkg.precompile()
Precompiling project...
  ✗ BinaryProvider
  ✗ Molly
  0 dependencies successfully precompiled in 4 seconds. 152 already precompiled.

ERROR: The following 1 direct dependency failed to precompile:

Molly [aa0f7f06-fcc0-5ec4-a7f3-a573f33f9c4c]

Failed to precompile Molly [aa0f7f06-fcc0-5ec4-a7f3-a573f33f9c4c] to /Users/UserName/.julia/compiled/v1.8/Molly/jl_eOF7rH.
ERROR: LoadError: InitError: UndefVarError: libchemfiles not defined
Stacktrace:
  [1] chfl_version
    @ ~/.julia/packages/Chemfiles/QiqBE/src/generated/cdef.jl:19 [inlined]
  [2] version
    @ ~/.julia/packages/Chemfiles/QiqBE/src/Chemfiles.jl:36 [inlined]
  [3] __init__()
    @ Chemfiles ~/.julia/packages/Chemfiles/QiqBE/src/Chemfiles.jl:145
  [4] _include_from_serialized(pkg::Base.PkgId, path::String, depmods::Vector{Any})
    @ Base ./loading.jl:831
  [5] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt64)
    @ Base ./loading.jl:1039
  [6] _require(pkg::Base.PkgId)
    @ Base ./loading.jl:1315
  [7] _require_prelocked(uuidkey::Base.PkgId)
    @ Base ./loading.jl:1200
  [8] macro expansion
    @ ./loading.jl:1180 [inlined]
  [9] macro expansion
    @ ./lock.jl:223 [inlined]
 [10] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1144
 [11] include
    @ ./Base.jl:422 [inlined]
 [12] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::Nothing)
    @ Base ./loading.jl:1554
 [13] top-level scope
    @ stdin:1
during initialization of module Chemfiles
in expression starting at /Users/UserName/.julia/packages/Molly/GBoZA/src/Molly.jl:1
in expression starting at stdin:1
Stacktrace:
 [1] pkgerror(msg::String)
   @ Pkg.Types /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/Types.jl:67
 [2] precompile(ctx::Pkg.Types.Context, pkgs::Vector{String}; internal_call::Bool, strict::Bool, warn_loaded::Bool, already_instantiated::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Pkg.API /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1427
 [3] precompile
   @ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1058 [inlined]
 [4] #precompile#225
   @ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1057 [inlined]
 [5] precompile (repeats 2 times)
   @ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Pkg/src/API.jl:1057 [inlined]
 [6] top-level scope
   @ REPL[1]:1

I think it would be useful to discuss here the design changes required for working on GPUs.

@jgreener64, could you detail a bit on what's needed for working on GPUs?
I was wondering if KernelAbstractions.jl could help.

I'm not really an expert with how Markdown works or how it parses links but a lot of links in the docs seem to head to https://github.com/JuliaMolSim/Molly.jl/blob/master/docs/src/@ref , when clicking on them it leads to the 404 GitHub page. The links that head to Wikipedia and stuff obviously work but everything that points back to something else in Molly.jl docs seems broken. I don't really know what @ref is supposed to do so I decided to flag an issue instead of fixing it.

I believe there is a sign mistake in function potential_cutoff(cutoff::ShiftedForceCutoff, r2, inter, params) at cutoffs.jl at line 92,
potential(inter, r2, invr2, params) - (r - rc) * fc - potential(inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params) should read potential(inter, r2, invr2, params) + (r - rc) * fc - potential(inter, cutoff.sqdist_cutoff, cutoff.inv_sqdist_cutoff, params). Indeed the correction should subtract the derivative of the potential, not the force. The (exagerated) effect can be seen on the following pictures (the "linear correction" curves):

Current version

Correct version

I also implemented a spline interpolating cutoff, I can open a PR for it if you are interested.

I think releasing a new version would be appropriate at this point. There are some changes to the acceleration and force functions and it would be easier to have a new version for the development in the ParticleAccelerations / NBodySummations future package.

I think registering the package in the JuliaMolSim registry at least could be helpful for the developing packages that need Molly for their tests (such as the upcoming ParticleAccelerations.jl).

I'm not sure how Pkg would deal with a package with the same UUID and name in both General and a private registry.

I've looked at the performance of a standard liquid argon simulation and I'd like to propose some small design changes.

• The interactions could be stored in tuples or named tuples if key based indexing is required. This should be coupled with adding another type parameter (for each) in the structure definition. The idea of this change is to avoid the presence of abstract types in the acceleration loop. With the current design type inference is not able to obtain concrete types and thus there are some allocations that could be avoided.
• The force! functions should return nothing since with the current design they are type unstable as the if branches at the beginning can return nothing. Since these are in-place functions, returning the forces does not seem particularly useful (and it is not used as far as I can tell). Due to the small union optimizations this is not a big performance problem, but it would be nice to improve the situation.

If you agree, I can make PRs to address this changes.

I am interested in differential learning for problems like mean_min_distance. First, great job in developing the package!
How can I use multiple steps from a single simulation to calculate mean_min_distance? when I use logger data, the gradient is nothing, which makes sense. Any help or tutorial is highly appreciated.

I could not get the tests to work on my GTX970 GPU. Seems like there is an issue with

function DistanceVecNeighborFinder(;
                                nb_matrix,
                                matrix_14=falses(size(nb_matrix)),
                                n_steps=10,
                                dist_cutoff)
    n_atoms = size(nb_matrix, 1)
    if isa(nb_matrix, CuArray)
        is = cu(hcat([collect(1:n_atoms) for i in 1:n_atoms]...))
        js = cu(permutedims(is, (2, 1)))
        m14 = cu(matrix_14)
    else
        is = hcat([collect(1:n_atoms) for i in 1:n_atoms]...)
        js = permutedims(is, (2, 1))
        m14 = matrix_14
    end
    return DistanceVecNeighborFinder{typeof(dist_cutoff), typeof(nb_matrix), typ
eof(is)}(
            nb_matrix, m14, n_steps, dist_cutoff, is, js)
end

Specifically when called in test/protein.jl

I think permuteddims doesn't work on a CuArray, so I tried keeping it as an array, but eventually ran into an issue with turning is into an array for the DistanceVecNeighborFinder. I tried a bunch of different variations, so I'll just leave the unchanged error here:

OpenMM protein comparison: Error During Test at /home/leios/projects/Molly.jl/test/protein.jl:54
  Got exception outside of a @test
  MethodError: no method matching iterate(::Nothing)
  Closest candidates are:
    iterate(::Union{LinRange, StepRangeLen}) at ~/builds/julia-1.7.1/share/julia/base/range.jl:826
    iterate(::Union{LinRange, StepRangeLen}, ::Integer) at ~/builds/julia-1.7.1/share/julia/base/range.jl:826
    iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at ~/builds/julia-1.7.1/share/julia/base/dict.jl:695
    ...
  Stacktrace:
    [1] indexed_iterate(I::Nothing, i::Int64)
      @ Base ./tuple.jl:92
    [2] CUDA.MemoryInfo()
      @ CUDA ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:155
    [3] OutOfGPUMemoryError (repeats 2 times)
      @ ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:199 [inlined]
    [4] throw_api_error(res::CUDA.cudaError_enum)
      @ CUDA ~/.julia/packages/CUDA/VWaZ6/lib/cudadrv/error.jl:89
    [5] macro expansion
      @ ~/.julia/packages/CUDA/VWaZ6/lib/cudadrv/error.jl:101 [inlined]
    [6] cuMemAlloc_v2(dptr::Base.RefValue{CuPtr{Nothing}}, bytesize::Int64)
      @ CUDA ~/.julia/packages/CUDA/VWaZ6/lib/utils/call.jl:26
    [7] #alloc#1
      @ ~/.julia/packages/CUDA/VWaZ6/lib/cudadrv/memory.jl:86 [inlined]
    [8] macro expansion
      @ ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:41 [inlined]
    [9] macro expansion
      @ ./timing.jl:299 [inlined]
   [10] actual_alloc(bytes::Int64; async::Bool, stream::CuStream)
      @ CUDA ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:39
   [11] macro expansion
      @ ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:224 [inlined]
   [12] macro expansion
      @ ./timing.jl:299 [inlined]
   [13] #_alloc#204
      @ ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:305 [inlined]
   [14] #alloc#203
      @ ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:291 [inlined]
   [15] alloc
      @ ~/.julia/packages/CUDA/VWaZ6/src/pool.jl:287 [inlined]
   [16] CuArray{Int64, 2, CUDA.Mem.DeviceBuffer}(#unused#::UndefInitializer, dims::Tuple{Int64, Int64})
      @ CUDA ~/.julia/packages/CUDA/VWaZ6/src/array.jl:42
   [17] similar
      @ ~/.julia/packages/CUDA/VWaZ6/src/array.jl:164 [inlined]
   [18] permutedims(B::CuArray{Int64, 2, CUDA.Mem.DeviceBuffer}, perm::Tuple{Int64, Int64})
      @ Base ./multidimensional.jl:1503
   [19] DistanceVecNeighborFinder(; nb_matrix::CuArray{Bool, 2, CUDA.Mem.DeviceBuffer}, matrix_14::CuArray{Bool, 2, CUDA.Mem.DeviceBuffer}, n_steps::Int64, dist_cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}})
      @ Molly ~/projects/Molly.jl/src/neighbors.jl:134
   [20] System(coord_file::String, force_field::OpenMMForceField{Float64, Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, Quantity{Float64, 𝐋^2 𝐌 𝐍^-1 𝐓^-2, Unitful.FreeUnits{(kJ, mol^-1), 𝐋^2 𝐌 𝐍^-1 𝐓^-2, nothing}}, Quantity{Float64, 𝐌 𝐍^-1 𝐓^-2, Unitful.FreeUnits{(kJ, nm^-2, mol^-1), 𝐌 𝐍^-1 𝐓^-2, nothing}}}; velocities::CuArray{SVector{3, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(nm, ps^-1), 𝐋 𝐓^-1, nothing}}}, 1, CUDA.Mem.DeviceBuffer}, box_size::Nothing, loggers::Dict{Any, Any}, units::Bool, gpu::Bool, gpu_diff_safe::Bool, dist_cutoff::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, nl_dist::Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, rename_terminal_res::Bool)
      @ Molly ~/projects/Molly.jl/src/setup.jl:678
   [21] macro expansion
      @ ~/projects/Molly.jl/test/protein.jl:156 [inlined]
   [22] macro expansion
      @ ~/builds/julia-1.7.1/share/julia/stdlib/v1.7/Test/src/Test.jl:1283 [inlined]
   [23] top-level scope
      @ ~/projects/Molly.jl/test/protein.jl:55
   [24] include(fname::String)
      @ Base.MainInclude ./client.jl:451
   [25] top-level scope
      @ ~/projects/Molly.jl/test/runtests.jl:71
   [26] include(fname::String)
      @ Base.MainInclude ./client.jl:451
   [27] top-level scope
      @ none:6
   [28] eval
      @ ./boot.jl:373 [inlined]
   [29] exec_options(opts::Base.JLOptions)
      @ Base ./client.jl:268
   [30] _start()
      @ Base ./client.jl:495

This could be related to #16 , but I felt it was different enough to warrant a separate issue.

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_units functions 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)

I was thinking that it may be useful to consider an option to toggle progress monitoring, as during benchmarks it may not always be useful.
An other thing that may improve user experience would be to switch to ProgressLogging.jl as this could show progress in IDEs (such as Juno) and also the terminal.

I'm excited that Molecular dynamics is implemented in julia so I'd like to test it!
But I struggled with finding inplicit units that Molly uses and I found them in GROMACS website.

I know Molly targets GROMACS file formats and units, but I don't want to find Molly documents in GROMACS website.
Also clear unit descriptiion would be helpful when non-GROMACS users (actually it's me!) try to convert their non-GROMACS file formats suitable for Molly.
Or redirecting GROMACS unit description page will save user's time
https://manual.gromacs.org/documentation/current/reference-manual/topologies/topology-file-formats.html

Thank you for your effort.

The docs link is broken!

I have been watching Molly.jl and its fast development for a while and really like it. Recently decided to use it for my new research project and first test results are promising. See a simple model of molecular motor walking along a polymer:

https://nxcloud.omid.land/s/czXWqnM8KZFFxDK

However Molly.jl still lacks some essential features of a proper MD package, like FENE bonds and cosine angles. I already implemented them for my own research, if nobody else is working on those currently, I can do some code styling and make a pull request to add them to Molly.

Hello,

I ran the code at the beginning of the docs. I get a bounds error:
BoundsError: attempt to access 0-element Vector{Any} at index [1]
because at line 676 of types.jl the system tries to assign each atom a symbol:
AtomsBase.atomic_symbol(s::Union{System, ReplicaSystem}, i::Integer) = Symbol(s.atoms_data[i].element)
However, atom_data is not a required parameter in the System struct so something seems wrong here to me.

Am I doing something wrong? The entire code is below:

using Molly

n_atoms = 100
atom_mass = 10.0u"u"
atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]

boundary = CubicBoundary(2.0u"nm", 2.0u"nm", 2.0u"nm") # Periodic boundary conditions
coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm") # Random placement without clashing

temp = 100.0u"K"
velocities = [velocity(atom_mass, temp) for i in 1:n_atoms]

pairwise_inters = (LennardJones(),) # Don't forget the trailing comma!

sys = System(
    atoms=atoms,
    pairwise_inters=pairwise_inters,
    coords=coords,
    velocities=velocities,
    boundary=boundary,
    loggers=(
        temp=TemperatureLogger(10),
        coords=CoordinateLogger(10),
    ),
)

simulator = VelocityVerlet(
    dt=0.002u"ps",
    coupling=AndersenThermostat(temp, 1.0u"ps"),
)

simulate!(sys, simulator, 1_000)

For example in the first example, place_atoms should be placeatoms. And AndersenThermostat has wrong number of variables, etc...

n_atoms = 100
box_size = SVector(2.0, 2.0, 2.0)u"nm"
temp = 298.0u"K"
atom_mass = 10.0u"u"

atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]
coords = place_atoms(n_atoms, box_size, 0.3u"nm")
velocities = [velocity(atom_mass, temp) for i in 1:n_atoms]
general_inters = (LennardJones(),)
simulator = VelocityVerlet(dt=0.002u"ps", coupling=AndersenThermostat(temp, 1.0u"ps"))

sys = System(
atoms=atoms,
general_inters=general_inters,
coords=coords,
velocities=velocities,
box_size=box_size,
loggers=Dict("temp" => TemperatureLogger(100)),
)

simulate!(sys, simulator, 10_000)

Just a remark that I noted that the latest AtomsBase interface (0.3) is not properly supported by Molly. E.g. hasatomkey(sys, :testkey) on a Molly system gives an infinite stacktrace.

We should take a look at this and integrate with AtomsBaseTesting in the CI of Molly.

Would it please be possible to add an example of neural network as the simulator?

Also, how can we make a nice Gif?

There are a couple things that the System implementation is missing for it to be fully interchangable with other AtomsBase implementations.

1. atomic_symbol(::System, ::Integer) is not implemented
2. getindex(::System, ::Integer) returns a struct with insufficient fields to be useful in the way that iterating over other AbstractSystem is useful. Specifically this is a problem in my Atomistic and DFTK integration code where I have to have a special case for a System versus an AbstractSystem. This could be resolved simply by having getindex return an AtomView (though any current usage of getindex would have to be modified).

To get a sense of the functionality that Molly provides, I followed the first example in the documentation. At the end of this example I use the visualize function to use Makie to export the coordinates to a gif. However, the gif that I get scrolls as it animates.

I'm not sure if this is an issue with Molly or something underlying with Makie, but does anyone have any thoughts on this?

Following the simulating a gas example,
[Well written documentation, by the way !]
the visualize step fails:

visualize(s.loggers["coords"], box_size, "sim_lj.gif")
WARNING: both Makie and CUDA export "record"; uses of it in module Molly must be qualified

ERROR: UndefVarError: record not defined
Stacktrace:
 [1] visualize(::CoordinateLogger{SArray{Tuple{3},Float64,1,3}}, ::Float64, ::String; connections::Array{Tuple{Int64,Int64},1}, connection_frames::Array{BitArray{1},1}, trails::Int64, framerate::Int64, color::Symbol, connection_color::Symbol, markersize::Float64, linewidth::Float64, transparency::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/ederag/.julia/environments/v1.5/dev/Molly/src/makie.jl:77
 [2] visualize(::CoordinateLogger{SArray{Tuple{3},Float64,1,3}}, ::Float64, ::String) at /home/ederag/.julia/environments/v1.5/dev/Molly/src/makie.jl:19
 [3] top-level scope at REPL[3]:1

The kinetic_energy function takes only one argument—neighbors should be removed.

I tried running the example code for the protein.

After a copy and paste of the code (fixing the neighbor -> neighbour typo) and loading the result into VMD I was surprised seeing the atoms jumping around like crazy.

Indeed a look at the temperature logger confirmed that the temperature raised from 260 to 6000 K after the first 500 steps.

I was using Molly for a simulation and initially was using the incorrect units (kg and J). The result of this was that the simulation just hung indefinitely. Initially I thought that it was just taking a long time but after a couple hours of no progress I realized my mistake. Obviously having numbers that are off by 20+ orders of magnitude is an edge case that isn't expected, but perhaps there could be some sort of check to ideally indicate that the input is obviously unphysical or at the very least actually terminate.

Long term, it would be very nice to have integration with a unit package such as Unitful.

Hi! I have implemented some custom pairwise interactions that use the DistanceNeighborFinder. The implementations work on the CPU. I then tried to use the interactions on the GPU. I followed the examples on the Molly documentation, casting the relevant arrays as CuArray and ensuring all the interactions evaluated to true when using the isbitstype function. However, I get an error that is related to the CUDA.jl package:

Error: LoadError: InvalidIRError: compiling MethodInstance for Molly.pairwise_force_kernel!(... ... ...)
Reason: unsupported call to an unknown function (call to ijl_get_nth_field_checked) 

The error only appears when I use three (or more) pairwise interactions together:

pairwise_interactions = (
    Interaction1(
        use_neighbors=true,
    ),
    Interaction2(
        use_neighbors=true,
    ),
    Interaction3(
        use_neighbors=true,
    ),
)

The CUDA kernel is compiled successfully if I use any pair of interactions (interaction1 + interaction2, interaction1 + interaction3, etc.). By looking at the stack trace, I tracked down the error back to these lines 36-39 of the src/cuda.jl file:

f = force_gpu(inters[1], dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
for inter in inters[2:end]
    f += force_gpu(inter, dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special)
end

I still need to familiarize myself more with the CUDA.jl package, but maybe the instructions generated when using tuples of pairwise interactions are not supported yet for generating CUDA kernels. I tried different versions of the same lines of code, hoping to find an alternative representation that is compatible with CUDA.jl, and found that the following way makes the error disappear:

f = sum(force_gpu(inter, dr, coord_i, coord_j, atoms[i], atoms[j], boundary, special) for inter in inters)

By introducing this change, the CI passes all the tests successfully. I have not yet checked the performance implications, but I can try if there are input files for benchmarking the GPU code.

Finally, I can submit a PR if this change looks ok.

UPDATE: I found that using other combinations of pairwise interactions makes the error reappear. I will further investigate the cause...

Hi everyone, I was testing Molly in a new computational resource. I've just tried to reproduce the main example (the fluid in a Lennard Jones potential). An important information is: I was running on Julia 1.8.1.

All was going good,

using Molly

n_atoms = 100
boundary = CubicBoundary(2.0u"nm", 2.0u"nm", 2.0u"nm")
temp = 298.0u"K"
atom_mass = 10.0u"u"

atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]

when the code crashed because I've tried to run this line:

coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm")

I got the following error:

ERROR: MethodError: no method matching place_atoms(::Int64, ::CubicBoundary{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}; min_dist=0.3 nm)
Closest candidates are:
  place_atoms(::Integer, ::Any, ::Any; max_attempts) at ~/.julia/packages/Molly/RD5GY/src/setup.jl:22 got unsupported keyword argument "min_dist"
Stacktrace:
 [1] top-level scope
   @ REPL[11]:1
 [2] top-level scope
   @ ~/.julia/packages/CUDA/DfvRa/src/initialization.jl:52

As the error said, I took a look in the file setup.jl (line 25).

As you can see, the method place_atoms receives an argument of type Integer, but n_atoms is an Int64. If you go to the Julia REPL, you will see that Int64 is a subtype of Integer, but they are not equal.

julia> Int64 == Integer
false
julia> Int64 <: Integer
true

The solution is to change the syntax of the function:

function place_atoms(n_atoms::T ...) where T <: Integer
...
end

I can help solving it by sending a PR.

All my best regards
Leticia Madureira

Dear Molly developers,

I have been working on ion implantation simulation for a while and want to develop some new algorithms to deal with the microscopic collisions. Then I came across Molly.jl, in particular its gravity examples. I have two questions now, regarding the velocity-dependent forces.

Question 1. How difficult it is?

If I want to modify the function interface of force() in allow for arguments representing (lab frame or relative) velocities, how much effort will I make to get the code running with simulations using the Verlet algorithm?

Question 2. Is it necessary?

Since there isn't a block-diagram sketch of the overall design of Molly.jl, I only grasped a portion of it by digging into the call sequence. I realised that the force() function is designed not to accepting velocity as input argument. It is a deliberate choice? I mean, was it decided at the very beginning of the project that the package don't allow for velocity-dependent forces?

Consider the simulation of galaxy dynamics, the particles are moving in effective dust clouds, experiencing drag forces. I guess, in this case, velocity dependent force are absolutely necessary to make realistic simulations.

In my case of ion implantation simulation, the projectile and the targets are ions. They interact via screened Coulomb potential and a dissipation force depending the relative velocity. The dissipation force has significant impact on the scattering behaviour. In the regime when the projectile's velocity is comparable to the velocity of thermal motions of the targets, one must use molecular dynamics. The velocity-dependent dissipation force is necessary to reproduce the trajectories of the projectile and targets. That is basically why I open this issue to discuss the possibility.

I am an experienced Julia programmer and want to invest on this beautiful project of Molly.jl. Hope to have your opinions :-)

I was just looking through the Lennard-Jones force calculations and noticed that every time you calculate the force between two atoms, you are calculating the mixing rules

σ = sqrt(atom_i.σ * atom_j.σ) 
ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

Given the simplicity of the force calculation, these two square roots likely are adding a considerable amount of flops to the simulation. Typically simulation codes will precalculate the mixing rules and then simply look up the necessary sigma and epsilon, i.e.,

σ = vdw_table.σ(atom_i.type,  atom_j.type) 
ϵ = vdw_table.ϵ(atom_i.type,  atom_j.type)

I do not know if type is available, but, incorporating it into the atom struct so that you can use this parameter lookup table, will likely speed up the calculation a reasonable amount. The type would be an integer, so it should be isbits true and so shouldn't interfere with CUDA/GPU. This will also make the structures lighter since only type needs to be saved per atom, and not epsilon and sigma (also, other force-fields, such as EXP-6, have more parameters so you save even more with them).

Also, usually the size parameters use the Lorenz mixing rule i.e., are

σ = (atom_i.σ + atom_j.σ) / 2

It would be good to add the feature allowing the user to select the mixing rules to use for each parameter type.

I have in mind two possible target systems to use as an example for the membrane barostat:

• DOPC bilayer using the CHARMM force field (all-atom)
• DOPC bilayer using the MARTINI 3 force field (coarse-grained)

Each of these examples requires further development:

• For the CHARMM force field: PME, cmap interaction, table potential, CHARMM .prm parser
• For the MARTINI 3 force field: GROMACS .top parser, GROMACS .gro parser

The easiest one is implementing a MARTINI 3 example. I can work on the required parsers and send a PR. @jgreener64, what do you think?

On the other hand, I am currently using the barostat with my own coarse-grained model. As I mentioned to you before, I am trying to improve it by using Molly's differentiable simulation. I can add an example based on my results after the publication.

Hi,
I tried implementing a custom coupling function analogue to the documentation but did not get it to work.
The issue seems to be, that the apply_coupling method is not found, I get a "ERROR: LoadError: MethodError: no method matching apply_coupling!".
Can you help me with that? Did I miss something there?
Thanks!

Full error message:

ERROR: LoadError: MethodError: no method matching apply_coupling!(::System{3, false, Float64, Vector{Atom{Float64, Quantity{Float64, 𝐌, Unit
ful.FreeUnits{(u,), 𝐌, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, Quantity{Float64, 𝐋^2 𝐌 𝐍^-1 𝐓^-2, Unitful.Fre
eUnits{(kJ, mol^-1), 𝐋^2 𝐌 𝐍^-1 𝐓^-2, nothing}}}}, Vector{SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}, CubicBoun
dary{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}, Vector{SVector{3, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(nm, ps^-1),
 𝐋 𝐓^-1, nothing}}}}, Vector{Any}, Nothing, Tuple{LennardJones{false, NoCutoff, Int64, Int64, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), 𝐋 𝐌 𝐍^-
1 𝐓^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), 𝐋^2 𝐌 𝐍^-1 𝐓^-2, nothing}}}, Tuple{}, Tuple{}, Tuple{}, NoNeighborFinder, NamedTuple{(:tem
p, :coords), Tuple{GeneralObservableLogger{Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}, typeof(Molly.temperature_wrapper)}, Ge
neralObservableLogger{Vector{SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}, typeof(Molly.coordinates_wrapper)}}}, 
Quantity{Float64, 𝐋^2 𝐌 𝚯^-1 𝐓^-2, Unitful.FreeUnits{(kJ, K^-1), 𝐋^2 𝐌 𝚯^-1 𝐓^-2, nothing}}, Unitful.FreeUnits{(kJ, nm^-1, mol^-1), 𝐋 𝐌 𝐍^-1
 𝐓^-2, nothing}, Unitful.FreeUnits{(kJ, mol^-1), 𝐋^2 𝐌 𝐍^-1 𝐓^-2, nothing}, Vector{Quantity{Float64, 𝐌, Unitful.FreeUnits{(u,), 𝐌, nothing}}
}}, ::MyCoupler, ::VelocityVerlet{Quantity{Float64, 𝐓, Unitful.FreeUnits{(ps,), 𝐓, nothing}}, Tuple{AndersenThermostat{Quantity{Float64, 𝚯, 
Unitful.FreeUnits{(K,), 𝚯, nothing}}, Quantity{Float64, 𝐓, Unitful.FreeUnits{(ps,), 𝐓, nothing}}}, MonteCarloBarostat{Float64, Quantity{Floa
t64, 𝐌 𝐋^-1 𝐓^-2, Unitful.FreeUnits{(bar,), 𝐌 𝐋^-1 𝐓^-2, nothing}}, Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, nothing}}, Quantity{Floa
t64, 𝐋^3, Unitful.FreeUnits{(nm^3,), 𝐋^3, nothing}}}, MyCoupler}}, ::Nothing, ::Int64; n_threads::Int64)

Closest candidates are:
  apply_coupling!(::Any, ::Union{Tuple, NamedTuple}, ::Any, ::Any, ::Any; kwargs...)
   @ Molly ~/miniconda3/share/julia/packages/Molly/SRXWU/src/coupling.jl:23
  apply_coupling!(::Any, ::RescaleThermostat, ::Any, ::Any, ::Integer; n_threads)
   @ Molly ~/miniconda3/share/julia/packages/Molly/SRXWU/src/coupling.jl:100
  apply_coupling!(::Any, ::BerendsenThermostat, ::Any, ::Any, ::Integer; n_threads)
   @ Molly ~/miniconda3/share/julia/packages/Molly/SRXWU/src/coupling.jl:124
  ...

I tried to provide a minimal (not working) example from code blocks in the documentation:

using Molly

struct MyCoupler
    # Any properties, e.g. a target temperature or coupling constant
end

function apply_coupling!(sys, coupling::MyCoupler, sim, neighbors, step_n;
    n_threads=Threads.nthreads())
# Do something to the simulation, e.g. scale the velocities
# Return whether the coupling has invalidated the currently stored forces,
#   for example by changing the coordinates
recompute_forces = false
return recompute_forces
end
custom_coupler = MyCoupler()

n_atoms = 100
atom_mass = 10.0u"u"
atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]
boundary = CubicBoundary(2.0u"nm") # Periodic boundary conditions with a 2 nm cube
coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm") # Random placement without clashing

temp = 100.0u"K"
velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms]

pairwise_inters = (LennardJones(),) # Don't forget the trailing comma!

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=pairwise_inters,
)

simulator = VelocityVerlet(dt=0.001u"ps", coupling=custom_coupler)

simulate!(sys, simulator, 1_000)

I am using Molly.jl for my doctoral project and want to add a citation.
Do you have a preferred way of citation?
Have you considered adding a "CITATION" file?
(https://docs.github.com/en/repositories/managing-your-repositorys-settings-and-features/customizing-your-repository/about-citation-files)
Thanks a lot!

Given a System I want to obtain a new System without mutation. Therefore I thought of using Accessors.jl.
However, it turns out that Accessors.set cannot find the constructor of System.

I then tried myself

sys = [... construct some system, e.g. from the example of Molly's README.md]
fields = map(fieldnames(System)) do n getfield(sys, n) end
System(fields...) # ERROR: MethodError: no method matching
System{map(typeof, fields)...}(fields...)  # ERROR: MethodError: no method matching

From what I understand the constructor should be able to be found, right?
Can someone confirm this problem? Is this a problem with Molly?

This issue is a roadmap for Molly.jl development. Feel free to discuss things here or submit a PR. Bear in mind that significant refactoring will probably occur as the package develops.

Low hanging fruit

Want to get involved? These issues might be the place to start:

• Speedups to the code. Both on the level of individual functions and the overall algorithm. [Getting there]
• Look over the code to find mistakes in implementation. There will be some because...
• Test coverage is minimal. We need tests for the 'core' functions that are unlikely to change their return values, e.g. the force functions. [Done]

Bigger things to tackle

• Energy minimisation by gradient descent [Done].
• Add other forcefields. [Partly done]
• Add other thermostats/ensembles. [Partly done]
• Get it working on the GPU. [Done]
• Line by line checking of how a mature MD package works to see where we have missed things and to get an idea of how to organise the codebase in future. [Partly done]
• Neighbour list in a cubic lattice, and use of group lists. [Partly done]
• Setup of MD simulations, e.g. adding hydrogen, solvent. [Partly done]
• Trajectory/topology file format readers and writers. The current topology reader is not very robust. [Partly done]

Projects that will require more planning or will cause significant changes

• Abstracting the forcefield from the integration to allow different forcefields and integrators to be combined. [Done]
• Allow definition of custom forcefields. [Done]
• Parallelisation. [Done]

Okay, so this is less of a roadmap and more of a list of ideas to work on. But most or all of the above would be required to make this a package suitable for general use, so in a sense this is the roadmap.

Hi,
I'm trying to use Molly.jl to reproduce some of the work I have done in Gromacs, however, it raises error in System building stage.
Specifically, I use the following Gromacs commands to prepare the pdb file(with water box and ions):

gmx pdb2gmx -ff amber99sb -f input.pdb -water spce -ignh -o pro.pdb #-ff forcefiled; -f input file
gmx editconf -f pro.pdb -o pro-box.pdb -d 1.0 -bt octahedron #Add box
gmx solvate -cp pro-box.pdb -cs spc216.gro -p topol.top -o pro-sol.gro #Add solvate
gmx grompp -f em-steep.mdp -c pro-sol.gro -p topol.top -o ion.tpr -maxwarn 1 
echo 14 | gmx genion -s ion.tpr -neutral -conc 0.15 -p topol.top -o pro-ion.pdb 

And use the following Molly.jl code to generate the system:

using Molly
ff = MolecularForceField("ff99SB.xml", "spce_standard.xml")
sys = System("pro-ion.pdb", ff)

the ff99SB.xml and spce_standard.xml files are found from openmm source code.

This gives me the following error:

ERROR: LoadError: KeyError: key "HB1" not found

and a lot of warning:

┌ Warning: PDB reader: found unexpected, non-standard atom 'CD1' in residue 'ILE' (resid 527)
└ @ Chemfiles ~/.julia/packages/Chemfiles/NoLfC/src/misc.jl:29
┌ Warning: PDB reader: found unexpected, non-standard atom 'CD1' in residue 'ILE' (resid 527)
└ @ Chemfiles ~/.julia/packages/Chemfiles/NoLfC/src/misc.jl:29
┌ Warning: PDB reader: found unexpected, non-standard atom 'CD1' in residue 'ILE' (resid 527)
└ @ Chemfiles ~/.julia/packages/Chemfiles/NoLfC/src/misc.jl:29
┌ Warning: PDB reader: found unexpected, non-standard atom 'O' in residue 'LYS' (resid 528)
└ @ Chemfiles ~/.julia/packages/Chemfiles/NoLfC/src/misc.jl:29
┌ Warning: PDB reader: found unexpected, non-standard atom 'O' in residue 'LYS' (resid 528)

...

I know the process of combining openmm, Gromacs and Molly.jl toolchain is kind of error prone. However, since Molly.jl currently does not have toolset for processing pdb files(creating water box, adding ions, etc), this is the best I come up with...

Do you have any clues about what's going on here? Any suggestions on pdb preparation for Molly.jl?

At the minute the GPU/differentiable path is Zygote-compatible and hence uses non-mutating broadcasted operations. This works, but is rather slow and very GPU memory-intensive.

Long term the plan is to switch to Enzyme-compatible GPU kernels to calculate and sum the forces using the neighbour list. This will be much faster both with and without gradients, and should help us move towards the speeds of existing MD software. These kernels could be used as part of the general interaction interface as is, or another interface could emerge to use them. Enzyme and Zygote can be used together, so it should be possible to replace the force summation alone and retain the functionality of the package.

One consideration is how general such kernels should be. A general pairwise force summation kernel for user-defined force functions would be useful for Lennard-Jones and Coulomb interactions, and hence would be sufficient for macromolecular simulation. Other more specialised multi-body kernels could live in Molly or elsewhere depending on how generic they are.

Another concern is how the neighbour list is best stored (calculation of the neighbour list can also be GPU accelerated but that is a somewhat separate issue).

Something to bear in mind is the extension from using one to multiple GPUs for the same simulation. It is probably best to start with one GPU and go from there.

This issue is to track and discuss this development. @leios

Useful links:

• A N-body Lennard-Jones simulator in Julia: https://gist.github.com/AlfTetzlaff/880ad7bf4e1486d8a051580a6262491c
• NVIDIA N-body discussion: https://developer.nvidia.com/gpugems/gpugems3/part-v-physics-simulation/chapter-31-fast-n-body-simulation-cuda
• CellListMap.jl: https://github.com/m3g/CellListMap.jl

Right now the interaction in the library are passed down individually as keyword arguments pairwise_inters, specific_inters_list and general_inters. I have a suggenstion that maybe these can be wrapped in a structure like:

const struct Interactions{PW, SF, GN}
    pairwise_inters::PW
    specific_inters_list::SF
    general_inters::GN
end

This can be helpful in grouping together interactions that are meant for a particular simulation/system without having to worry about x3 variables that carry this information. The usefulness of this is not that much when we are talking about simulation with a single type of interactions, but becomes apparent in cases such as Hamiltonian-REMD (which I am trying to implement) where all the replicas have different force fields.

Although we can always use three different vectors to store the different variations of three interactions whenever we want to do such a thing, it seems a bit non-structured (pun intended) way of doing things, when the i'th elements are clearly meant to be used together.

Still, this will be unnecessary if such a situation is rare in practice (which I am not sure about) and also this will increase complexity a bit for simple simulations like with simple LennardJones potential where you will have to do:

inters = Interactions(pairwise_inters=(LennardJones(..), ))

...

Here a workaround can be to write simple convenience functions that set this up for the user:

LennardJonesInteraction(...) = Interactions(pairwise_inters=(LennardJones(..), ))

Please give your views on this.

I noticed that the current tests only check if the simulation doesn't error, but not if the results are wrong. Since I'm refactoring several internal components, I think I should be careful not to change the results. To this end I was thinking to check for some general properties of the simulation to ensure that it is correct.
I was thinking of checking for the energy of the system. Do you have any other ideas?

The error might be on my side because I am new to Julia, but the following line has worked for me thus far and in the case of Molly it produces an error:

julia> import Pkg; Pkg.add("Molly");
ERROR: LoadError: UndefVarError: dualize not defined
Stacktrace:       
 [1] include(mod::Module, _path::String)
   @ Base .\Base.jl:419
 [2] include(x::String)
   @ Molly C:\Users\idab\.julia\packages\Molly\6TwMF\src\Molly.jl:1
 [3] top-level scope
   @ C:\Users\idab\.julia\packages\Molly\6TwMF\src\Molly.jl:62     
 [4] include
   @ .\Base.jl:419 [inlined]
 [5] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
   @ Base .\loading.jl:1554
 [6] top-level scope
   @ stdin:1
in expression starting at C:\Users\idab\.julia\packages\Molly\6TwMF\src\zygote.jl:5
in expression starting at C:\Users\idab\.julia\packages\Molly\6TwMF\src\Molly.jl:1
in expression starting at stdin:1
ERROR: Failed to precompile Molly [aa0f7f06-fcc0-5ec4-a7f3-a573f33f9c4c] to C:\Users\idab\.julia\compiled\v1.8\Molly\jl_985.tmp.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
    @ Base .\loading.jl:1707
  [3] compilecache
    @ .\loading.jl:1651 [inlined]
  [4] _require(pkg::Base.PkgId)
    @ Base .\loading.jl:1337
  [5] _require_prelocked(uuidkey::Base.PkgId)
    @ Base .\loading.jl:1200
  [6] macro expansion
    @ .\loading.jl:1180 [inlined]
  [7] macro expansion
    @ .\lock.jl:223 [inlined]
  [8] require(into::Module, mod::Symbol)
    @ Base .\loading.jl:1144
  [9] eval
    @ .\boot.jl:368 [inlined]
 [10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1428
 [11] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base .\essentials.jl:729
 [12] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base .\essentials.jl:726
 [13] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\eval.jl:233
 [14] (::VSCodeServer.var"#66#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\eval.jl:157
 [15] withpath(f::VSCodeServer.var"#66#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\repl.jl:249
 [16] (::VSCodeServer.var"#65#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\eval.jl:155
 [17] hideprompt(f::VSCodeServer.var"#65#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [18] (::VSCodeServer.var"#64#68"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\eval.jl:126
 [19] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:511
 [20] with_logger
    @ .\logging.jl:623 [inlined]
 [21] (::VSCodeServer.var"#63#67"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\eval.jl:225
 [22] #invokelatest#2
    @ .\essentials.jl:729 [inlined]
 [23] invokelatest(::Any)
    @ Base .\essentials.jl:726
 [24] macro expansion
    @ c:\Users\idab\.vscode\extensions\julialang.language-julia-1.38.2\scripts\packages\VSCodeServer\src\eval.jl:34 [inlined]
 [25] (::VSCodeServer.var"#61#62")()
    @ VSCodeServer .\task.jl:484

Some additional info:

julia> versioninfo()
Julia Version 1.8.5
Commit 17cfb8e65e (2023-01-08 06:45 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 12 × Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
  Threads: 1 on 12 virtual cores
Environment:
  JULIA_EDITOR = code

Any interest in making a LAMMPS feature parity list similar to what Flux has with PyTorch? Or honestly just an updated list of missing features I think would be good. I definitely have my own ideas, and I think it would be nice to have a publicly available list of features Molly has and does not have. Maybe separate it into two lists: one to match parity with LAMMPS and another that is features only in Molly.