Simplify Molly interface

Project.toml

@@ -10,6 +10,7 @@ AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
 FLoops = "cc61a311-1640-44b5-9fba-1b764f453329"
 Folds = "41a02a25-b8f0-4f67-bc48-60067656b558"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NeighbourLists = "2fcf5ba9-9ed4-57cf-b73f-ff513e316b9c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

docs/src/molly.md

@@ -14,51 +14,38 @@ There are couple of notes that need to be understood.
 using Molly
 using ACEmd
 
-using AtomsBase
 using ExtXYZ
 using Unitful
 
 
-# Load ACE data and initial structure
+# Load ACE model and initial structure
 fname_ace = joinpath(pkgdir(ACEmd), "data", "TiAl.json")
 fname_xyz = joinpath(pkgdir(ACEmd), "data", "TiAl-big.xyz")
 
-data = FastSystem(ExtXYZ.Atoms(read_frame(fname_xyz)))
+data = ExtXYZ.Atoms(read_frame(fname_xyz))
 pot = load_ace_model(fname_ace)
 
-# Prepare data to Molly compatible format
-atoms = [Molly.Atom( index=i, mass=atomic_mass(data, i) ) for i in 1:length(data) ]
-
-boundary = begin
-    box = bounding_box(data)
-    CubicBoundary(box[1][1], box[2][2], box[3][3])
-end
-
-atom_data = [ (; :Z=>z,:element=>s)  for (z,s) in zip(atomic_number(data), atomic_symbol(data))  ]
+# Pack data to Molly compatible format
+# data is AtomsBase compatible structure
+sys = Molly.System(data, pot)
 
 # Set up temperature and velocities
 temp = 298.0u"K"
-velocities = [random_velocity(m, temp) for m in atomic_mass(data)]
+vel = random_velocities(sys, temp)
+
+# Add initial velocities and loggers
+sys = Molly.System(
+    sys;
+    velocities = vel,
+    loggers=(temp=TemperatureLogger(100),)
+)
 
 # Set up simulator
 simulator = VelocityVerlet(
     dt=1.0u"fs",
     coupling=AndersenThermostat(temp, 1.0u"ps"),
 )
 
-# Set up Molly system
-sys = System(
-           atoms=atoms,
-           atoms_data = atom_data,
-           coords=position(data),
-           velocities=velocities,
-           general_inters = (pot,),
-           boundary=boundary,
-           loggers=(temp=TemperatureLogger(100),),
-           energy_units=u"eV",
-           force_units=u"eV/Å",
-       )
-
 
 # Perform MD
 simulate!(sys, simulator, 1000)

ext/ACE_Molly_ext.jl

@@ -2,6 +2,8 @@ module ACE_Molly_ext
 
 import ACE1: AtomicNumber
 using ACEmd
+using AtomsBase
+using LinearAlgebra: isdiag
 using Molly
 
 
@@ -31,4 +33,48 @@ function ACEmd._atomic_number(sys::Molly.System, i)
     return AtomicNumber( sys.atoms_data[i].Z )
 end
 
+function Molly.System(
+        sys::AbstractSystem,
+        pot::ACEpotential;
+        energy_units = pot.energy_unit,
+        force_units = pot.energy_unit/pot.length_unit,
+        velocity_units = pot.length_unit/u"fs",
+        kwargs...
+    )
+    @assert dimension(energy_units) == dimension(u"J") "energy unit has wrong dimenstions"
+    @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) ]
+
+    boundary = begin
+        box = bounding_box(sys)
+        if isdiag( hcat(box...) )
+           tmp = CubicBoundary(box[1][1], box[2][2], box[3][3])
+        else
+           tmp = TriclinicBoundary(box...)
+        end
+        tmp
+    end
+
+    atom_data = [ (; :Z=>z,:element=>s)  for (z,s) in zip(atomic_number(sys), atomic_symbol(sys))  ]
+
+    return Molly.System(
+        atoms=atoms,
+        atoms_data = atom_data,
+        coords= map(sys) do r
+            SVector(position(r)...)
+        end,
+        general_inters = (pot,),
+        boundary=boundary,
+        energy_units=energy_units,
+        force_units=force_units,
+        velocities = map(sys) do a
+            tmp = SVector(velocity(a)...)
+            uconvert.(velocity_units, tmp)
+        end,
+        kwargs...
+    )
+end
+
 end # Module

test/runtests.jl

@@ -84,24 +84,16 @@ end
 
 @testset "Molly support" begin
     pot = load_ace_model(fname_ace)
-    data = FastSystem(ExtXYZ.Atoms(read_frame(fname_xyz)))
+    data = ExtXYZ.Atoms(read_frame(fname_xyz))
 
-    atoms = [Molly.Atom( index=i, mass=AtomsBase.atomic_mass(data,i) ) for i in 1:length(data) ]
-    boundary = begin
-        box = bounding_box(data)
-        CubicBoundary(box[1][1], box[2][2], box[3][3])
-    end
-    adata = [ (; :Z=>z,:element=>s)  for (z,s) in zip(AtomsBase.atomic_number(data), AtomsBase.atomic_symbol(data))  ]
-
-    sys = System(
-           atoms=atoms,
-           atoms_data = adata,
-           coords=position(data),
-           general_inters = (pot,),
-           boundary=boundary,
-           energy_units=u"eV",
-           force_units=u"eV/Å",
-       )
+    sys = Molly.System(data, pot)
+
     @test ace_energy(pot, data)  ≈ Molly.potential_energy(sys)
     @test all( ace_forces(pot, data) .≈ Molly.forces(sys) )
+
+    simulator = VelocityVerlet(
+        dt=1.0u"fs",
+        coupling=AndersenThermostat(300u"K", 1.0u"ps"),
+    )
+    simulate!(sys, simulator, 10)
 end