From e3fa8b8c25e2828520b673928dfda35e74f75db6 Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Thu, 4 Jan 2024 10:22:34 +0100 Subject: [PATCH] Initial commit: Add functions for Diatomic molecules #7 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Includes functions from DiatomicBasics.jl in a submodule "Structure" - Interaction with NQCDynamics.jl Simulation objects in an extension. - Moved ase interface to an extension. - Version bump to 0.3.0 since extensions require Julia ≥1.9 --- Project.toml | 8 +- ext/NQCDynamicsBindings.jl | 62 ++++++++++ src/io/ase.jl => ext/aseInterfaceExt.jl | 6 +- src/NQCBase.jl | 7 +- src/structures.jl | 145 ++++++++++++++++++++++++ 5 files changed, 221 insertions(+), 7 deletions(-) create mode 100644 ext/NQCDynamicsBindings.jl rename src/io/ase.jl => ext/aseInterfaceExt.jl (95%) create mode 100644 src/structures.jl diff --git a/Project.toml b/Project.toml index c3d0aa3..bbb7acb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NQCBase" uuid = "78c76ebc-5665-4934-b512-82d81b5cbfb7" authors = ["James Gardner and contributors"] -version = "0.2.1" +version = "0.3.0" [deps] AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a" @@ -24,7 +24,11 @@ Requires = "1" StaticArraysCore = "1" Unitful = "1" UnitfulAtomic = "1" -julia = "1" +julia = "≥1.9" + +[extensions] +aseInterfaceExt = "PyCall" +NQCDynamicsBindings = "NQCDynamics" [extras] AtomsBaseTesting = "ed7c10db-df7e-4efa-a7be-4f4190f7f227" diff --git a/ext/NQCDynamicsBindings.jl b/ext/NQCDynamicsBindings.jl new file mode 100644 index 0000000..175f8e8 --- /dev/null +++ b/ext/NQCDynamicsBindings.jl @@ -0,0 +1,62 @@ +module NQCDynamicsBindings + +using NQCBase,NQCDynamics + +""" + distance(config::Any, i1, i2) + +Interatomic distance in Angstrom for DynamicsVariables. +""" +function NQCBase.Structure.distance(config::Any, i1::int_or_index, i2::int_or_index) + return NQCBase.Structure.distance(pos,i1,i2) +end + +""" + NQCBase.Structure.fractional_mass(sim::NQCDynamics.AbstractSimulation, index1::Int, index2::Int) + +Returns m1/(m1+m2) and m2/(m1+m2) as a vector. +""" +function NQCBase.Structure.fractional_mass(sim::NQCDynamics.AbstractSimulation, index1::Int, index2::Int) + return NQCBase.Structure.fractional_mass(sim.atoms,index1,index2) +end + +function NQCBase.Structure.reduced_mass(sim::NQCDynamics.AbstractSimulation, index1::Int, index2::Int) + return NQCBase.Structure.reduced_mass(sim.atoms,index1,index2) +end + +function NQCBase.Structure.minimum_distance_translation(config::Matrix, ind1::Int, ind2::Int, simulation::NQCDynamics.AbstractSimulation;cutoff::Int=50) + return NQCBase.Structure.minimum_distance_translation(config,ind1,ind2,simulation.cell;cutoff=cutoff) +end + +function NQCBase.Structure.minimum_distance_translation(config::Any, ind1::Int, ind2::Int, simulation::NQCDynamics.AbstractSimulation;cutoff::Int=50) + return NQCBase.Structure.minimum_distance_translation(NQCDynamics.get_positions(config),ind1,ind2,simulation.cell;cutoff=cutoff) +end + +function NQCBase.Structure.pbc_distance(config::Matrix, ind1, ind2, sim::NQCDynamics.AbstractSimulation; args...) + return NQCBase.Structure.pbc_distance(config,ind1,ind2,sim.cell; args...) +end + +function NQCBase.Structure.pbc_distance(config::Any, ind1, ind2, sim::NQCDynamics.AbstractSimulation; args...) + return NQCBase.Structure.pbc_distance(NQCDynamics.get_positions(config),ind1,ind2,sim.cell; args...) +end + +function NQCBase.Structure.pbc_center_of_mass(config::Matrix, ind1, ind2, sim::NQCDynamics.AbstractSimulation; args...) + return NQCBase.Structure.pbc_center_of_mass(config,ind1,ind2,sim.cell, sim.atoms; args...) +end + +function NQCBase.Structure.pbc_center_of_mass(config::Any, ind1, ind2, sim::NQCDynamics.AbstractSimulation; args...) + return NQCBase.Structure.pbc_center_of_mass(NQCDynamics.get_positions(config),ind1,ind2,sim.cell, sim.atoms; args...) +end + +function NQCBase.Structure.velocity_center_of_mass(config::Matrix, ind1, ind2, sim::NQCDynamics.AbstractSimulation) + return NQCBase.Structure.velocity_center_of_mass(config,ind1,ind2, sim.atoms) +end + +function NQCBase.Structure.velocity_center_of_mass(config::Any, ind1, ind2, sim::NQCDynamics.AbstractSimulation) + return NQCBase.Structure.velocity_center_of_mass(NQCDynamics.get_velocities(config),ind1,ind2, sim.atoms) +end + + + + +end \ No newline at end of file diff --git a/src/io/ase.jl b/ext/aseInterfaceExt.jl similarity index 95% rename from src/io/ase.jl rename to ext/aseInterfaceExt.jl index d03fcdb..628c1b4 100644 --- a/src/io/ase.jl +++ b/ext/aseInterfaceExt.jl @@ -1,5 +1,7 @@ +module aseInterfaceExt -using .PyCall +using NQCBase +using PyCall using Unitful, UnitfulAtomic export convert_from_ase_atoms @@ -40,3 +42,5 @@ function Cell(ase_atoms::PyObject) return PeriodicCell{Float64}(austrip.(ase_atoms.cell.array'u"Å"), ase_atoms.pbc) end end + +end \ No newline at end of file diff --git a/src/NQCBase.jl b/src/NQCBase.jl index efb58a5..14a740d 100644 --- a/src/NQCBase.jl +++ b/src/NQCBase.jl @@ -3,17 +3,16 @@ module NQCBase using PeriodicTable using Unitful, UnitfulAtomic using Requires +using Reexport : @reexport include("unit_conversions.jl") include("atoms.jl") include("cells.jl") include("io/extxyz.jl") -function __init__() - @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" @eval include("io/ase.jl") -end - include("atoms_base.jl") +include("structures.jl") +export Structure export Cell export System export Trajectory diff --git a/src/structures.jl b/src/structures.jl new file mode 100644 index 0000000..ff3ed1c --- /dev/null +++ b/src/structures.jl @@ -0,0 +1,145 @@ +module Structure + +using NQCBase: Atoms, InfiniteCell, PeriodicCell, au_to_ang +using LinearAlgebra +using Unitful +using UnitfulAtomic + +int_or_index=Union{Int, CartesianIndex} + +""" + distance(config::Matrix, i1, i2) + Interatomic distance in Angstrom for a position Matrix. +TBW +""" +function distance(config::Matrix, i1::int_or_index, i2::int_or_index) + return au_to_ang(norm(config[:,i1].-config[:,i2]))*u"Å" +end + +function distance(v1::Matrix, v2::Matrix) + return au_to_ang(norm(v1-v2))*u"Å" +end + +""" + fractional_mass(sim::Simulation, index1::Int, index2::Int) + +Returns m1/(m1+m2) and m2/(m1+m2) as a vector. +""" +function fractional_mass(atoms::Atoms, index1::int_or_index, index2::int_or_index) + return atoms.masses[[index1, index2]]./sum(atoms.masses[[index1, index2]]) +end + + +""" + reduced_mass(sim::Simulation, index1::Int, index2::Int) + +Returns the reduced mass of the diatomic in atomic unit. +""" +function reduced_mass(atoms::Atoms, index1::Int, index2::Int) + return (atoms.masses[index1]*atoms.masses[index2])/sum(atoms.masses[[index1, index2]]) +end + +#ToDo NQCDynamics aware version Simulation -> Atoms + + +#ToDo Variant with Simulation -> PeriodicCell +""" + minimum_distance_translation(config::Matrix, ind1::Int, ind2::Int, simulation::NQCDynamics.AbstractSimulation) + +Outputs a translation vector to move config[:,ind2] such that the closest distance between `ind1` and `ind2` is reached. +**The search of neighbouring unit cells will expand until the cutoff. If configurations are already subject to a minimum image convention, `cutoff=1` reduces unnecessary overhead. ** +""" +function minimum_distance_translation(config::Matrix, ind1::Int, ind2::Int, cell::PeriodicCell;cutoff::Int=50) + @debug "Cutoff was set to $(cutoff)" + if cutoff==0 || isa(cell, InfiniteCell) + return [0.0,0.0,0.0] + end + translations=[] + min_distance=[] + distance_converged=false + search=1 + while !(distance_converged || search==cutoff+1) + translations=Iterators.map(x->[x...],Iterators.product(-search:search, -search:search, -search:search)) + distances=map(x->norm(config[:,ind2].-config[:,ind1]+cell.vectors*x), translations) + push!(min_distance, (min(distances...), argmin(distances))) + distance_converged=length(min_distance)>1 ? min_distance[end][1]≈min_distance[end-1][1] : false + search+=1 + end + @debug "PBC distance check ended after $(search-1) iterations." + if search==cutoff+1 + @debug "PBC distance check terminated due to cutoff. If you aren't PBC-wrapping positions, check your cutoff is sufficiently high. " + end + @debug "Translation of choice is $(collect(translations)[min_distance[end][2]])" + return cell.vectors*collect(translations)[min_distance[end][2]] +end + +function minimum_distance_translation(config::Matrix, ind1::Int, ind2::Int, cell::InfiniteCell;cutoff::Int=50) + config +end + + +""" + angle_between(v1::Vector, v2::Vector) + +Returns the angle between two vectors **in º** based on the scalar product definition. +""" +function angle_between(v1::Vector, v2::Vector) + return rad2deg(acos(dot(v1,v2)/prod(norm, [v1, v2]))) +end + +#ToDo Versions for Any --> Matrix +#ToDo Versions for Simulation --> Cell, Masses + +""" + pbc_distance(config::Matrix, ind1::Int, ind2::Int, simulation::NQCDynamics.AbstractSimulation) + +**Returns in Angstom, not in Bohr - Check units.** + +Calculates the distance between two atoms, including a check if the copy of the second atom in any neighbouring unit cell is closer. +This variant is designed for trajectories where cell boundary wrapping has been used. +""" +function pbc_distance(config::Matrix, ind1::Int, ind2::Int, cell::PeriodicCell;args...) + return auconvert(u"Å",norm(config[:, ind1]-config[:,ind2]-minimum_distance_translation(config, ind1, ind2, cell;args...))) +end + +""" + pbc_center_of_mass(config::Matrix, ind1::Int, ind2::Int, cell::PeriodicCell, atoms::Atoms;args...) + +Generates center of mass coordinates for two atoms, including a check if the copy of the second atom in any neighbouring unit cell is closer. +""" +function pbc_center_of_mass(config::Matrix, ind1::Int, ind2::Int, cell::PeriodicCell, atoms::Atoms;args...) + pbc_translation=minimum_distance_translation(config, ind1, ind2, cell;args...) + return (atoms.masses[ind1]*config[:,ind1]+atoms.masses[ind2]*(config[:,ind2]+pbc_translation))/sum(atoms.masses[[ind1,ind2]]) +end + +""" + center_of_mass(config::Matrix, ind1::Int, ind2::Int, atoms::Atoms) + +Generates center of mass coordinates for two atoms. +""" +function center_of_mass(config::Matrix, ind1::Int, ind2::Int, atoms::Atoms) + return (atoms.masses[ind1]*config[:,ind1]+atoms.masses[ind2]*config[:,ind2])/sum(atoms.masses[[ind1,ind2]]) +end + +""" + velocity_center_of_mass(config::Matrix, ind1::Int, ind2::Int, simulation::NQCDynamics.AbstractSimulation) + + `sum(m_i*v_i)/sum(m_i)` +TBW +""" +function velocity_center_of_mass(config::Matrix, ind1::Int, ind2::Int, atoms::Atoms) + return center_of_mass(config, ind1, ind2, atoms) +end + +struct OutputSubsetKineticEnergy + indices::Vector{Int} +end +function (kinetic_energy::OutputSubsetKineticEnergy)(sol, i) + return map(x->DynamicsUtils.classical_kinetic_energy(sol.prob.p.atoms.masses[kinetic_energy.indices], x), [DynamicsUtils.get_velocities(i) for i in sol.u]) +end + + +export minimum_distance_translation, pbc_distance, distance, angle_between, pbc_center_of_mass, velocity_center_of_mass, center_of_mass, OutputSubsetKineticEnergy, reduced_mass, fractional_mass + + +end \ No newline at end of file