remove AtomsBase charge unit

Project.toml

@@ -33,6 +33,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 UnitfulChainRules = "f31437dd-25a7-4345-875f-756556e6935d"
 UnsafeAtomicsLLVM = "d80eeb9a-aca5-4d75-85e5-170c8b632249"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
@@ -63,6 +64,7 @@ Requires = "1"
 SimpleCrystals = "0.1.6"
 StaticArrays = "1.5.17"
 Unitful = "1"
+UnitfulAtomic = "1"
 UnitfulChainRules = "0.1.2"
 UnsafeAtomicsLLVM = "0.1"
 Zygote = "0.6.56"

docs/src/differentiable.md

@@ -43,7 +43,7 @@ end
 ```
 Now we can set up and run the simulation in a similar way to that described in the [Molly documentation](@ref).
 The difference is that we wrap the simulation in a `loss` function.
-This returns a single value that we want to obtain gradients with respect to, in this case the value of the above function at the end of the simulation.
+This returns a single value that we want to obtain gradients with respect to, in this case the difference between the value of the above function at the end of the simulation and a target distance.
 The `Zygote.ignore()` block allows us to ignore code for the purposes of obtaining gradients; you could add the [`visualize`](@ref) function there for example.
 ```julia
 using Zygote

src/Molly.jl

@@ -27,6 +27,7 @@ using PeriodicTable
 using Requires
 using SimpleCrystals
 using Unitful
+using UnitfulAtomic
 using UnitfulChainRules
 using UnsafeAtomicsLLVM
 using Zygote

src/energy.jl

@@ -47,7 +47,6 @@ function temperature(sys)
     end
 end
 
-
 """
     potential_energy(system, neighbors=nothing; n_threads=Threads.nthreads())
 

src/force.jl

@@ -10,7 +10,6 @@ export
     SpecificForce4Atoms,
     forces
 
-
 """
     accelerations(system, neighbors=nothing; n_threads=Threads.nthreads())
 

src/types.jl

@@ -969,12 +969,11 @@ AtomsBase.atomic_mass(s::Union{System, ReplicaSystem}, i::Integer) = mass(s.atom
 
 function Base.getindex(sys::Union{System, ReplicaSystem}, i, x::Symbol)
     atomsbase_keys = (:position, :velocity, :atomic_mass, :atomic_number)
-    custom_keys = (:charge,)
     if hasatomkey(sys, x)
         if x in atomsbase_keys
             return getproperty(AtomsBase, x)(sys, i)
-        elseif x in custom_keys
-            return getproperty(Molly, x)(sys, i)
+        elseif x == :charge
+            return getproperty(Molly, x)(sys, i) * u"e_au"
         end
     else
         throw(KeyError("key $x not present in the system"))
@@ -1091,7 +1090,7 @@ function System(sys::AbstractSystem{D}, energy_units, force_units) where D
     for (i, atom) in enumerate(sys)
         atoms[i] = Atom(
             index=i,
-            charge=get(atom, :charge, 0.0),
+            charge=ustrip(get(atom, :charge, 0.0)), # Remove e unit
             mass=atomic_mass(atom),
             σ=(0.0 * length_unit),
             ϵ=(0.0 * energy_units),