Update Molly support to v0.19 (#63)

* basic changes for Molly support

* add symbols to molly data

* Fix compat for ACE1 and ACE1x

* restrict to ACE1 v0.12.1

* update Molly docs

Project.toml

@@ -31,15 +31,15 @@ ACE_Molly_ext = "Molly"
 ACEmd_ACEfit_ext = "ACEfit"
 
 [compat]
-ACE1 = "0.11.5, 0.12"
+ACE1 = "0.12.1"
 ACE1x = "0.1.8"
 ACEbase = "0.4.3"
 ACEfit = "0.1.4"
 AtomsBase = "0.3"
 AtomsCalculators = "0.1.1"
 ChunkSplitters = "2"
 Folds = "0.2"
-Molly = "0.15, 0.16, 0.17, 0.18"
+Molly = "0.19"
 NeighbourLists = "0.5"
 Reexport = "1"
 StaticArrays = "1"

docs/src/molly.md

@@ -3,10 +3,7 @@
 ACE support for [Molly](https://github.com/JuliaMolSim/Molly.jl) is loaded with Julia package extensions.
 To use you need to use at least v1.9 of Julia. You also need to have both Molly and ACEapi added to the environment you are using, for extension to be loaded.
 
-There are couple of notes that need to be understood.
-
-- Default unit for energy in Molly is kJ/mol. You need to change this to anything that is not per mole (same for force units)
-- As of writing this Molly does not have fully compatible AtomsBase interface. So building `System` structure is a bit clunky. Most of all, `atoms_data` needs to have defined `element` for atomic symbol and `Z` for nuclear charge. As these are needed by ACE.
+You need to understand that by default Molly uses kJ/mol for energy. This is not compatible with ACEmd and AtomsCalculators, so you need to change the energy unit. ACEmd adds convenience functions to do this for you and at the same time initialize the system.
 
 ## Example
 
@@ -25,8 +22,8 @@ fname_xyz = joinpath(pkgdir(ACEmd), "data", "TiAl-big.xyz")
 data = ExtXYZ.Atoms(read_frame(fname_xyz))
 pot = load_ace_model(fname_ace)
 
-# Pack data to Molly compatible format
-# data is AtomsBase compatible structure
+# Insert the potential to Molly and change the units
+# to reflect this. Also initializes system structure.
 sys = Molly.System(data, pot)
 
 # Set up temperature and velocities

ext/ACE_Molly_ext.jl

@@ -6,33 +6,8 @@ using AtomsBase
 using LinearAlgebra: isdiag
 using Molly
 
-
-function Molly.forces(
-        acp::ACEpotential,
-        sys,
-        neighbors=nothing;
-        n_threads=Threads.nthreads(),
-        executor=ThreadedEx()
-    )
-    return ace_forces(acp, sys; executor=executor, ntasks=n_threads) 
-end
-
-
-function Molly.potential_energy(
-        acp::ACEpotential,
-        sys,
-        neighbors=nothing;
-        n_threads=nothing,
-        executor=ThreadedEx()
-    )
-    return ace_energy(acp, sys; executor=executor) 
-end
 
 
-function ACEmd._atomic_number(sys::Molly.System, i) 
-    return AtomicNumber( sys.atoms_data[i].Z )
-end
-
 function Molly.System(
         sys::AbstractSystem,
         pot::ACEpotential;
@@ -45,7 +20,8 @@ function Molly.System(
     @assert dimension(force_units) == dimension(u"kg*m/s^2")  "force unit has wrong dimenstions"
     @assert dimension(velocity_units) == dimension(u"m/s")  "velocity unit has wrong dimenstions"
 
-    atoms = [Molly.Atom( index=i, mass=atomic_mass(sys, i) ) for i in 1:length(sys) ]
+    atoms = [Molly.Atom(; index=i, mass=atomic_mass(sys, i) ) for i in 1:length(sys) ]
+    atoms_data = [ Molly.AtomData(; element=String(atomic_symbol(sys,i))) for i in 1:length(sys)]
 
     boundary = begin
         box = bounding_box(sys)
@@ -56,12 +32,10 @@ function Molly.System(
         end
         tmp
     end
-
-    atom_data = [ (; :Z=>z,:element=>s)  for (z,s) in zip(atomic_number(sys), atomic_symbol(sys))  ]
 
-    return Molly.System(
+    return Molly.System(;
         atoms=atoms,
-        atoms_data = atom_data,
+        atoms_data=atoms_data,
         coords= map(sys) do r
             SVector(position(r)...)
         end,

src/utils.jl

@@ -106,4 +106,4 @@ end
 
 
 _atomic_number(at::ACE1.Atoms, i)          = at.Z[i]
-_atomic_number(ab::AbstractSystem, i) = ACE1.AtomicNumber(AtomsBase.atomic_number(ab,i))
+_atomic_number(ab, i) = ACE1.AtomicNumber(AtomsBase.atomic_number(ab,i))