Complete AtomsBase support

docs/make.jl

@@ -4,6 +4,7 @@ using PotentialCalculation, PotentialCalculation.psi4
 makedocs(sitename="PotentialCalculation.jl",
          pages=["Home" => "index.md",
                 "Install" => "install.md",
+                "Building Molecules" => "structure_creation.md"
                 "Usage" => "use.md",
                 "References" => "reference.md"]
 

docs/src/structure_creation.md

@@ -0,0 +1,33 @@
+# Building Molecules
+
+Input structures are defined by [AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl),
+but internally calculations are done with a structure called [`Cluster`](@ref).
+It is [AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl) compatible and
+you can convert back and forth and use it where AtomsBase structures can be used.
+Atoms base is exported by PotentialCalculation, so there is no need to load it.
+
+
+While you can create [`Cluster`](@ref) directly, it is recommended that you use
+[AtomsBase](https://github.com/JuliaMolSim/AtomsBase.jl) to create structures.
+You can e.g. use [AtomsIO](https://github.com/mfherbst/AtomsIO.jl)
+to load data molecule or just create molecules by hand on the fly.
+
+## Example
+
+```julia
+using PotentialCalculation
+using AtomsIO
+
+# Load from file with AtomsIO
+system = load_system("molecule.xyz")
+
+# Convert to Cluster
+csystem = Cluster(system)
+
+# Convert back to AtomsBase standard structure
+fsystem = FlexibleSystem(csystem)
+
+# Create structure by hand
+N2 = isolated_system( [Atom(:N, [1., 0., 0.].*u"Å"), Atom(:N, [0., 0., 0.].*u"Å")] )
+```
+

docs/src/use.md

@@ -30,7 +30,10 @@ mp2 = Calculator(
       )
 
 # Creating argon atom
-Ar = Cluster(zeros(3), AtomOnlySymbol("Ar"))
+Ar = Cluster( Atom(:Ar, zeros(3)u"Å")  )
+# Also
+# N2 = Cluster( Atom(:N, [1., 0., 0.].*u"Å"), Atom(:N, [0., 0., 0.].*u"Å") )
+# Or any other AtomsBase compatible struct can be used
 
 # File where other molecule is taken
 trajfile="Some trajectory.xyz"

src/PotentialCalculation.jl

@@ -2,6 +2,7 @@ module PotentialCalculation
 
 using Reexport
 
+@reexport using AtomsBase
 @reexport using Unitful
 @reexport using UnitfulAtomic
 

src/SubModules/clusters.jl

@@ -6,7 +6,6 @@ export center_cluster!
 export center_coordinates
 export Cluster
 export cluster_angle
-export ClusterNoSymbols
 export dihedral_angle
 export distances
 export move!
@@ -30,33 +29,9 @@ import Base.==, Base.+
 
 
 
-abstract type AbstractCluster end
+abstract type AbstractCluster <: AbstractSystem{3} end
 abstract type AbstractClusterWithSymbols{T} <: AbstractCluster where T<:AbstractAtom end
 
-"""
-    ClusterNoSymbols <: AbstractCluster
-
-Basic cluster has only location of atoms.
-
-# Fields
-- `xyz::Array{Float64,2}` : location of atoms in 3d space, first index is x, y, z coordinate
-"""
-mutable struct ClusterNoSymbols <: AbstractCluster
-    xyz::Array{Float64,2}
-    function ClusterNoSymbols(xyz::AbstractArray{<:Real,2})
-        if size(xyz,1) != 3
-            throw(DimensionMismatch("ClusterNoSymbols - xyz has wrong dimensions"))
-        end
-        new(xyz)
-    end
-    function ClusterNoSymbols(xyz::AbstractArray{<:Real,1})
-        if length(xyz) != 3
-            throw(DimensionMismatch("ClusterNoSymbols - atoms has wrong dimensions length=$(size(xyz))"))
-        end
-        new(reshape(xyz,3,1))
-    end
-end
-
 
 """
     Cluster{T} <: AbstractClusterWithSymbols{T} where T<:AbstractAtom
@@ -111,26 +86,23 @@ function (+)(c1::Cluster{T}, c2::Cluster{T}) where T
     return Cluster(hcat(c1.xyz,c2.xyz),vcat(c1.atoms,c2.atoms))
 end
 
-function (+)(c1::ClusterNoSymbols, c2::ClusterNoSymbols)
-    return ClusterNoSymbols(hcat(c1.xyz,c2.xyz))
-end
-
 
-function Base.getindex(C::Cluster, i::Int)
+function (C::Cluster)(i::Int)
     return Cluster(C.xyz[:,i], [C.atoms[i]])
 end
 
-function Base.getindex(C::ClusterNoSymbols, i::Int)
-    return ClusterNoSymbols(C.xyz[:,i])
+
+function Base.getindex(c::Cluster, i::Int)
+    Atom( atomic_symbol(c,i), position(c,i) )
 end
 
 
-function Base.getindex(C::Cluster, i::AbstractUnitRange)
-    return Cluster(C.xyz[:,i], C.atoms[i])
+function Base.getindex(c::Cluster, s::Symbol)
+    [ x for x in c if atomic_symbol(x) == s ]
 end
 
-function Base.getindex(C::ClusterNoSymbols, i::AbstractUnitRange)
-    return ClusterNoSymbols(C.xyz[:,i])
+function (C::Cluster)(i::AbstractUnitRange)
+    return Cluster(C.xyz[:,i], C.atoms[i])
 end
 
 
@@ -206,7 +178,7 @@ Unitful.uconvert(::typeof(u"Å"), r::Real) = r*u"Å"
 Moves cluster by `r`
 """
 function move!(c::AbstractCluster,r)
-    tmp = @. uconvert(u"Å", r) |> ustrip
+    tmp = ustrip.(u"Å", r)
     c.xyz .+= tmp
 end
 
@@ -330,6 +302,14 @@ end
 
 ## AtomsBase support
 
+
+AtomsBase.species_type(::Cluster) = Atom
+
+AtomsBase.atomkeys(::Cluster) = (:atomic_symbol, :position)
+
+Base.keys(::Cluster) = ()
+
+
 function AtomsBase.bounding_box(::AbstractCluster)
     a = SVector{3}( [Inf, 0., 0.] .* u"bohr" )
     b = SVector{3}( [0., Inf, 0.] .* u"bohr" )
@@ -343,21 +323,37 @@ function AtomsBase.boundary_conditions(::AbstractCluster)
 end
 
 
-function Cluster(sys::AtomsBase.FlexibleSystem)
+function Cluster(sys::AbstractSystem)
+    @assert :atomic_symbol in atomkeys(sys)
+    @assert :position in atomkeys(sys)
     a = AtomOnlySymbol.( atomic_symbol(sys) )
     pos = map( position(sys) ) do r
         ustrip.(u"Å", r)
     end
     return Cluster(hcat(pos...), a)
 end
 
+function Cluster(aarray::Atom...)
+    xyz = ustrip.(u"Å", hcat( position.(aarray)... ) )
+    a = (AtomOnlySymbol ∘ String ∘ atomic_symbol).( collect(aarray)  )
+    return Cluster(xyz, a)
+end
+
+Cluster(aarray::AbstractArray{Atom}) = Cluster(aarray...)
+
 
 function AtomsBase.FlexibleSystem(c::Cluster; kwargs...)
-    a = [
-        Symbol(c.atoms[i].id) => SVector{3}( c.xyz[:,i] .*u"Å")
-        for i in 1:length(c)
-    ]
-    return isolated_system(a; kwargs...)
+    return isolated_system(collect(c); kwargs...)
+end
+
+
+function AtomsBase.position(c::Cluster, i::Int)
+    return SVector( c.xyz[:,i]... ) .* u"Å"
+end
+
+function AtomsBase.atomic_symbol(c::Cluster,i::Int)
+    return Symbol(c.atoms[i].id)
 end
 
+
 end #module Clusters

src/SubModules/restarttools.jl

@@ -24,6 +24,7 @@ using ..Calculators
 using ..Clusters
 using ..Fileaccess
 using ..Sample
+using AtomsBase
 using Distributed
 using ProgressMeter
 using Unitful
@@ -172,7 +173,7 @@ function continue_calculation(
     iscalculated = deepcopy(iscal)
     @async for i in eachindex(data["Points"])
         if ! iscal[i]
-            x = (calculator,data["Points"][i][rc1], data["Points"][i][rc2])
+            x = ( calculator,data["Points"][i](rc1), data["Points"][i](rc2) )
             put!(jobs, (i,x) )
         end
     end
@@ -249,14 +250,21 @@ function calculate_potential(fname::AbstractString, calculator::Calculator;
     data = load_jld_data(fname)
     @info "File $(fname) loaded - calculating with different method"
     flush(stdout)
+
+    # Split data to different clusters for the calculation
     lc1 = length(data["cluster1"])
-    c1_points = map( x-> x[1:lc1], data["Points"])
-    c2_points = map( x-> x[lc1+1:end], data["Points"])
+    lc2 = length(data["cluster2"])
+    rc1 =  1:lc1           #UnitRange for cluster1
+    rc2 =  lc1+1:lc1+lc2   #UnitRange for cluster2
+    c1_points = map( c-> c(rc1), data["Points"])
+    c2_points = map( c-> c(rc2), data["Points"])
 
+    # Set progressbar max value
     if pbar
         prog = Progress(length(c1_points) ,dt=1, desc="Calculating points:")
     end
 
+    # Create channels and schedule calculation
     jobs = RemoteChannel(()->Channel(2*nworkers()))
     results = RemoteChannel(()->Channel(2*nworkers()))
 
@@ -269,6 +277,7 @@ function calculate_potential(fname::AbstractString, calculator::Calculator;
         remote_do(parallel_work, p, jobs, results, _calculate_points)
     end
 
+    # Collect results
     energy = similar(c1_points, Float64)
     iscalculated = falses(size(energy))
     n = 0
@@ -501,24 +510,7 @@ function calculate_adaptive_sample_inputs(inputs::AbstractArray{InputAdaptiveSam
  end
 
 
-"""
-    create_inputs(cluster1, cluster2, calculator::Calculator;  kwargs...)
-
-
-
-Creates inputs for the calculations. Use [`calculate_adaptive_sample_inputs`](@ref) on
-results to do the calculation itself.
-
-Inputs for the molecules (named cluster1/2) can be loaded from xyz-file or given by as
-[`Cluster`](@ref) structures. If given as file names, then `nsamples` points are randomly
-picked from the trajectories for calculation. Total number of lines sampled is thus
-`nlines` times `nsamples`.
-
-Parallelisation is done over `nsamples`, thus it is recommended for it to be multiple of
-number workers processes.
-
-# Arguments
-- `cluster1` : something that can be interpreted as [`Cluster`](@ref) or an Array of it
+"""AbstractClusterthing that can be interpreted as [`Cluster`](@ref) or an Array of it
 - `cluster2` : something that can be interpreted as [`Cluster`](@ref) or an Array of it
 - `calculator::Calculator` : [`Calculator`](@ref) used in sampling
 
@@ -611,8 +603,8 @@ function create_inputs(cluster1,
             )
 end
 
-function create_inputs(cluster1::AbstractCluster,
-                cluster2::AbstractCluster,
+function create_inputs(cluster1::AbstractSystem,
+                cluster2::AbstractSystem,
                 calculator::Calculator;
                 nlines=1,
                 nsamples=2,
@@ -629,8 +621,8 @@ function create_inputs(cluster1::AbstractCluster,
                                   sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
 end
 
-function create_inputs(cluster1::AbstractCluster,
-                cluster2::AbstractArray{<:AbstractCluster},
+function create_inputs(cluster1::AbstractSystem,
+                cluster2::AbstractArray{<:AbstractSystem},
                 calculator::Calculator;
                 nlines=1,
                 nsamples=2,
@@ -648,8 +640,8 @@ function create_inputs(cluster1::AbstractCluster,
 end
 
 
-function create_inputs(cluster1::AbstractArray{<:AbstractCluster},
-                cluster2::AbstractCluster,
+function create_inputs(cluster1::AbstractArray{<:AbstractSystem},
+                cluster2::AbstractSystem,
                 calculator::Calculator;
                 nlines=1,
                 nsamples=2,
@@ -666,8 +658,8 @@ function create_inputs(cluster1::AbstractArray{<:AbstractCluster},
                                   sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
 end
 
-function create_inputs(cluster1::AbstractArray{<:AbstractCluster},
-                cluster2::AbstractArray{<:AbstractCluster},
+function create_inputs(cluster1::AbstractArray{<:AbstractSystem},
+                cluster2::AbstractArray{<:AbstractSystem},
                 calculator::Calculator;
                 nlines=1,
                 nsamples=2,