diff --git a/docs/make.jl b/docs/make.jl index 9163d85..df50516 100644 --- a/docs/make.jl +++ b/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"] diff --git a/docs/src/structure_creation.md b/docs/src/structure_creation.md new file mode 100644 index 0000000..c111a35 --- /dev/null +++ b/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"Å")] ) +``` + diff --git a/docs/src/use.md b/docs/src/use.md index eb0f84e..41ee9dc 100644 --- a/docs/src/use.md +++ b/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" diff --git a/src/PotentialCalculation.jl b/src/PotentialCalculation.jl index 40d0456..4ff9e47 100644 --- a/src/PotentialCalculation.jl +++ b/src/PotentialCalculation.jl @@ -2,6 +2,7 @@ module PotentialCalculation using Reexport +@reexport using AtomsBase @reexport using Unitful @reexport using UnitfulAtomic diff --git a/src/SubModules/clusters.jl b/src/SubModules/clusters.jl index 0f053a8..4731541 100644 --- a/src/SubModules/clusters.jl +++ b/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,7 +323,9 @@ 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) @@ -351,13 +333,27 @@ function Cluster(sys::AtomsBase.FlexibleSystem) 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 diff --git a/src/SubModules/restarttools.jl b/src/SubModules/restarttools.jl index dff8a12..286e4de 100644 --- a/src/SubModules/restarttools.jl +++ b/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, diff --git a/src/SubModules/sample.jl b/src/SubModules/sample.jl index 4928144..ed3380d 100644 --- a/src/SubModules/sample.jl +++ b/src/SubModules/sample.jl @@ -1,5 +1,6 @@ module Sample +using AtomsBase using LinearAlgebra using ProgressMeter using Rotations @@ -46,7 +47,7 @@ end """ - line_sampler(cl1::Cluster, cl2::Cluster; npoints=10, mindis=2.0, maxdis=9.0) + line_sampler(cl1::AbstractSystem, cl2::AbstractSystem; npoints=10, mindis=2.0, maxdis=9.0) Samples cluster position in a line like fashion with even intervals. @@ -57,9 +58,9 @@ Sampling is done by first generating random direction and placing `cl2` in that direction from `cl2` at `mindis` and then picking `npoints` in even distances to `maxdis`. """ -function line_sampler(cl1::Cluster, cl2::Cluster; npoints=10, mindis=2.0u"Å", maxdis=9.0u"Å") - c1 = deepcopy(cl1) - c2 = deepcopy(cl2) +function line_sampler(cl1::AbstractSystem, cl2::AbstractSystem; npoints=10, mindis=2.0u"Å", maxdis=9.0u"Å") + c1 = Cluster(cl1) + c2 = Cluster(cl2) center_cluster!(c1) center_cluster!(c2) @@ -78,7 +79,7 @@ end """ - adaptive_line_sampler(cal::Calculator, cl1::Cluster, cl2::Cluster, max_e=0; unit="cm-1", + adaptive_line_sampler(cal::Calculator, cl1::AbstractSystem, cl2::AbstractSystem, max_e=0; unit="cm-1", npoints=10, maxdis=9.0, sstep=0.1, startdistance=3.0, pchannel=undef) Calculates potential on a line type distances that is suitable for visualization. @@ -97,14 +98,14 @@ has energy lower than `max_e` - `id=""` : extra identificaltion to help to do parallel computations - `pchannel=undef` : (Remote)Channel where progess information is added """ -function adaptive_line_sampler(cal::Calculator, cl1::Cluster, cl2::Cluster, max_e=0u"cm^-1"; +function adaptive_line_sampler(cal::Calculator, cl1::AbstractSystem, cl2::AbstractSystem, max_e=0u"cm^-1"; npoints=10, maxdis=9.0u"Å", sstep=0.1, startdistance=3.0u"Å", basename="base", id="", pchannel=undef) @debug "Starting adaptive_line_sampler" # take deepcopys so that originals are not moved - c1 = deepcopy(cl1) - c2 = deepcopy(cl2) + c1 = Cluster(cl1) + c2 = Cluster(cl2) center_cluster!(c1) center_cluster!(c2) @@ -191,11 +192,11 @@ end Structure used to help use of adaptive line samplers # Arguments -- `cal::Calculator` : [`calculator`](@ref) used in calculations -- `cl1::Cluster` : first cluster -- `cl2::Cluster` : second cluster -- `nlines` : number of lines calculated -- `max_e` : maximum energy for configuration +- `cal::Calculator` : [`calculator`](@ref) used in calculations +- `cl1::AbstractSystem` : first cluster +- `cl2::AbstractSystem` : second cluster +- `nlines` : number of lines calculated +- `max_e` : maximum energy for configuration # Keywords - `npoints` : number of points in line @@ -218,17 +219,17 @@ mutable struct InputAdaptiveSampler maxdis sstep startdistance - function InputAdaptiveSampler(cal::Calculator, cl1::Cluster, cl2::Cluster, + function InputAdaptiveSampler(cal::Calculator, cl1::AbstractSystem, cl2::AbstractSystem, nlines, max_e=0u"cm^-1"; npoints=10, maxdis=9.0u"Å", sstep=0.1u"Å", startdistance=3.0u"Å") - new(cal, cl1, cl2, nlines, max_e,npoints,maxdis,sstep,startdistance) + new(cal, Cluster(cl1), Cluster(cl2), nlines, max_e,npoints,maxdis,sstep,startdistance) end end """ -sample_multiple_adaptive_lines(cal::Calculator, cl1::Cluster, cl2::Cluster, nlines, max_e=0; +sample_multiple_adaptive_lines(cal::Calculator, cl1::AbstractSystem, cl2::AbstractSystem, nlines, max_e=0; unit="cm-1", npoints=10, maxdis=9.0, sstep=0.1, startdistance=3.0, basename="base", id="") @@ -237,18 +238,20 @@ Distances are sampled in line like fashion in random directions with even distan Sampling of line is done with [`adaptive_line_sampler`](@ref). # Arguments -- `cal::Calculator` : [`calculator`](@ref) used in calculations -- `cl1::Cluster` : first cluster -- `cl2::Cluster` : second cluster -- `nlines` : number of lines calculated -- `max_e=0` : maximum energy for configuration -- `npoints=10` : number of points in line -- `maxdis=9.0` : maximum cluster distance -- `sstep=0.1` : search step used by [`adaptive_line_sampler`](@ref) -- `startdistance=3.0` : distance where search is started -- `basename="base"` : prefix for temporary files in calculations -- `id=""` : extra identificaltion to help to do parallel computations -- `pchannel=undef` : (Remote)Channel where progess information is added +- `cal::Calculator` : [`calculator`](@ref) used in calculations +- `cl1::AbstractSystem` : first cluster +- `cl2::AbstractSystem` : second cluster +- `nlines` : number of lines calculated +- `max_e=0` : maximum energy for configuration + +# Keywords +- `npoints=10` : number of points in line +- `maxdis=9.0` : maximum cluster distance +- `sstep=0.1` : search step used by [`adaptive_line_sampler`](@ref) +- `startdistance=3.0` : distance where search is started +- `basename="base"` : prefix for temporary files in calculations +- `id=""` : extra identificaltion to help to do parallel computations +- `pchannel=undef` : (Remote)Channel where progess information is added # Returns [`Dict`](@ref) with keys @@ -256,11 +259,11 @@ Sampling of line is done with [`adaptive_line_sampler`](@ref). * "Points" : array of [`Cluster`](@ref) representing points where energy was calculated * "Mindis" : array of minimum distances in each line """ -function sample_multiple_adaptive_lines(cal::Calculator, cl1::Cluster, cl2::Cluster, nlines, max_e=0u"cm^-1"; +function sample_multiple_adaptive_lines(cal::Calculator, cl1::AbstractSystem, cl2::AbstractSystem, nlines, max_e=0u"cm^-1"; npoints=10, maxdis=9.0u"Å", sstep=0.1u"Å", startdistance=3.0, basename="base", id="", pchannel=undef) - c1 = deepcopy(cl1) - c2 = deepcopy(cl2) + c1 = Cluster(cl1) + c2 = Cluster(cl2) out = Vector{Dict}(undef, nlines) rtmp = Float64[] sr = startdistance diff --git a/test/test_clusters.jl b/test/test_clusters.jl index f28a017..d209a2b 100644 --- a/test/test_clusters.jl +++ b/test/test_clusters.jl @@ -15,52 +15,35 @@ using AtomsBase -7.2001330967 0.3718768293 -0.0703451879]', AtomOnlySymbol.(["C", "O", "O", "H", "H"]) ) - Ar = Cluster(rand(3), AtomOnlySymbol("Ar")) - - nf = ClusterNoSymbols(rand(3,5)) - nar = ClusterNoSymbols(rand(3)) - - @test length(nf) == length(formic_acid) + Ar = Cluster( Atom(:Ar, rand(3)u"Å" ) ) t1 = formic_acid[5] t2 = formic_acid[end] - @test t1.atoms == t2.atoms && t1.xyz == t2.xyz - - tn = nf[3] + @test position(t1) == position(t2) - @test all(tn.xyz .== nf.xyz[:,3]) + @test length(formic_acid(1:2)) == 2 - @test length(formic_acid[1:2]) == 2 - @test length(nf[2:4]) == 3 - - @test length(nf+nar) == 6 @test length(formic_acid + Ar) == 6 distances(formic_acid, Ar) - distances(formic_acid, nf) distances(formic_acid, 1:2, 3:4) - distances(formic_acid, 1, nf, 2) - distances(formic_acid, 1:2, nf, 2:4) + distances(formic_acid, 1:2, formic_acid, 2:4) - show(devnull, nf) print(devnull, formic_acid) - move!(nf, [10,0,0]) - - for x in (formic_acid, nf) - cc = center_coordinates(x) - center_cluster!(x) - @test all(abs.(center_coordinates(x)) .< 1E-12u"Å" .* ones(size(cc))) - end + + cc = center_coordinates(formic_acid) + center_cluster!(formic_acid) + @test all(abs.(center_coordinates(formic_acid)) .< 1E-12u"Å" .* ones(size(cc))) l = distances(formic_acid, 2,4) a = cluster_angle(formic_acid, 1,2,4) d = dihedral_angle(formic_acid, 3,1,2,4) - cluster_angle(formic_acid, 1,2, nf, 3) + cluster_angle(formic_acid, 1,2, formic_acid, 3) @test l != distances(formic_acid, 1,2) @test a != cluster_angle(formic_acid, 1,2,3) @@ -74,17 +57,6 @@ using AtomsBase @test l ≈ tl && a ≈ ta && d ≈ td end - distances(nf, 2,4) - cluster_angle(nf, 1,2,4) - dihedral_angle(nf, 3,1,2,4) - - - - @test_throws DimensionMismatch ClusterNoSymbols(rand(4)) - @test_throws DimensionMismatch ClusterNoSymbols(rand(2)) - @test_throws DimensionMismatch ClusterNoSymbols(rand(2,3)) - @test_throws DimensionMismatch ClusterNoSymbols(rand(4,2)) - @test_throws DimensionMismatch Cluster(rand(2), [AtomOnlySymbol("H")]) @test_throws DimensionMismatch Cluster(rand(2), AtomOnlySymbol("H")) @test_throws DimensionMismatch Cluster(rand(3),AtomOnlySymbol.(["H", "O"])) @@ -92,10 +64,10 @@ using AtomsBase @test_throws DimensionMismatch Cluster(rand(4,2),AtomOnlySymbol.(["H", "O"])) # AtomsBase - fa = FlexibleSystem(formic_acid) - pf = FlexibleSystem(formic_acid; e=1) - @test haskey(pf, :e) + fa = FlexibleSystem(formic_acid; e=1) + pf = Cluster(collect(formic_acid)) + @test haskey(fa, :e) cfa = Cluster(fa) - @test cfa.atoms == formic_acid.atoms - @test cfa.xyz == formic_acid.xyz + @test all( atomic_symbol(cfa) .== atomic_symbol(formic_acid) ) + @test all( position(cfa) .== position(formic_acid) ) end diff --git a/test/test_restarttools.jl b/test/test_restarttools.jl index fbdc041..6112d00 100644 --- a/test/test_restarttools.jl +++ b/test/test_restarttools.jl @@ -20,8 +20,8 @@ formic_acid=Cluster( -7.2001330967 0.3718768293 -0.0703451879]', AtomOnlySymbol.(["C", "O", "O", "H", "H"]) ) -Ar = Cluster(rand(3), AtomOnlySymbol.(["Ar"])) -N2 = Cluster([[0.0 1.0]; [0.0 0.0]; [0.0 0.0]], AtomOnlySymbol.(["N", "N"])) +Ar = Cluster( Atom(:Ar, rand(3)u"Å" ) ) +N2 = isolated_system( [Atom(:N, [1., 0., 0.].*u"Å"), Atom(:N, [0., 0., 0.].*u"Å")] ) open(xyzname,"w") do io print_xyz(io, formic_acid) @@ -31,7 +31,7 @@ pbar=true testrestarts = false -if Sys.which("orca") != nothing +if Sys.which("orca") != nothing && Sys.which("orca_scf") != nothing @info "Orca binary found. Testing ORCA." @testset "Orca" begin ca = Calculator("blyp d3bj TIGHTSCF", "def2-svp", Orca())