diff --git a/Project.toml b/Project.toml index fc4e86d69..c6920369b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,9 @@ name = "NQCDynamics" uuid = "36248dfb-79eb-4f4d-ab9c-e29ea5f33e14" authors = ["James "] -version = "0.14.0" + +version = "0.15.0" + [deps] AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d" @@ -46,6 +48,9 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" +[weakdeps] +MACEModels = "29bafc4c-4f38-420f-a153-8a5a5e0bd5f6" + [compat] AdvancedHMC = "0.5, 0.6" AdvancedMH = "0.8" @@ -58,7 +63,6 @@ FastBroadcast = "0.1, 0.2" FastLapackInterface = "1, 2" HDF5 = "0.15, 0.16, 0.17" Interpolations = "0.13, 0.14, 0.15" -MKL_jll = "2023" MuladdMacro = "0.2" NQCBase = "0.2" NQCDistributions = "0.1" @@ -86,7 +90,7 @@ UnPack = "1" UnicodePlots = "2, 3" Unitful = "1" UnitfulAtomic = "1" -julia = "1.7" +julia = "≥1.9" [extras] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" diff --git a/docs/Project.toml b/docs/Project.toml index 1a502da20..fa07961d5 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -11,13 +11,13 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MACEModels = "29bafc4c-4f38-420f-a153-8a5a5e0bd5f6" NNInterfaces = "07253886-9aba-4a5d-b3d5-29b8b100a664" NQCBase = "78c76ebc-5665-4934-b512-82d81b5cbfb7" NQCDistributions = "657014a1-bd25-4aa2-92cb-cbcede0310ad" NQCDynamics = "36248dfb-79eb-4f4d-ab9c-e29ea5f33e14" NQCModels = "c814dc9f-a51f-4eaf-877f-82eda4edad48" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" RingPolymerArrays = "81cd853e-d334-476c-b156-f95acc8a32dc" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/docs/make.jl b/docs/make.jl index a27290e31..0616a7af8 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,6 @@ using Documenter using DocumenterCitations -using NQCBase, NQCModels, NQCDistributions, NQCDynamics +using NQCBase, NQCModels, NQCDistributions, NQCDynamics, MACEModels using CubeLDFAModel, NNInterfaces bib = CitationBibliography(joinpath(@__DIR__, "references.bib")) @@ -15,12 +15,14 @@ end @time makedocs(; plugins=[bib], sitename="NQCDynamics.jl", - modules=[NQCDynamics, NQCDistributions, NQCModels, NQCBase, CubeLDFAModel], + modules=[NQCDynamics, NQCDistributions, NQCModels, NQCBase, CubeLDFAModel, MACEModels], doctest=false, format=Documenter.HTML( prettyurls=get(ENV, "CI", nothing) == "true", canonical="https://nqcd.github.io/NQCDynamics.jl/stable/", assets=["assets/favicon.ico", "assets/citations.css"], + size_threshold = 5000*1024, + example_size_threshold = 8000*1024, ), authors="James Gardner and contributors.", pages=[ @@ -31,10 +33,11 @@ end "Saving and loading" => "saving_loading.md" "NQCModels.jl" => Any[ "NQCModels/overview.md" + "NQCModels/combining_models.md" "NQCModels/analyticmodels.md" - "NQCModels/ase.md" - "NQCModels/neuralnetworkmodels.md" + "NQCModels/fullsizemodels.md" "NQCModels/frictionmodels.md" + "NQCModels/neuralnetworkmodels.md" ] "NQCDistributions.jl" => Any[ "NQCDistributions/overview.md" diff --git a/docs/references.bib b/docs/references.bib index 5172d9fda..089e1ee55 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -361,4 +361,18 @@ @article{Korol2020 URL = {https://doi.org/10.1063/1.5134810}, eprint = {https://doi.org/10.1063/1.5134810} } - +@article{Scholz2019, + title = {Vibrational Response and Motion of Carbon Monoxide on {{Cu}}(100) Driven by Femtosecond Laser Pulses: {{Molecular}} Dynamics with Electronic Friction}, + shorttitle = {Vibrational Response and Motion of Carbon Monoxide on {{Cu}}(100) Driven by Femtosecond Laser Pulses}, + author = {Scholz, Robert and Lindner, Steven and Lon{\v c}ari{\'c}, Ivor and Tremblay, Jean Christophe and Juaristi, J. I. and Alducin, M. and Saalfrank, Peter}, + year = {2019}, + month = dec, + journal = {Phys. Rev. B}, + volume = {100}, + number = {24}, + pages = {245431}, + issn = {2469-9950, 2469-9969}, + doi = {10.1103/PhysRevB.100.245431}, + urldate = {2023-04-10}, + langid = {english}, +} diff --git a/docs/src/NQCModels/combining_models.md b/docs/src/NQCModels/combining_models.md new file mode 100644 index 000000000..58f9edb35 --- /dev/null +++ b/docs/src/NQCModels/combining_models.md @@ -0,0 +1,47 @@ +```@setup logging +@info "Expanding src/NQCModels/combining_models.md..." +start_time = time() +``` + +# Composing multiple models + +## Why? + +For some non-adiabatic dynamics methods, we would like to calculate certain quantities using a number of different approximations or different machine learning model interfaces. +However, a model for a potential energy surface does not and should not be able to provide e.g. an electronic friction tensor. +Instead, we can use the modular nature of NQCModels.jl to combine different components to a non-adiabatic dynamics method into a single model that can be used to set up a [`Simulation`](@ref) + +For example, in order to run [Molecular Dynamics with Electronic friction](@ref mdef-dynamics), we need to use a `Model` which supports the following methods: +- `potential` to determine the total energy +- `derivative` to determine interatomic forces +- `friction` to determine the electronic friction tensor + +However, most models for potential energy surfaces only provide `potential` and `derivative`, whereas `friction!` is provided by an `ElectronicFrictionProvider`. + +To combine different models together for a system, we can use [`Subsystem`](@ref)s to define which atoms in the system a given model should be applied to, and a [`CompositeModel`](@ref) to combine the different models into one. +If different [`Subsystem`](@ref)s should experience different effective temperatures, a [`Thermostat`](@ref) can be provided when initialising a [`Simulation`](@ref). + +![An overview of the different components used to combine models in NQCModels.jl](../assets/compositemodels/struct-explainer.svg) + +### Example: Combining an `AdiabaticModel` with `ElectronicFrictionProvider`s + +```@example compositemodels +using NQCDynamics + +atoms = vcat([:H, :H], [:Cu for _ in 1:54]) # create example surface structure +positions = rand(3,56) # Placeholder for atomic positions +pes_model = Subsystem(Free(3)) # PES model with potential and derivative +friction_model = Subsystem(RandomFriction(3)) # Friction model supporting friction! only + +println("PES model: \n", pes_model, "\nFriction model:", "friction_model") + +complete_model = CompositeModel(pes_model, friction_model) +``` + +As shown above, we have combined the `RandomFriction` provider with a model potential to give a total model which can be used in a simulation. + +```@example compositemodels +println("Potential: ", NQCModels.potential(complete_model, positions)) +println("Derivative: ", NQCModels.derivative(complete_model, positions)) +println("Friction: ", NQCModels.friction(complete_model, positions)) +``` diff --git a/docs/src/NQCModels/ase.md b/docs/src/NQCModels/fullsizemodels.md similarity index 54% rename from docs/src/NQCModels/ase.md rename to docs/src/NQCModels/fullsizemodels.md index 354101ed5..7f0bf8dca 100644 --- a/docs/src/NQCModels/ase.md +++ b/docs/src/NQCModels/fullsizemodels.md @@ -1,8 +1,11 @@ ```@setup logging -@info "Expanding src/NQCModels/ase.md..." +@info "Expanding src/NQCModels/fullsizemodels.md..." start_time = time() ``` -# ASE interface + +# Full dimensional model library + +## ASE interface The easiest way to obtain potentials and forces from established codes is to use the interfaces implemented in [ASE](https://wiki.fysik.dtu.dk/ase/). @@ -10,6 +13,7 @@ use the interfaces implemented in [ASE](https://wiki.fysik.dtu.dk/ase/). We provide the [`AdiabaticASEModel`](@ref) which wraps an ASE atoms object and its associated calculator to implement the required [`potential`](@ref) and [`derivative`](@ref) functions. +Several examples for connecting common machine-learning interatomic potentials to NQCModels.jl through the ASE interface are shown in the [MLIP examples](@ref) section. !!! note @@ -21,11 +25,12 @@ associated calculator to implement the required [`potential`](@ref) and pre-installed Python or install its own. Refer to the [PyCall README](https://github.com/JuliaPy/PyCall.jl) for installation and configuration instructions. - -## Example + +### Example First, it is necessary to import `ase` and create the `ase.Atoms` object and attach the desired calculator. This works exactly as in Python: + ```@example ase using PythonCall: pyimport, pylist using ASEconvert @@ -39,24 +44,63 @@ nothing # hide Next, the [`AdiabaticASEModel`](@ref) is created by passing the `ase.Atoms` object directly to the model: + ```@repl ase using NQCModels model = AdiabaticASEModel(h2) ``` + Now the model can be used in the same way as any of the previously introduced analytic models. + ```@repl ase potential(model, rand(3, 2)) derivative(model, rand(3, 2)) ``` -!!! tip +!!! tip In theory, this should work with any of the ASE calculators that correctly implement the `get_potential_energy` and `get_forces` functions. For instance, you can use [SchNetPack (SPK)](https://github.com/atomistic-machine-learning/schnetpack) by passing their ASE calculator to the `AdiabaticASEModel`. - Take a look at [Neural network models](@ref) to learn more. + Take a look at [MLIP examples](@ref) to learn more. + +## MACE interface + +[MACEModels.jl](https://github.com/NQCD/MACEModels.jl) is an interface to the [MACE](https://github.com/ACEsuit/mace) code. The package attempts to improve calculation speeds by directly interfacing with `MACE`, instead of going through `ase`. + +More information on how to use MACE models is available in the [MACEModels.jl](@ref) API documentation. + +!!! tip + + To make the installation of a python environment to execute MACE in optional, the `MACEModel` is provided in a separate package, which needs to be installed with `]add MACEModels`. + +### Example + +```julia +using NQCModels, PyCall, NQCDynamics.jl, MACEModels + +ensemble_paths = [ + "model01.model", + "model02.model", + "model03.model", +] + +structure = Atoms([:N, :H, :H, :H]) # Define the structure +cell = InfiniteCell() # Set periodicity (none in this case) + +# Create the MACE model +model = MACEModel( + atoms, # Atoms (used for neighbour lists) + cell, # Cell (used for neighbour lists) + ensemble_paths; # Vector containing paths to one or more models + device = "cpu", # Device (or vector of devices) to use + default_dtype = Float32, # Data type of the models (Float32 or Float64) + batch_size = 1, # Number of structures to evaluate at once (improves overall throughput) +) +``` + ```@setup logging runtime = round(time() - start_time; digits=2) @info "...done after $runtime s." diff --git a/docs/src/NQCModels/neuralnetworkmodels.md b/docs/src/NQCModels/neuralnetworkmodels.md index e8011c54a..1f2f0e4c7 100644 --- a/docs/src/NQCModels/neuralnetworkmodels.md +++ b/docs/src/NQCModels/neuralnetworkmodels.md @@ -2,7 +2,10 @@ @info "Expanding src/NQCModels/neuralnetworkmodels.md..." start_time = time() ``` -# Neural network models + +# MLIP Examples + +## SchNet (SchNetPack) models Using the [ASE interface](@ref) we can directly use models trained using [SchNetPack](https://github.com/atomistic-machine-learning/schnetpack). @@ -23,8 +26,9 @@ provide the relative path. First we load the model into an `ase` calculator and attach it to our diatomic hydrogen molecule. + ```julia -using PyCall +using PythonCall ase = pyimport("ase") spkutils = pyimport("schnetpack.utils") @@ -40,10 +44,11 @@ h2.set_calculator(calc) We can obtain the energies and forces from `ase` directly in the usual way, converting them to atomic units using [UnitfulAtomic](https://github.com/sostock/UnitfulAtomic.jl). + ```julia-repl using Unitful, UnitfulAtomic; -austrip(h2.get_total_energy() * u"eV") -austrip.(h2.get_forces() .* u"eV/Å") +austrip(pyconvert(Float64,h2.get_total_energy()) * u"eV") +austrip.(pyconvert(Matrix{Float64}, h2.get_forces()) .* u"eV/Å") ``` !!! warning @@ -54,6 +59,7 @@ austrip.(h2.get_forces() .* u"eV/Å") Then, we can convert the ASE output into the format used in NQCModels, which makes it possible to use the SchNet model e.g. for molecular dynamics calculations within NQCDynamics.jl: + ```julia-repl using NQCModels; model = AdiabaticASEModel(h2); @@ -63,6 +69,34 @@ r = [0 0; 0 0; 0 ustrip(auconvert(0.74u"Å"))] potential(model, r) derivative(model, r) ``` + +## MACE models + +!!! warning + + A [more direct interface](https://github.com/NQCD/MACEModels.jl) for the use of MACE models is now available. + + The following example is still valid but the new interface is recommended for new projects. ([Example here](@ref mdef-ttm-example)) + +The following example shows how to connect a trained [MACE](https://github.com/ACEsuit/mace) model to NQCDynamics.jl using the [ASE interface](@ref). + +```julia +using PythonCall, NQCModels + +ase_io = pyimport("ase.io") # make sure ase is installed in your environment first +mace_module = pyimport("mace.calculators") # make sure mace is installed in your environment first + +structure = ase_io.read("starting_structure.xyz") +mace_calculator = mace_module.MACECalculator( + "path/to/model.model", + device = "cpu", # or cuda to run on GPU + default_dtype = "float32" # if this was set during training as well +) + +structure.set_calculator(mace_calculator) +model = AdiabaticASEModel(structure) +``` + ```@setup logging runtime = round(time() - start_time; digits=2) @info "...done after $runtime s." diff --git a/docs/src/NQCModels/overview.md b/docs/src/NQCModels/overview.md index 8e4567cc8..ecc6c344f 100644 --- a/docs/src/NQCModels/overview.md +++ b/docs/src/NQCModels/overview.md @@ -2,13 +2,20 @@ @info "Expanding src/NQCModels/overview.md..." start_time = time() ``` + # NQCModels.jl +!!! details "Overview of all model types currently implemented." + + ```@raw html + + ``` + To perform nonadiabatic molecular dynamics simulations, it is necessary to define the system Hamiltonian. For simple models, this often comes in the form of small matrix in the diabatic representation but equally the electronic Hamiltonian could be obtained directly -from *ab initio* electronic structure theory. +from _ab initio_ electronic structure theory. `NQCModels.jl` is a package that aims to provide a common interface for defining these models that is flexible enough to allow for a wide range @@ -30,9 +37,11 @@ and [`DiabaticModel`](@ref NQCModels.DiabaticModels.DiabaticModel). The [`AdiabaticModel`](@ref NQCModels.AdiabaticModels.AdiabaticModel) is used for adiabatic dynamics, providing only the potential and force used in classical mechanics. -The [`DiabaticModel`](@ref NQCModels.DiabaticModels.DiabaticModel) is used for nonadiabatic dynamics, +The [`DiabaticModel`](@ref NQCModels.DiabaticModels.DiabaticModel) is used for nonadiabatic dynamics, where the potential is instead a `Hermitian` matrix. +## Using `Model`s + In the [Getting started](@ref) section we briefly touched on how the [`AdiabaticModel`](@ref NQCModels.AdiabaticModels.AdiabaticModel) works and introduced one of the included models. @@ -87,25 +96,18 @@ with respect to each degree of freedom. In this case, the `Matrix` has `size = (1, 2)`, but it should be clear how this can extend to arbitrary numbers of atoms and degrees of freedom for complex models. -The models currently available can be seen in type tree of the -[`Model`](@ref NQCModels.Model) below. -The leaves of the tree are the concrete models, whereas each branch is one of the abstract -types. -Each of these models can be seen in the [Analytic model library](@ref) and -many shall return later when we investigate the dynamics methods. +## Included models -```@example -import AbstractTrees # hide -import InteractiveUtils: subtypes # hide -import NQCModels: Model # hide -AbstractTrees.children(x::Type) = subtypes(x) # hide -AbstractTrees.print_tree(Model) # hide -``` +`NQCModels.jl` includes both analytical models as well as interfaces to connect to full-dimensional potential energy surfaces. +An overview of the implemented analytical models can be seen in the [Analytic model library](@ref) and +many shall return later when we investigate the dynamics methods. +An overview of interfaces to full-dimensional potential energy surfaces can be seen in the [Full dimensional model library](@ref). !!! note "Contributing new models" To learn more about NQCModels.jl and learn how to implement new models, visit the [developer documentation](@ref devdocs-model). + ```@setup logging runtime = round(time() - start_time; digits=2) @info "...done after $runtime s." diff --git a/docs/src/api/NQCDynamics/analysis.md b/docs/src/api/NQCDynamics/analysis.md index a0e057699..bedf21907 100644 --- a/docs/src/api/NQCDynamics/analysis.md +++ b/docs/src/api/NQCDynamics/analysis.md @@ -5,10 +5,10 @@ start_time = time() # Analysis ```@autodocs -Modules=[NQCDynamics.Analysis, NQCDynamics.Analysis.Diatomic] +Modules=[NQCDynamics.Analysis, NQCDynamics.Analysis.Diatomic, NQCDynamics.Analysis.RigidRotator, NQCDynamics.Analysis.Postprocess] ``` ```@setup logging runtime = round(time() - start_time; digits=2) @info "...done after $runtime s." -``` \ No newline at end of file +``` diff --git a/docs/src/api/NQCDynamics/structure.md b/docs/src/api/NQCDynamics/structure.md index 585f3703e..a2022d147 100644 --- a/docs/src/api/NQCDynamics/structure.md +++ b/docs/src/api/NQCDynamics/structure.md @@ -2,20 +2,21 @@ @info "Expanding src/api/NQCDynamics/structure.md..." start_time = time() ``` + # Structure -This submodule contains utility functions to analyse and modify atomic structure data, such as interatomic distances, centres of mass, both with and without support for periodic boundary conditions. -These functions can be used to build more sophisticated output functions, or for basic analysis of simulation results in post. +This submodule contains utility functions to analyse and modify atomic structure data, such as interatomic distances, centres of mass, both with and without support for periodic boundary conditions. + +These functions can be used to build more sophisticated output functions, or for basic analysis of simulation results in post. This module **doesn't contain**: + - Basic definitions of atomic structures (e.g. `Atoms`, `PeriodicCell`, ...). These are defined in [`NQCBase`](@ref). - Functions to generate atomic structures. These should be added to [`NQCDynamics.InitialConditions`](@ref). -- Analysis functions for specific systems (e.g. molecules on surfaces). These should be added to [`NQCDynamics.Analysis`](@ref). - +- Analysis functions for specific systems (e.g. molecules on surfaces). These should be added to [`NQCDynamics.Analysis`](@ref). ## Method reference - ```@autodocs Modules=[NQCDynamics.Structure] ``` diff --git a/docs/src/api/NQCModels/mace.md b/docs/src/api/NQCModels/mace.md new file mode 100644 index 000000000..c81a5ed65 --- /dev/null +++ b/docs/src/api/NQCModels/mace.md @@ -0,0 +1,16 @@ +```@setup logging +@info "Expanding src/api/NQCModels/mace.md..." +start_time = time() +``` + +# MACEModels.jl + +```@autodocs +Modules = [MACEModels] +Private = true +``` + +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/assets/compositemodels/struct-explainer.svg b/docs/src/assets/compositemodels/struct-explainer.svg new file mode 100644 index 000000000..3f1a32774 --- /dev/null +++ b/docs/src/assets/compositemodels/struct-explainer.svg @@ -0,0 +1,8 @@ + + + + + + + + CompositeModelSubsystemmodel: NQCModels.modelindices: Vector{Int}Subsystemmodel: NQCModels.ElectronicFrictionProviderindices: Vector{Int}Subsystemmodel: NQCModels.ElectronicFrictionProviderindices: Vector{Int}SubsystemNQCModels.Model(which model to use)indices::Vector{Int}(which atoms the modelshould affect)Thermostattemperature::FunctionNumber(what temperature to set)indices::Vector{Int}(which atoms to apply thetemperature to) \ No newline at end of file diff --git a/docs/src/assets/figures/Model_types.svg b/docs/src/assets/figures/Model_types.svg new file mode 100644 index 000000000..50dae2142 --- /dev/null +++ b/docs/src/assets/figures/Model_types.svg @@ -0,0 +1,2203 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/mdef-ttm-example/2TM_graphic.mp4 b/docs/src/assets/mdef-ttm-example/2TM_graphic.mp4 new file mode 100644 index 000000000..4d555e0a0 Binary files /dev/null and b/docs/src/assets/mdef-ttm-example/2TM_graphic.mp4 differ diff --git a/docs/src/assets/mdef-ttm-example/Thermostat_Strategies.svg b/docs/src/assets/mdef-ttm-example/Thermostat_Strategies.svg new file mode 100644 index 000000000..5ec9a51c8 --- /dev/null +++ b/docs/src/assets/mdef-ttm-example/Thermostat_Strategies.svg @@ -0,0 +1,586 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/dynamicssimulations/dynamicsmethods/mdef.md b/docs/src/dynamicssimulations/dynamicsmethods/mdef.md index 9b4719e56..d4f19ddb8 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/mdef.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/mdef.md @@ -8,7 +8,6 @@ start_time = time() A set of fundamental and technologically relevant chemical processes (surface scattering, dissociative chemisorption, surface diffusion, recombinative desorption, etc.) are often catalyzed at the metal surface of several late transition metals (Au, Ag, Cu, Pt, Pd, Rh, etc). These metallic surfaces, unlike other surfaces, are characterized by a dense manifold of electronic states at the Fermi level, which produce continuous conduction and valence bands without a band gap. A theoretical description of the chemical processes at these metal surfaces is often challenging due to the Born-Oppenheimer (BO) approximation no longer being valid. With the breakdown of the Born-Oppenheimer approximation, nonadiabatic effects have to be considered to describe, e.g., the energy exchange that can take place between adsorbate and substrate degrees of freedom (DOF). - A fully quantum dynamical approach of this complex scenario is currently unfeasible and the gas-surface reaction dynamics are often described using quasi-classical methods where nuclear motion is described classically. Molecular dynamics with electronic friction (MDEF) is one of main methods used to deal with the nonadiabaticity in gas-surface chemical reactions. MDEF has been widely employed to decribe and simulate the nuclear dynamics in several molecular systems. It is a theoretical model based on a ground-state Langevin equation of motion which introduces nonadiabatic effects by using frictional and stochastic forces. This approach was originally introduced by Head-Gordon and Tully and the nonadiabatic effects can be included through different electronic friction models (see section below, LDFA and TDPT). The nuclear coordinates of the adsorbate atoms evolve as follows: @@ -101,7 +100,7 @@ Local density friction approximation (LDFA) is a theoretical model which describ of each atom. In our current LDFA implementation, a set of pre-calculated electronic friction coefficients (``\eta_{e,i}``) computed at different Wigner-Seitz radius (``r_s``) are used to fit and get an analytical expression to connect any ``r_s`` values with an single electronic friction coefficient by means of -cubic Spline functions. The Wigner-Sietz radius is connected to the metal substrate electron density by the following equation, +cubic Spline functions. The Wigner-Seitz radius is connected to the metal substrate electron density by the following equation, ```math r_s(\rho) = (\frac{3}{4\pi \rho (\mathbf{r_{i}})})^{1/3} @@ -125,3 +124,25 @@ View the [friction models page](@ref models-friction) to learn about how this ca runtime = round(time() - start_time; digits=2) @info "...done after $runtime s." ``` + +## Use of Composite models for phonon thermostatting + +```julia +# PES applies to all atoms +pes_subsystem = Subsystem(pes_model) + +# Electronic friction applies to adsorbate atoms 55,56 +electronic_friction = Subsystem(adsorbate_friction_model, indices=[55,56]) + +# Combine models and generate Simulation +combined_model = CompositeModel(pes_subsystem, electronic_friction) + +sim_T_el_only = Simulation{MDEF}( + atoms, + combined_model; + temperature = T_el_function, + cell=cell +) +``` + +See the [example page](@ref mdef-ttm-example) for a longer explanation on how to compose multiple models. diff --git a/docs/src/dynamicssimulations/dynamicssimulations.md b/docs/src/dynamicssimulations/dynamicssimulations.md index be5d4d953..9a406fd4c 100644 --- a/docs/src/dynamicssimulations/dynamicssimulations.md +++ b/docs/src/dynamicssimulations/dynamicssimulations.md @@ -32,6 +32,28 @@ sim = Simulation{Ehrenfest}(atoms, TullyModelOne(); temperature=0, cell=Infinite Here we have initialised the simulation parameters, including the default temperature and cell explicitly. `sim` takes the place of the `p` parameter seen throughout [DifferentialEquations](https://diffeq.sciml.ai/stable/). +### Simulation Temperature + +The temperature of a simulation isn't used by some dynamics methods, but can always be included. +Methods such as [Langevin dynamics](@ref langevin-dynamics) or [MDEF](@ref mdef-dynamics) use temperature to determine the magnitude of random flucturations to preserve detailed balance. + +The temperature for a simulation can be provided in three ways: + +**Fixed temperature** + +This can be a number, either in atomic units or a Unitful quantity (e.g. `300u"K"`). + +**Temperature function** + +For time-dependent temperatures, the temperature can be set to a function. This function must take the time in Unitful `u"fs"` as an input and return a temperature in Unitful `u"K"`. + +**Different temperatures per atom** + +If different parts of the simulated system should be affected by different effective temperatures, a `Vector` of [`Temperature`](@ref)s can be passed as the temperature argument. +However, only one temperature should be applied to each atom in the system. + +### DynamicsVariables + For [DifferentialEquations](https://diffeq.sciml.ai/stable/) to allow for a wide variety of solvers, the input arrays [`DynamicsVariables`](@ref) must be [`AbstractArray`](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)s. In nonadiabatic dynamics simulations, we usually have different groups of variables that behave in particular ways. @@ -57,6 +79,8 @@ This helps to ensure each method follows a similar workflow, making it easy to s The output of this function takes the place of the `u` argument seen throughout [DifferentialEquations](https://diffeq.sciml.ai/stable/). +### Running dynamics + With both the [`Simulation`](@ref) and [`DynamicsVariables`](@ref) in hand, the central function is [`run_dynamics`](@ref) which allows us to perform a single dynamics trajectory. [`run_dynamics`](@ref) takes the simulation parameters `sim` and the initial conditions `u0`, along with a time span `tspan` diff --git a/docs/src/examples/mdef_multithermostat.md b/docs/src/examples/mdef_multithermostat.md new file mode 100644 index 000000000..78f918f67 --- /dev/null +++ b/docs/src/examples/mdef_multithermostat.md @@ -0,0 +1,199 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/mdef.md..." +start_time = time() +``` +# [MDEF with the Two-temperature model and separate Electron and Phonon Thermostats](@id mdef-ttm-example) + +!!! note + + This page focuses on using the model combination tools in NQCDynamics.jl and NQCModels.jl to perform MDEF simulations with separate thermostats for different parts of the system. + The page on [Molecular Dynamics with Electronic Friction](@ref mdef-dynamics) contains more information on MDEF as a method, and is a better starting point to setting up simulations. + +## Intro🦆tion +On this page, we will look at the example of simulating the response of hydrogen adsorbed to a Cu(111) surface to an ultrafast laser pulse using [Molecular Dynamics with Electronic Friction](@ref mdef-dynamics) and the Two-temperature model. + +As detailed before, the nuclear coordinates in MDEF evolve as follows: +```math +\mathbf{M}\ddot{\mathbf{R}} = - \nabla_R V(\mathbf{R}) + \mathbf{F}(t) - \Gamma(\mathbf{R}) \dot{\mathbf{R}} +``` +The first term on the right hand side of the equation corresponds to a conservative force associated with the potential energy surface (PES) as in the adiabatic case. The third term is the friction force and it comes from multiplication between the electronic friction object (``\Gamma(\mathbf{R})``) and the velocity. +Finally, the second term is a temperature and friction-dependent stochastic force which ensures detailed balance. + +## The two-temperature model + +!!! info + + Temperature profiles for the Two-temperature model can be generated using the [LightMatter.jl](https://github.com/maurergroup/LightMatter.jl) package. (Will be released soon) + + +The two-temperature model (TTM) is a system of two coupled differential equations that describe the time evolution of the electronic and lattice temperatures, $T_{el}$ and $T_{ph}$. In this work, lateral heat transport across the metal surface is neglected, allowing for a simplified one-dimensional representation of the two-temperature model. The TTM is given by the following set of equations: + +**Electronic heat capacity** +```math +C_{el}\frac{\partial}{\partial t}T_{el}=\nabla(\kappa\nabla T_{el})-g(T_{el}-T_{ph})+S(z,t) +``` + +**Phononic heat capacity** +```math +C_{ph}\frac{\partial}{\partial t}T_{ph}=g(T_{el}-T_{ph}) +``` + +**Laser source term** +```math +S(z,t)=F\cdot\frac{\exp{\frac{-z}{\xi}}\cdot\exp{\frac{-(t-t_0)^2}{2r^2}}}{\xi\cdot\sqrt{2\pi}\tau}; + \tau=\frac{\mathrm{FWHM}}{2\sqrt{2\ln{2}}} +``` + +The model is parametrised with material-specific parameters obtained from ab-initio calculations or experiments. As an example, the material-specific parameters for bulk crystalline Cu are shown here ([Scholz2019](@cite)): + +| Parameter | Explanation | Value | +|-----------------------|---------------------------------------------------------------|-----------------------------------------------| +| ``g`` | Electron-phonon coupling constant | ``1\times 10^{-17}~\mathrm{W~m^{-3}~K^{-1}}`` | +| ``\xi`` | Laser penetration depth for ``\lambda=400~\mathrm{nm}`` | ``\frac{1}{\xi}=14.9~\mathrm{nm}`` | +| ``\gamma_{el}`` | Scaling constant for the electronic heat capacity ``C_{el}``. | ``98.0~\mathrm{J~m^{-3}~K^{-2}}`` | +| ``\theta_D`` | Debye temperature of Cu. | ``343~\mathrm{K}`` | +| ``N`` | Atom density of bulk Cu | ``8.5\times 10^{28}`` | +| ``\kappa_{RT}`` | Thermal conductivity at room temperature | ``401.0~\mathrm{W^{-1}~K^{-1}}`` | +| ``\mathrm{FWHM}`` | Full width half maximum of the laser pulse | ``150~\mathrm{fs}`` | + + + + +The electronic heat capacity $C_{el}(T_{el})$ was determined using a linear scaling relation: + +```math +C_{el}(T_{el})=\gamma_{el}T_{el} +``` + +The temperature dependence of the thermal conductivity of Cu was modelled using the following relation: + +```math +\kappa(T_{el}, T_{ph})=\kappa_{RT}\frac{T_{el}}{T_{ph}} +``` + +The phononic heat capacity was calculated using the Debye model with $N$ and $\theta_D$ values given above: + +```math +C_{ph}(T_{ph})=9nk_B\left(\frac{T_{ph}}{\theta_D}\right)^3\int_0^{\frac{\theta_D}{T_{ph}}}\frac{x^4e^x}{(e^x-1)^2}~dx +``` + + +An example of the TTM progression for a $80~\mathrm{J~m^{-2}}$ absorbed fluence is shown in the video below. + +![Temperature progressions, energy distributions and laser power for a TTM progression.](../assets/mdef-ttm-example/2TM_graphic.mp4) + +In our implementation of the TTM, the temperature progressions are saved to CSV files, which need to be interpolated to provide a continuous temperature function to use for our `Simulation`s. + +```julia +using BSplineKit +using Unitful, UnitfulAtomic +using CSV, DataFrames + +""" +This function reads a CSV file for a TTM solution containing: +1. The time in ps +2. Electronic temperature +3. Phononic temperature + +Depending on which index is specified, a continuous temperature function for the electronic or phononic temperature is returned. +""" +function T_function_from_file(file::String, index::Int=2) + TTM_file = CSV.read(file, DataFrame) + T_spline = interpolate(TTM_file.Time, TTM_file[:, index], BSplineOrder(4)) # is a cubic spline + T_extrapolation = extrapolate(T_spline, Smooth()) #! Don't use to go earlier than the first point! + T_function(time_ps) = T_extrapolation(ustrip(u"ps", time_ps)) * u"K" + return T_function +end + +T_el_function = T_function_from_file("ttm-file.csv", 2) +T_ph_function = T_function_from_file("ttm-file.csv", 3) +``` + +## Combining MDEF and the TTM +The electron and phonon temperatures are propagated as described above, yielding temperature profiles that are applied to MDEF simulations through the ``\mathcal{W}`` term in the Langevin equation. +``T_{el}`` is applied directly to the adsorbate H atoms, while ``\gamma`` as a function of H positions is determined using one of the MLIPs detailed above. +``T_{ph}`` is coupled to the surface Cu atoms using a position-independent coupling constant. + +![Different thermostatting strategies for our system](../assets/mdef-ttm-example/Thermostat_Strategies.svg) + +As shown above, ``T_{ph}`` can enter our system in different ways. We will use [`Subsystem`](@ref)s to generate a [`CompositeModel`](@ref) and [`Simulation`](@ref) for these different strategies. + +```julia +using NQCDynamics +using Unitful, UnitfulAtomic +using MACEModels + +# Import starting structure +atoms, position, cell = read_extxyz("starting_structure.xyz") + +# Load PES model +pes_model = MACEModel( + atoms, + cell, + ["mace-model.model"]; + device="cpu", + default_dtype=Float32, + mobile_atoms::Vector{Int}=18:56, # Keep the top 4 layers mobile, freeze the bottom 2. +) + +# Load adsorbate and surface friction models +adsorbate_friction_model = RandomFriction(3) +surface_friction_model = ConstantFriction(γ=0.2) +``` + +Now that all models are defined, we can combine them in different ways for a simulation: + +### Electron thermostat only + +```julia +# PES applies to all atoms +pes_subsystem = Subsystem(pes_model) + +# Electronic friction applies to adsorbate atoms 55,56 +electronic_friction = Subsystem(adsorbate_friction_model, indices=[55,56]) + +# Combine models and generate Simulation +combined_model = CompositeModel(pes_subsystem, electronic_friction) + +sim_T_el_only = Simulation{MDEF}( + atoms, + combined_model; + temperature = T_el_function, + cell=cell +) +``` + +### Phonon thermostatting top four layers +```julia +# PES applies to all atoms +pes_subsystem = Subsystem(pes_model) + +# Electronic friction applies to adsorbate atoms 55,56 +electronic_friction = Subsystem(adsorbate_friction_model, indices=[55,56]) + +# Phononic friction applies to surface atoms 18:54 +phononic_friction = Subsystem(surface_friction_model, indices=18:54) + +# Combine models and generate Simulation +combined_model = CompositeModel( + pes_subsystem, + electronic_friction, + phononic_friction +) + +# Create Temperature objects to apply T_el to adsorbates, T_ph to surface +# Manually specified indices +electron_thermostat = Temperature(T_el_function, indices = [55,56]) +# Inherit indices from a Subsystem +phonon_thermostat = Temperature(T_ph_function, phononic_friction) + +sim_T_el_only = Simulation{MDEF}( + atoms, + combined_model; + temperature = [ + electron_thermostat, + phonon_thermostat + ], + cell=cell +) +``` diff --git a/docs/src/examples/reactive_scattering.md b/docs/src/examples/reactive_scattering.md index d91f48a27..78c26d8da 100644 --- a/docs/src/examples/reactive_scattering.md +++ b/docs/src/examples/reactive_scattering.md @@ -32,7 +32,7 @@ surface (`height` keyword) with a normal incidence. As usual, all quantities default to atomic units. Here we use [Unitful](https://painterqubits.github.io/Unitful.jl/stable/) to input the translational energy and height using different units, where they are later converted internally. -```@example h2scatter +```julia @example h2scatter using NQCDynamics using NNInterfaces using Unitful @@ -79,7 +79,7 @@ terminate. This means we can set our time span relatively long since we expect most simulations to terminate before reaching the time limit. -```@example h2scatter +```julia @example h2scatter using Statistics: mean using LinearAlgebra: norm @@ -117,7 +117,7 @@ This takes a `.cube` file containing the electron density and will provide the f during the dynamics. Here we have given the new model our model from above, which will provide the forces. -```@example h2scatter +```julia @example h2scatter using CubeLDFAModel model = LDFAModel(model, "../assets/friction/test.cube", atoms, friction_atoms=[1,2], cell) ``` diff --git a/ext/MACEModelsext.jl b/ext/MACEModelsext.jl new file mode 100644 index 000000000..91c16b9d3 --- /dev/null +++ b/ext/MACEModelsext.jl @@ -0,0 +1,41 @@ +module MACEModelsext +""" +To take advantage of batch evaluation in ring-polymer MD simulations, we need to overload potential +and derivative evaluation in the ring polymer calculators in NQCDynamics. +""" + +using MACEModels +using NQCDynamics: Calculators + +function Calculators.evaluate_potential!(calc::RingPolymerAdiabaticCalculator{T,MACEModel}, R::AbstractArray{S,3}) where {T,S} + calc.stats[:potential] += 1 + potential_vector = @views [R[:, :, i] for i in axes(R, 3)] + calc.potential = NQCModels.potential(calc.model, calc.model.atoms, potential_vector, calc.model.cell) + return nothing +end + +function Calculators.evaluate_derivative!(calc::RingPolymerAdiabaticCalculator{T,MACEModel}, R::AbstractArray{S,3}) where {T,S} + calc.stats[:derivative] += 1 + derivative_vector = @views [R[:, :, i] for i in axes(R, 3)] + NQCModels.derivative!(calc.model, calc.model.atoms, calc.derivative, derivative_vector, calc.model.cell) + return nothing +end + +function Calculators.evaluate_potential!(calc::RingPolymerFrictionCalculator{T,MACEModel}, R::AbstractArray{S,3}) where {T,S} + @debug "RPMDEF accelerated potential evaluation" + calc.stats[:potential] += 1 + potential_vector = @views [R[:, :, i] for i in axes(R, 3)] + calc.potential = NQCModels.potential(calc.model, calc.model.atoms, potential_vector, calc.model.cell) + return nothing +end + +function Calculators.evaluate_derivative!(calc::RingPolymerFrictionCalculator{T,MACEModel}, R::AbstractArray{S,3}) where {T,S} + @debug "RPMDEF accelerated derivative evaluation" + calc.stats[:derivative] += 1 + derivative_vector = @views [R[:, :, i] for i in axes(R, 3)] + NQCModels.derivative!(calc.model, calc.model.atoms, calc.derivative, derivative_vector, calc.model.cell) + return nothing +end + + +end diff --git a/src/Analysis/Analysis.jl b/src/Analysis/Analysis.jl index 888e044ac..9d374ac69 100644 --- a/src/Analysis/Analysis.jl +++ b/src/Analysis/Analysis.jl @@ -4,15 +4,16 @@ Analysis functions common enough to be included in the main package. module Analysis using Reexport: @reexport -# Diatomic analysis functions @alexsp32 +# Diatomic analysis functions include("diatomic.jl") export Diatomic -# Rebinding of quantise_diatomic under Analysis. -using NQCDynamics: InitialConditions -function quantise_diatomic(sim, v, r; args...) - InitialConditions.QuantisedDiatomic.quantise_diatomic(sim, v, r; args...) -end -export quantise_diatomic +# Allow re-processing of simulation data. +include("postprocess.jl") +export Postprocess + +# Rigid rotator energies +include("rigid_rotator.jl") +export RigidRotator -end \ No newline at end of file +end diff --git a/src/Analysis/diatomic.jl b/src/Analysis/diatomic.jl index 3bba9e8fc..e8a74601a 100644 --- a/src/Analysis/diatomic.jl +++ b/src/Analysis/diatomic.jl @@ -1,5 +1,5 @@ """ -Analysis functions for surface chemistry of diatomic molecules. +Analysis functions for surface chemistry of diatomic molecules. """ module Diatomic using NQCDynamics: AbstractSimulation, Simulation, get_positions, Structure @@ -9,24 +9,37 @@ using LinearAlgebra using Statistics # Surface normal projection -surface_normal_height(x::AbstractVector, surface_normal::AbstractVector=[0,0,1])=norm(x.*normalize(surface_normal)) +surface_normal_height(x::AbstractVector, surface_normal::AbstractVector=[0, 0, 1]) = norm(x .* normalize(surface_normal)) """ surface_distance_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_distance_threshold=5.0*u"Å") Checks that the diatomic molecule is at least `surface_distance_threshold` away from the highest substrate atom in the simulation. """ -function surface_distance_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_distance_threshold=5.0*u"Å") - molecule_position=surface_normal_height(Structure.pbc_center_of_mass(x, indices..., simulation)) # molecule position wrt surface normal - substrate_heights=surface_normal_height.(eachcol(get_positions(x)[:, symdiff(1:length(simulation.atoms.masses), indices)])) # substrate positions along surface normal - highest_z=max(substrate_heights[substrate_heights.≤molecule_position]...) # Ignore substrate above molecule. - #? surface_distance_threshold must be < vacuum above surface and > distance between substrate layers. - if au_to_ang(molecule_position-highest_z)*u"Å"≥surface_distance_threshold - @debug "Surface distance condition evaluated true." - return true - else - return false - end +function surface_distance_condition( + x::AbstractArray, # DynamicsVariables + indices::Vector{Int}, # Which indices are the diatomic species? All others are counted as substrate. + simulation::AbstractSimulation; + surface_distance_threshold=austrip(5.0 * u"Å"), # Condition is true when the centre of mass of the diatomic is farther than this distance from the highest other atom. + surface_normal=Float64[0, 0, 1] # Unit direction orthogonal to the surface +) + molecule_position = surface_normal_height( + Structure.pbc_center_of_mass(x, indices..., simulation), + surface_normal + ) # Height of the molecule in the surface normal direction + + # Get the height of all substrate atoms in surface normal direction. + substrate_heights = [surface_normal_height(get_positions(x)[:, substrate_atom_id]) for substrate_atom_id in symdiff(1:length(simulation.atoms.masses), indices, surface_normal)] + + # Ignore substrate above molecule in case PBC wrapping puts one above the diatomic + highest_z = max(substrate_heights[substrate_heights.≤molecule_position]...) + + if molecule_position - highest_z ≥ surface_distance_threshold + @debug "Surface distance condition evaluated true." + return true + else + return false + end end """ @@ -34,62 +47,86 @@ end Evaluates true if the centre of mass velocity vector of the diatomic molecule points to the surface. """ -function com_velocity_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_normal::AbstractVector=[0,0,1]) - if dot(Structure.velocity_center_of_mass(x, indices..., simulation), normalize(surface_normal))>0 - return false - else - return true - end +function com_velocity_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; surface_normal::AbstractVector=[0, 0, 1]) + if dot(Structure.velocity_center_of_mass(x, indices..., simulation), normalize(surface_normal)) > 0 + return false + else + return true + end end """ - get_desorption_frame(trajectory::Vector, diatomic_indices::Vector{Int}, simulation::AbstractSimulation;surface_normal::Vector=[0,0,1], surface_distance_threshold=5.0*u"Å") + close_approach_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; threshold = 1.5u"Å") -Determines the index in a trajectory where surface desorption begins. + Evaluate true if the diatomic bond length is below `threshold`. +""" +function close_approach_condition(x::AbstractArray, indices::Vector{Int}, simulation::AbstractSimulation; threshold = 1.5u"Å") + if Structure.pbc_distance(x, indices..., simulation) ≤ threshold + return true + else + return false + end +end + +""" + get_desorption_frame(trajectory::AbstractVector, diatomic_indices::Vector{Int}, simulation::AbstractSimulation; surface_normal::Vector=[0, 0, 1], surface_distance_threshold=5.0 * u"Å", fallback_distance_threshold = 1.5u"Å") + +Determines the index in a trajectory where surface desorption begins. This is evaluated using two conditions: -1. In the trajectory, the diatomic must be `surface_distance_threshold` or further away from the highest other atom. (In `surface_normal` direction). +1. In the trajectory, the diatomic must be `surface_distance_threshold` or further away from the highest other atom. (In `surface_normal` direction). + +2. Desorption begins at the turning point of the centre of mass velocity component along `surface_normal`, indicating overall movement away from the surface. -2. Desorption begins at the turning point of the centre of mass velocity component along `surface_normal`, indicating overall movement away from the surface. +If the second condition is never reached (can happen for particularly quick desorptions), the `fallback_distance_threshold` is used to find the last point where the +diatomic bond length is above the given value and saves from that point onwards. """ -function get_desorption_frame(trajectory::AbstractVector, diatomic_indices::Vector{Int}, simulation::AbstractSimulation;surface_normal::Vector = [0,0,1], surface_distance_threshold = 5.0*u"Å") - desorbed_frame=findfirst(surface_distance_condition.(trajectory, Ref(diatomic_indices), Ref(simulation); surface_distance_threshold=surface_distance_threshold)) - +function get_desorption_frame(trajectory::AbstractVector, diatomic_indices::Vector{Int}, simulation::AbstractSimulation; surface_normal::Vector=[0, 0, 1], surface_distance_threshold=5.0 * u"Å", fallback_distance_threshold = 1.5u"Å") + desorbed_frame = findfirst(surface_distance_condition.(trajectory, Ref(diatomic_indices), Ref(simulation); surface_distance_threshold=surface_distance_threshold)) + if isa(desorbed_frame, Nothing) @debug "No desorption event found." return nothing else @debug "Desorption observed in snapshot $(desorbed_frame)" - leaving_surface_frame=findlast(com_velocity_condition.(view(trajectory, 1:desorbed_frame), Ref(diatomic_indices), Ref(simulation); surface_normal=surface_normal)) #ToDo testing views for memory efficiency, need to test time penalty. Also need to test if running on everything and findfirst-ing the Bool array is quicker. - if leaving_surface_frame<0 - @warn "Desorption event detected, but no direction change detected." - return nothing + leaving_surface_frame = findlast(com_velocity_condition.(view(trajectory, 1:desorbed_frame), Ref(diatomic_indices), Ref(simulation); surface_normal=surface_normal)) #ToDo testing views for memory efficiency, need to test time penalty. Also need to test if running on everything and findfirst-ing the Bool array is quicker. + if isnothing(leaving_surface_frame) + @debug "Centre of mass velocity criterion was never met. Falling back to distance threshold." + leaving_surface_frame = findlast(close_approach_condition.(view(trajectory, 1:desorbed_frame), Ref(diatomic_indices), Ref(simulation); threshold = fallback_distance_threshold)) + if isnothing(leaving_surface_frame) + @warn begin + println("H-H distance threshold was never met. Something is wrong in desorption detection logic - Returning entire trajectory for debugging.") + end + return 1 + else + return leaving_surface_frame + end else return leaving_surface_frame end end end -function get_desorption_angle(trajectory::AbstractVector, indices::Vector{Int}, simulation::AbstractSimulation; surface_normal = [0,0,1], surface_distance_threshold = 5.0*u"Å") - # First determine where the reaction occurred on the surface. - desorption_frame = get_desorption_frame(trajectory, indices, simulation;surface_distance_threshold = surface_distance_threshold, surface_normal = surface_normal) +function get_desorption_angle(trajectory::AbstractVector, indices::Vector{Int}, simulation::AbstractSimulation; surface_normal=[0, 0, 1], surface_distance_threshold=5.0 * u"Å") + # First determine where the reaction occurred on the surface. + desorption_frame = get_desorption_frame(trajectory, indices, simulation; surface_distance_threshold=surface_distance_threshold, surface_normal=surface_normal) if isa(desorption_frame, Nothing) @debug "No desorption event detected in trajectory" return nothing end @debug "Desorption frame: $(desorption_frame)" - # Determine the average centre of mass velocity to decrease error due to vibration and rotation orthogonal to true translational component. - com_velocities = zeros(eltype(trajectory[1]), length(surface_normal), length(trajectory)-desorption_frame) + # Determine the average centre of mass velocity to decrease error due to vibration and rotation orthogonal to true translational component. + com_velocities = zeros(eltype(trajectory[1]), length(surface_normal), length(trajectory) - desorption_frame) for i in 1:length(trajectory)-desorption_frame com_velocities[:, i] .= Structure.velocity_center_of_mass(trajectory[i+desorption_frame], indices[1], indices[2], simulation) end - average_velocity=mean(com_velocities;dims=2) + average_velocity = mean(com_velocities; dims=2) # Now convert into an angle by arccos((a•b)/(|a|*|b|)) return Structure.angle_between(vec(average_velocity), surface_normal) end -# 6×6 Cartesian to internal coordinate transformation. +# 6×6 Cartesian to internal coordinate transformation. function transform_to_internal_coordinates(to_transform::Matrix, config::Matrix, index1::Int, index2::Int, sim::Simulation) U = transform_U(config, index1, index2, sim) # Generate transformation matrix return transpose(U) * to_transform * U @@ -101,7 +138,7 @@ function transform_from_internal_coordinates(to_transform::Matrix, config::Matri end function transform_r(config, index1::Int, index2::Int) - return sqrt(sum(mapslices(x->(x[2]-x[1])^2, config[:, [index1, index2]];dims=2))) + return sqrt(sum(mapslices(x -> (x[2] - x[1])^2, config[:, [index1, index2]]; dims=2))) end function transform_r1(config, index1::Int, index2::Int) @@ -116,31 +153,31 @@ Builds diatomic Cartesian to internal coordinate transformation matrix as descri """ function transform_U(config::Matrix, index1::Int, index2::Int, sim::Simulation) masses = Structure.fractional_mass(sim, index1, index2) - config[:,index2] .+= Structure.minimum_distance_translation(config, index1, index2, sim) # PBC wrap positions for correct transformation. + config[:, index2] .+= Structure.minimum_distance_translation(config, index1, index2, sim) # PBC wrap positions for correct transformation. r = transform_r(config, index1, index2) r1 = transform_r1(config, index1, index2) unity = LinearAlgebra.I(3) U_i1 = vcat( - mapslices(x->(x[1]-x[2])/r, config[:, [index1, index2]]; dims=2), - mapslices(x->(x[2]-x[1])/r, config[:, [index2, index1]]; dims=2), + mapslices(x -> (x[1] - x[2]) / r, config[:, [index1, index2]]; dims=2), + mapslices(x -> (x[2] - x[1]) / r, config[:, [index2, index1]]; dims=2), ) U_i2 = vcat( - mapslices(x->(x[2]-x[1])*(config[3, index2]-config[3, index1])/(r^2*r1), config[1:2, [index1, index2]]; dims=2), - -r1/r^2, - mapslices(x->(x[1]-x[2])*(config[3, index2]-config[3, index1])/(r^2*r1), config[1:2, [index2, index1]]; dims=2), - r1/r^2, + mapslices(x -> (x[2] - x[1]) * (config[3, index2] - config[3, index1]) / (r^2 * r1), config[1:2, [index1, index2]]; dims=2), + -r1 / r^2, + mapslices(x -> (x[1] - x[2]) * (config[3, index2] - config[3, index1]) / (r^2 * r1), config[1:2, [index2, index1]]; dims=2), + r1 / r^2, ) U_i3 = [ - -(config[2, index1]-config[2, index2]), - (config[1, index2]-config[1, index1]), + -(config[2, index1] - config[2, index2]), + (config[1, index2] - config[1, index1]), 0.0, - -(config[2, index2]-config[2, index1]), - (config[1, index1]-config[1, index2]), + -(config[2, index2] - config[2, index1]), + (config[1, index1] - config[1, index2]), 0.0, ] ./ (r1^2) - return hcat(U_i1, U_i2, U_i3, vcat(unity.*masses[1], unity.*masses[2])) + return hcat(U_i1, U_i2, U_i3, vcat(unity .* masses[1], unity .* masses[2])) end export get_desorption_frame, get_desorption_angle, transform_from_internal_coordinates, transform_to_internal_coordinates, transform_U -end \ No newline at end of file +end diff --git a/src/Analysis/postprocess.jl b/src/Analysis/postprocess.jl new file mode 100644 index 000000000..a44e445af --- /dev/null +++ b/src/Analysis/postprocess.jl @@ -0,0 +1,46 @@ +module Postprocess +using NQCDynamics: AbstractSimulation, Ensembles + +export apply_output_functions + +struct FakeProblem + p::AbstractSimulation +end + +struct FakeSolution + t::AbstractArray + u::AbstractArray + prob::FakeProblem +end + +""" + apply_output_functions(sol::FakeSolution, output_functions; savetime::Bool=true) + +Evaluates output functions on a DifferentialEquations.jl solution object or fake solution object +generated using a defined Simulation and a `DynamicsVariables` type. + +Basically equivalent to running `run_dynamics()` with the same output functions, but without +doing the dynamics simulation again. +""" +function apply_output_functions(sol, output_functions; savetime::Bool=true) + if !isa(output_functions, Tuple) + output_functions = (output_functions,) + end + output = Ensembles.EnsembleSaver(output_functions, savetime) + return output(sol, 1) +end + +""" + apply_output_functions(sim<:AbstractSimulation, u_type::AbstractVector, t_type::AbstractVector, output_functions; savetime::Bool=true) + +Evaluates output functions on a defined Simulation, a `DynamicsVariables` type output and a time-type output. + +Basically equivalent to running `run_dynamics()` with the same output functions, but without +doing the dynamics simulation again. +""" +function apply_output_functions(sim::AbstractSimulation, u_type::AbstractVector, t_type::AbstractVector, output_functions; savetime::Bool=true) + sol = FakeSolution(t_type, u_type, FakeProblem(sim)) + return apply_output_functions(sol, output_functions; savetime=savetime)[1] +end + +end diff --git a/src/Analysis/rigid_rotator.jl b/src/Analysis/rigid_rotator.jl new file mode 100644 index 000000000..62a32f2f7 --- /dev/null +++ b/src/Analysis/rigid_rotator.jl @@ -0,0 +1,36 @@ +module RigidRotator +using Unitful, UnitfulAtomic +using NQCDynamics: Structure, AbstractSimulation +using LinearAlgebra: norm + +""" + classical_translational_energy(config::Any, ind1::Union{Int, CartesianIndex}, ind2::Union{Int, CartesianIndex}, sim::Simulation) + +Returns the classical translational energy in Hartree +""" +function classical_translational_energy(config::Any, ind1::Union{Int,CartesianIndex}, ind2::Union{Int,CartesianIndex}, sim :: AbstractSimulation) + return sum(sim.atoms.masses[[ind1, ind2]]) / 2 * norm(velocity_center_of_mass(config, ind1, ind2, sim.atoms))^2 +end + +""" + classical_rotation_energy(J::Union{Int, CartesianIndex}, config::Any, ind1::Union{Int, CartesianIndex}, ind2::Union{Int, CartesianIndex}, sim::Simulation) + +Classical rotation energy of a rigid diatomic rotor in Hartree +""" +function classical_rotation_energy(J::Int, config::Any, ind1::Union{Int,CartesianIndex}, ind2::Union{Int,CartesianIndex}, sim :: AbstractSimulation) + I = reduced_mass(sim, ind1, ind2) * austrip(pbc_distance(config, ind1, ind2, sim))^2 # Classical moment of inertia in atomic units. + return J * (J + 1) / (2 * I) +end + +""" + harmonic_vibration_energy(ν::Union{Int, CartesianIndex}, k::Float, ind1::Union{Int, CartesianIndex}=1, ind2::Union{Int, CartesianIndex}=2, sim::Simulation) + +Vibrational energy of a harmonic oscillator with the force constant k and vibrational level ν. +""" +function harmonic_vibration_energy(ν::Union{Int}, k::Float64, ind1::Union{Int,CartesianIndex}, ind2::Union{Int,CartesianIndex}, sim :: AbstractSimulation) + return (ν + 1 / 2) * sqrt(k / reduced_mass(sim, ind1, ind2)) +end + +export classical_translational_energy, classical_rotation_energy, harmonic_vibration_energy + +end diff --git a/src/Calculators/Calculators.jl b/src/Calculators/Calculators.jl index e67e396c8..7e9e2beaf 100644 --- a/src/Calculators/Calculators.jl +++ b/src/Calculators/Calculators.jl @@ -19,7 +19,7 @@ using LinearAlgebra: LinearAlgebra, Hermitian, I, Eigen, tr using StaticArrays: SMatrix, SVector using RingPolymerArrays: get_centroid! -using NQCModels: NQCModels, Model, nstates, mobileatoms, dofs +using NQCModels: NQCModels, Model, nstates, mobileatoms, dofs, Subsystem, CompositeModel using NQCModels.AdiabaticModels: AdiabaticModel using NQCModels.DiabaticModels: DiabaticModel, DiabaticFrictionModel using NQCModels.FrictionModels: AdiabaticFrictionModel @@ -35,10 +35,10 @@ Each concrete calculator contains the `Model` and the fields to store the quanti obtained from the model. """ abstract type AbstractCalculator{T,M<:Model} end -abstract type AbstractAdiabaticCalculator{T,M<:AdiabaticModel} <: AbstractCalculator{T,M} end +abstract type AbstractAdiabaticCalculator{T,M<:Union{AdiabaticModel, CompositeModel}} <: AbstractCalculator{T,M} end abstract type AbstractDiabaticCalculator{T,M<:Union{DiabaticFrictionModel,DiabaticModel}} <: AbstractCalculator{T,M} end abstract type AbstractStaticDiabaticCalculator{T,M} <: AbstractDiabaticCalculator{T,M} end -abstract type AbstractFrictionCalculator{T,M<:AdiabaticFrictionModel} <: AbstractCalculator{T,M} end +abstract type AbstractFrictionCalculator{T,M<:Union{AdiabaticFrictionModel, CompositeModel}} <: AbstractCalculator{T,M} end NQCModels.nstates(calc::AbstractCalculator) = NQCModels.nstates(calc.model) NQCModels.eachstate(calc::AbstractCalculator) = NQCModels.eachstate(calc.model) diff --git a/src/Calculators/friction.jl b/src/Calculators/friction.jl index 084c20952..b8aaf597e 100644 --- a/src/Calculators/friction.jl +++ b/src/Calculators/friction.jl @@ -65,6 +65,19 @@ end function Calculator(model::AdiabaticFrictionModel, atoms::Integer, beads::Integer, t::Type{T}) where {T} RingPolymerFrictionCalculator{t}(model, atoms, beads) end +function Calculator(model::CompositeModel, atoms::Integer, t::Type{T}) where {T} + if any([!isa(s.model, AdiabaticModel) for s in NQCModels.get_pes_models(model.subsystems)]) + throw(ArgumentError("Currently, only CompositeModels using AdiabaticModels to supply a PES are supported. ")) + end + FrictionCalculator{t}(model, atoms) +end + +function Calculator(model::CompositeModel, atoms::Integer, beads::Integer, t::Type{T}) where {T} + if any([!isa(s.model, AdiabaticModel) for s in NQCModels.get_pes_models(model.subsystems)]) + throw(ArgumentError("Currently, only CompositeModels using AdiabaticModels to supply a PES are supported. ")) + end + RingPolymerFrictionCalculator{t}(model, atoms, beads) +end function evaluate_friction!(calc::AbstractFrictionCalculator, R::AbstractMatrix) calc.stats[:friction] += 1 diff --git a/src/DynamicsMethods/IntegrationAlgorithms/bcocb.jl b/src/DynamicsMethods/IntegrationAlgorithms/bcocb.jl index 024446535..9bea133b2 100644 --- a/src/DynamicsMethods/IntegrationAlgorithms/bcocb.jl +++ b/src/DynamicsMethods/IntegrationAlgorithms/bcocb.jl @@ -58,8 +58,8 @@ function FrictionCache(sim::RingPolymerSimulation{<:DynamicsMethods.ClassicalMet c1 = Diagonal(zeros(DoFs*atoms, DoFs*atoms)) c2 = zero(c1) - sqrtmass = 1 ./ sqrt.(repeat(sim.atoms.masses; inner=DoFs)) - σ = zero(sqrtmass) + sqrtmass = 1 ./ sqrt.(sim.atoms.masses) + σ = zero(repeat(sqrtmass; inner=DoFs)) MDEFCache(flatvtmp, tmp1, tmp2, gtmp, noise, c1, c2, sqrtmass, σ) end @@ -76,8 +76,8 @@ function FrictionCache(sim::RingPolymerSimulation{<:DynamicsMethods.ClassicalMet γ = [γ0, 2sqrt.(2sim.beads.normal_mode_springs[2:end])...] c1 = exp.(-dt.*γ) c2 = sqrt.(1 .- c1.^2) - sqrtmass = 1 ./ repeat(sqrt.(sim.atoms.masses'), ndofs(sim)) - σ = zero(sqrtmass) + sqrtmass = 1 ./ sqrt.(sim.atoms.masses') + σ = zero(repeat(sqrtmass, ndofs(sim))) LangevinCache(c1, c2, sqrtmass, σ) end diff --git a/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl b/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl index e3a1575ad..6aacb707d 100644 --- a/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl +++ b/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl @@ -130,7 +130,7 @@ end # O Λ = integrator.g(u2,p,t+dt*half) # friction tensor # noise strength: σ = square root of (temperature / mass) for each atom - σ = sqrt(get_temperature(p, t+dt*half)) ./ sqrt.(repeat(p.atoms.masses; inner=ndofs(p))) + σ = repeat(@. sqrt(get_temperature(p, t+dt*half) / p.atoms.masses);inner=ndofs(p)) # eigen decomposition of Λ γ, c = LAPACK.syev!('V', 'U', Λ) # symmetric eigen clamp!(γ, 0, Inf) diff --git a/src/DynamicsMethods/IntegrationAlgorithms/steps.jl b/src/DynamicsMethods/IntegrationAlgorithms/steps.jl index e4cd322f3..ccdd09352 100644 --- a/src/DynamicsMethods/IntegrationAlgorithms/steps.jl +++ b/src/DynamicsMethods/IntegrationAlgorithms/steps.jl @@ -43,7 +43,7 @@ function step_O!(cache, integrator) @unpack dutmp, flatdutmp, tmp1, tmp2, gtmp, noise, half, c1, c2 = cache Λ = gtmp - σ = sqrt(get_temperature(p, t+dt*half)) ./ sqrt.(repeat(p.atoms.masses; inner=ndofs(p))) + σ = repeat( sqrt.(get_temperature(p, t+dt*half) ./ p.atoms.masses);inner=ndofs(p)) @.. noise = σ*W.dW[:] / sqdt @@ -77,7 +77,7 @@ function step_O!(friction::MDEFCache, integrator, v, r, t) integrator.g(gtmp,r,p,t) Λ = gtmp - @.. σ = sqrt(get_ring_polymer_temperature(p, t)) * sqrtmass + σ = repeat(sqrt.(get_ring_polymer_temperature(p, t)) .* sqrtmass;inner=ndofs(p)) @views for i in axes(r, 3) @.. noise = σ * W.dW[:,:,i][:] / sqdt @@ -110,7 +110,7 @@ function step_O!(friction::LangevinCache, integrator, v::RingPolymerArray, r::Ri @unpack W, p, dt, sqdt = integrator @unpack c1, c2, sqrtmass, σ = friction - @.. σ = sqrt(get_ring_polymer_temperature(p, t)) * sqrtmass + σ = repeat(sqrt.(get_ring_polymer_temperature(p, t)) .* sqrtmass;inner=ndofs(p)) for i in axes(r, 3) for j in RingPolymerArrays.quantumindices(v) diff --git a/src/DynamicsOutputs.jl b/src/DynamicsOutputs.jl index 0d851f0f5..a293d2634 100644 --- a/src/DynamicsOutputs.jl +++ b/src/DynamicsOutputs.jl @@ -7,15 +7,16 @@ Defines a set of functions that can be used to calculate outputs for dynamics si module DynamicsOutputs using RingPolymerArrays: get_centroid -using Unitful: @u_str -using UnitfulAtomic: austrip +using Unitful +using UnitfulAtomic using LinearAlgebra: norm using ComponentArrays: ComponentVector using NQCDynamics: Estimators, DynamicsUtils, - Analysis + Analysis, + ndofs using NQCModels: NQCModels @@ -40,7 +41,7 @@ export OutputVelocity """ OutputCentroidVelocity(sol, i) -Output the velocity of the ring polymer centroid at each timestep during the trajectory. +Output the velocity of the ring polymer centroid at each timestep during the trajectory. """ OutputCentroidVelocity(sol, i) = [get_centroid(get_velocities(u)) for u in sol.u] export OutputCentroidVelocity @@ -96,10 +97,13 @@ Evaluate the classical kinetic energy at each timestep during the trajectory. OutputKineticEnergy(sol, i) = DynamicsUtils.classical_kinetic_energy.(sol.prob.p, sol.u) export OutputKineticEnergy +OutputFinalKineticEnergy(sol, i) = DynamicsUtils.classical_kinetic_energy(sol.prob.p, last(sol.u)) +export OutputFinalKineticEnergy + """ OutputSubsetKineticEnergy(sol, i) -Evaluate the classical kinetic energy of a subset of the entire system at each save step. +Evaluate the classical kinetic energy of a subset of the entire system at each save step. The subset is defined by `OutputSubsetKineticEnergy(indices)`. """ @@ -107,10 +111,23 @@ struct OutputSubsetKineticEnergy{T} indices::T end function (output::OutputSubsetKineticEnergy)(sol, i) - return map(x->DynamicsUtils.classical_kinetic_energy(sol.prob.p.atoms.masses[output.indices], x[:, output.indices]), [DynamicsUtils.get_velocities(i) for i in sol.u]) + return map(x -> DynamicsUtils.classical_kinetic_energy(sol.prob.p.atoms.masses[output.indices], x[:, output.indices]), [DynamicsUtils.get_velocities(i) for i in sol.u]) end export OutputSubsetKineticEnergy +""" +Evaluate the classical kinetic energy of a subset of the entire system at the end of the simulation. + +The subset is defined by `OutputSubsetKineticEnergy(indices)`. +""" +struct OutputFinalSubsetKineticEnergy + indices +end +function (output::OutputFinalSubsetKineticEnergy)(sol, i) + return DynamicsUtils.classical_kinetic_energy(sol.prob.p.atoms.masses[output.indices], DynamicsUtils.get_velocities(last(sol.u))[:, output.indices]) +end +export OutputFinalSubsetKineticEnergy + OutputSpringEnergy(sol, i) = DynamicsUtils.classical_spring_energy.(sol.prob.p, sol.u) export OutputSpringEnergy @@ -204,6 +221,13 @@ Output the end point of each trajectory. OutputFinal(sol, i) = last(sol.u) export OutputFinal +""" +Outputs the final time point of a trajectory. This is useful if simulations are terminated by a callback. +""" +OutputFinalTime(sol, i) = last(sol.t) +export OutputFinalTime + + """ Output the total number of surface hops during the trajectory """ @@ -232,24 +256,49 @@ export OutputDissociation function (output::OutputDissociation)(sol, i) R = DynamicsUtils.get_positions(last(sol.u)) - dissociated = norm(R[:,output.atom_indices[1]] .- R[:,output.atom_indices[2]]) > output.distance + dissociated = norm(R[:, output.atom_indices[1]] .- R[:, output.atom_indices[2]]) > output.distance return dissociated ? 1 : 0 end +""" +Outputs a 1 if the molecule has moved a certain distance from all atoms in the surface, 0 otherwise. + +Distance is defined by the projected distance from the highest non-adsorbate atom along the surface normal. +""" +struct OutputSurfaceDesorption + distance::Number + adsorbate_indices::Vector{Int} + surface_normal::AbstractVector + OutputSurfaceDesorption(distance, adsorbate_indices; surface_normal=[0, 0, 1]) = new(distance, adsorbate_indices, surface_normal) +end + +function (osd::OutputSurfaceDesorption)(sol, i) + desorption_frame = findfirst([Analysis.Diatomic.surface_distance_condition( + i, + osd.adsorbate_indices, + sol.prob.p; + surface_distance_threshold=osd.distance, + surface_normal=osd.surface_normal, + ) for i in sol.u] + ) + return desorption_frame === nothing ? 0 : 1 +end + +export OutputSurfaceDesorption + """ Output the vibrational and rotational quantum numbers of the final image. """ -struct OutputQuantisedDiatomic{S,H,V} - sim::S +struct OutputQuantisedDiatomic{H,V} height::H normal_vector::V end -OutputQuantisedDiatomic(sim; height=10, normal_vector=[0, 0, 1]) = OutputQuantisedDiatomic(sim, height, normal_vector) +OutputQuantisedDiatomic(; height=10, normal_vector=[0, 0, 1]) = OutputQuantisedDiatomic(height, normal_vector) export OutputQuantisedDiatomic function (output::OutputQuantisedDiatomic)(sol, i) - final = last(sol.u) - ν, J = QuantisedDiatomic.quantise_diatomic(output.sim, + final = last(sol.u) + ν, J = QuantisedDiatomic.quantise_diatomic(sol.prob.p, DynamicsUtils.get_velocities(final), DynamicsUtils.get_positions(final); height=output.height, normal_vector=output.normal_vector) return (ν, J) @@ -279,7 +328,7 @@ function (output::OutputStateResolvedScattering1D)(sol, i) transmission=zeros(NQCModels.nstates(output.sim)) ) x = DynamicsUtils.get_positions(final)[1] - if x > 0 # If final position past 0 then we count as transmission + if x > 0 # If final position past 0 then we count as transmission output.transmission .= populations else # If final position left of 0 then we count as reflection output.reflection .= populations @@ -297,15 +346,15 @@ struct OutputDesorptionAngle{I<:Vector{Int},S<:Vector{Float64},D} surface_distance_threshold::D end """ - OutputDesorptionAngle(indices; surface_normal = [0,0,1], surface_distance_threshold = 5.0u"Å") + OutputDesorptionAngle(indices; surface_normal=[0, 0, 1], surface_distance_threshold=austrip(5.0u"Å")) -Outputs the desorption angle in degrees (relative to the surface normal) if a desorption event is detected. +Outputs the desorption angle in degrees (relative to the surface normal) if a desorption event is detected. Use `surface_normal` to define the direction "away" from the surface. Most commonly, this would be in positive z direction. -A desorption is detected if the centre of mass of the molecule defined with `indices` is above `surface_distance_threshold` from the closest surface atom. -This is calculated with respect to `surface_normal` and will take into account periodic boundary conditions. +A desorption is detected if the centre of mass of the molecule defined with `indices` is above `surface_distance_threshold` from the closest surface atom. +This is calculated with respect to `surface_normal` and will take into account periodic boundary conditions. """ -OutputDesorptionAngle(indices; surface_normal = [0,0,1], surface_distance_threshold = 5.0u"Å") = OutputDesorptionAngle(indices, convert(Vector{Float64},surface_normal), surface_distance_threshold) +OutputDesorptionAngle(indices; surface_normal=[0, 0, 1], surface_distance_threshold=austrip(5.0u"Å")) = OutputDesorptionAngle(indices, convert(Vector{Float64}, surface_normal), surface_distance_threshold) export OutputDesorptionAngle """ @@ -317,25 +366,25 @@ function (output::OutputDesorptionAngle)(sol, i) return Analysis.Diatomic.get_desorption_angle(sol.u, output.indices, sol.prob.p; surface_normal=output.surface_normal, surface_distance_threshold=output.surface_distance_threshold) end -struct OutputDesorptionTrajectory{I<:Vector{Int}, N<:Vector{Float64}, D, F<:Int} +struct OutputDesorptionTrajectory{I<:Vector{Int},N<:Vector{Float64},D,F<:Int} indices::I surface_normal::N surface_distance_threshold::D extra_frames::F end """ - `OutputDesorptionTrajectory(indices; surface_normal = [0,0,1], surface_distance_threshold = 5.0u"Å", extra_frames = 0)` + OutputDesorptionTrajectory(indices; surface_normal=[0, 0, 1], surface_distance_threshold=austrip(5.0u"Å"), extra_frames=0) Like OutputDynamicsVariables, but only saves parts of the trajectory where desorption is occurring. -Use `surface_normal` to define the direction "away" from the surface. Most commonly, this would be in positive z direction. +Use `surface_normal` to define the direction "away" from the surface. Most commonly, this would be in positive z direction. -Use `extra_frames` to save additional steps before the desorption event begins. +Use `extra_frames` to save additional steps before the desorption event begins. -A desorption is detected if the centre of mass of the molecule defined with `indices` is above `surface_distance_threshold` from the closest surface atom. -This is calculated with respect to `surface_normal` and will take into account periodic boundary conditions. +A desorption is detected if the centre of mass of the molecule defined with `indices` is above `surface_distance_threshold` from the closest surface atom. +This is calculated with respect to `surface_normal` and will take into account periodic boundary conditions. """ -OutputDesorptionTrajectory(indices; surface_normal = [0,0,1], surface_distance_threshold = 5.0u"Å", extra_frames = 0) = OutputDesorptionTrajectory(indices, convert(Vector{Float64},surface_normal), surface_distance_threshold, extra_frames) +OutputDesorptionTrajectory(indices; surface_normal=[0, 0, 1], surface_distance_threshold=austrip(5.0u"Å"), extra_frames=0) = OutputDesorptionTrajectory(indices, convert(Vector{Float64}, surface_normal), surface_distance_threshold, extra_frames) """ (output::OutputDesorptionTrajectory)(sol, i) @@ -343,8 +392,104 @@ Only output parts of the trajectory where desorption is occurring. """ function (output::OutputDesorptionTrajectory)(sol, i) desorption_frame = Analysis.Diatomic.get_desorption_frame(sol.u, output.indices, sol.prob.p; surface_distance_threshold=output.surface_distance_threshold, surface_normal=output.surface_normal) - return isa(desorption_frame, Int) ? sol.u[desorption_frame-output.extra_frames:end] : nothing + if isnothing(desorption_frame) + return nothing + end + start_save_frame = desorption_frame - output.extra_frames + if start_save_frame < 0 + return sol.u + else + return sol.u[start_save_frame:desorption_frame] + end end export OutputDesorptionTrajectory +struct OutputDesorptionSnapshot{I<:Vector{Int},N<:Vector{Float64},D} + indices::I + surface_normal::N + surface_distance_threshold::D +end +""" + OutputDesorptionSnapshot(indices; surface_normal=[0, 0, 1], surface_distance_threshold=austrip(5.0u"Å")) + +Save DynamicsVariables where desorption starts. (Same conditions as for OutputDesorptionTrajectory) + +Use `surface_normal` to define the direction "away" from the surface. Most commonly, this would be in positive z direction. + +Use `extra_frames` to save additional steps before the desorption event begins. + +A desorption is detected if the centre of mass of the molecule defined with `indices` is above `surface_distance_threshold` from the closest surface atom. +This is calculated with respect to `surface_normal` and will take into account periodic boundary conditions. +""" +OutputDesorptionSnapshot(indices; surface_normal=[0, 0, 1], surface_distance_threshold=austrip(5.0u"Å")) = OutputDesorptionSnapshot(indices, convert(Vector{Float64}, surface_normal), surface_distance_threshold) +""" + (output::OutputDesorptionTrajectory)(sol, i) + +Only output trajectory snapshot where desorption begins. (Centre of mass velocity projected onto surface +normal changes sign) +""" +function (output::OutputDesorptionSnapshot)(sol, i) + desorption_frame = Analysis.Diatomic.get_desorption_frame(sol.u, output.indices, sol.prob.p; surface_distance_threshold=output.surface_distance_threshold, surface_normal=output.surface_normal) + return isa(desorption_frame, Int) ? sol.u[desorption_frame] : nothing +end +export OutputDesorptionSnapshot + +""" +Outputs the instantaneous temperature of the selected atoms **in K** in the system at every save point. + +Invoke with `OutputKineticTemperature(:)` for the entire system, or with `OutputKineticTemperature([1,2,3...])` for a subset of atoms. + +""" +struct OutputKineticTemperature + indices + OutputKineticTemperature(indices) = new(indices) +end + +function (out::OutputKineticTemperature)(sol, i) + # Allocate output vector + kinetic_energies = zeros(typeof(DynamicsUtils.classical_kinetic_energy(sol.prob.p, sol.u[1])), length(sol.u)) + # Determine number of atoms + if isa(out.indices, Colon) + n_atoms = length(sol.prob.p.atoms.masses) + elseif isa(out.indices, UnitRange) + n_atoms = length(collect(out.indices)) + else + n_atoms = length(out.indices) + end + # Calculate kinetic temperatures. + for snapshot in eachindex(sol.u) + kinetic_energy = DynamicsUtils.classical_kinetic_energy(sol.prob.p.atoms.masses[out.indices], DynamicsUtils.get_velocities(sol.u[snapshot])[:, out.indices]) + kinetic_energies[snapshot] = ustrip(uconvert(u"K", 2 * kinetic_energy * u"hartree/k_au") / ndofs(sol.prob.p) / n_atoms) + end + return kinetic_energies +end + +export OutputKineticTemperature + +""" + OutputNoise(sol, i) + +Outputs the noise generated by the integrator at each time step in the trajectory. + +**Note:** This requires using the `run_dynamics` command with `save_noise=true`, otherwise noise won't be accessible to this output function. +""" +function OutputNoise(sol, i) + return sol.W +end + +export OutputNoise + +""" + OutputEverything(sol,i) + +Outputs the full DifferentialEquations solution object. + +Storing this to disk is inefficient, but allows for full post-processing with +any of the functions defined in this module. +""" +OutputEverything(sol, i) = sol +OutputSolution(sol, i) = sol +export OutputEverything +export OutputSolution + end # module diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 812a39979..f7381b40d 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -21,7 +21,8 @@ module QuantisedDiatomic using QuadGK: QuadGK using Rotations: Rotations using LinearAlgebra: norm, normalize!, nullspace, cross, I -using UnitfulAtomic: austrip +using Unitful +using UnitfulAtomic using Distributions: Uniform using UnicodePlots: lineplot, lineplot!, DotCanvas, histogram using Optim: Optim @@ -37,6 +38,7 @@ using NQCModels.AdiabaticModels: AdiabaticModel export generate_configurations export quantise_diatomic +export energy_distribution_from_quantisation const TIMER = TimerOutputs.TimerOutput() @@ -49,11 +51,13 @@ struct EffectivePotential{JType,T,B,F} μ::T J::JType binding_curve::BindingCurve{T,B,F} + langer_modification::Bool end function (effective_potential::EffectivePotential)(r) - (; μ, J, binding_curve) = effective_potential - rotational = J * (J + 1) / (2μ * r^2) + (; μ, J, binding_curve, langer_modification) = effective_potential + L_squared = langer_modification ? (J + 1 / 2)^2 : J * (J + 1) + rotational = L_squared / (2μ * r^2) potential = binding_curve.fit(r) return rotational + potential @@ -92,7 +96,8 @@ function generate_configurations(sim, ν, J; direction=[0, 0, -1.0], atom_indices=[1, 2], r=zeros(size(sim)), - bond_lengths=0.5:0.01:5.0 + bond_lengths=0.5:0.01:5.0, + langer_modification=false ) μ = reduced_mass(masses(sim)[atom_indices]) @@ -103,7 +108,7 @@ function generate_configurations(sim, ν, J; binding_curve = calculate_binding_curve(bond_lengths, sim.calculator.model, environment) plot_binding_curve(binding_curve.bond_lengths, binding_curve.potential, binding_curve.fit) - V = EffectivePotential(μ, J, binding_curve) + V = EffectivePotential(μ, J, binding_curve, langer_modification) total_energy, bounds = find_total_energy(V, ν) @@ -118,7 +123,7 @@ function generate_configurations(sim, ν, J; end function generate_1D_vibrations(model::AdiabaticModel, μ::Real, ν::Integer; - samples=1000, bond_lengths=0.5:0.01:5.0 + samples=1000, bond_lengths=0.5:0.01:5.0, langer_modification=false ) J = 0 @@ -126,7 +131,7 @@ function generate_1D_vibrations(model::AdiabaticModel, μ::Real, ν::Integer; binding_curve = calculate_binding_curve(bond_lengths, model, environment) plot_binding_curve(binding_curve.bond_lengths, binding_curve.potential, binding_curve.fit) - V = EffectivePotential(μ, J, binding_curve) + V = EffectivePotential(μ, J, binding_curve, langer_modification) total_energy, bounds = find_total_energy(V, ν) @@ -253,14 +258,14 @@ function find_integral_bounds(total_energy::Real, V::EffectivePotential) return r₁, r₂ end -function binding_curve_from_structure(sim::Simulation, v::Matrix, r::Matrix; - bond_lengths=0.5:0.01:5.0, - height=10.0, - surface_normal=[0, 0, 1.0], - atom_indices=[1, 2], +function binding_curve_from_structure(sim::Simulation, v::Matrix, r::Matrix; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], args... ) - # Separate slab from molecule if necessary. + # Separate slab from molecule if necessary. r, slab = separate_slab_and_molecule(atom_indices, r) v, slab_v = separate_slab_and_molecule(atom_indices, v) environment = EvaluationEnvironment(atom_indices, size(sim), slab, austrip(height), surface_normal) @@ -272,11 +277,11 @@ function binding_curve_from_structure(sim::Simulation, v::Matrix, r::Matrix; end """ - quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; - bond_lengths=0.5:0.01:5.0, - height=10.0, - surface_normal=[0, 0, 1.0], - atom_indices=[1, 2], + quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], max_translation=1, show_timer=false, reset_timer=false ) @@ -284,120 +289,126 @@ end Quantise the vibrational and rotational degrees of freedom for the specified positions and velocities. -If the potential can be evaluated for the diatomic only, independent of position, -supplying a `Simulation` for just the diatomic will speed up evaluation. +If the potential can be evaluated for the diatomic only, independent of position, +supplying a `Simulation` for just the diatomic will speed up evaluation. When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`. If the potential is independent of centre of mass position, this has no effect. Otherwise, be sure to modify these parameters to give the intended behaviour. -If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms -will be used if positions are close to cell boundaries. -Set `max_translation` to the radius of surrounding unit cells to search. +If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms +will be used if positions are close to cell boundaries. +Set `max_translation` to the radius of surrounding unit cells to search. (e.g. 1 if positions are already wrapped around cell boundaries) """ -function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; - bond_lengths=0.5:0.01:5.0, - height=10.0, - surface_normal=[0, 0, 1.0], - atom_indices=[1, 2], +function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + output_energies=false, args... ) - binding_curve=binding_curve_from_structure(sim, v, r; bond_lengths=bond_lengths, height=height, surface_normal=surface_normal, atom_indices=atom_indices) - return quantise_diatomic(sim, v, r, binding_curve; height=height, surface_normal=surface_normal, atom_indices=atom_indices, args...) + binding_curve = binding_curve_from_structure(sim, v, r; bond_lengths=bond_lengths, height=height, surface_normal=surface_normal, atom_indices=atom_indices) + return quantise_diatomic(sim, v, r, binding_curve; height=height, surface_normal=surface_normal, atom_indices=atom_indices, output_energies=output_energies, args...) end """ quantise_diatomic(sim::Simulation, v::Vector{Matrix}, r::Vector{Matrix}; - bond_lengths=0.5:0.01:5.0, - height=10.0, - surface_normal=[0, 0, 1.0], - atom_indices=[1, 2], + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], show_timer=false, reset_timer=false ) -Quantise the vibrational and rotational degrees of freedom of multiple atomic configurations -given as a vector of velocity matrices and a vector of position matrices. +Quantise the vibrational and rotational degrees of freedom of multiple atomic configurations +given as a vector of velocity matrices and a vector of position matrices. -If the potential can be evaluated for the diatomic only, independent of position, -supplying a `Simulation` for just the diatomic will speed up evaluation. +If the potential can be evaluated for the diatomic only, independent of position, +supplying a `Simulation` for just the diatomic will speed up evaluation. When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`. If the potential is independent of centre of mass position, this has no effect. Otherwise, be sure to modify these parameters to give the intended behaviour. -If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms -will be used if positions are close to cell boundaries. -Set `max_translation` to the radius of surrounding unit cells to search. +If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms +will be used if positions are close to cell boundaries. +Set `max_translation` to the radius of surrounding unit cells to search. (e.g. 1 if positions are already wrapped around cell boundaries) Specify `show_timer=true` for performance timings of the EBK quantisation process and -`reset_timer=true` to see timings for each individual quantisation. +`reset_timer=true` to see timings for each individual quantisation. """ -function quantise_diatomic(sim::Simulation, v::Vector{<: Matrix{<: Any}}, r::Vector{<: Matrix{<: Any}}; - bond_lengths=0.5:0.01:5.0, - height=10.0, - surface_normal=[0, 0, 1.0], - atom_indices=[1, 2], +function quantise_diatomic(sim::Simulation, v::Vector{<:Matrix{<:Any}}, r::Vector{<:Matrix{<:Any}}; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + output_energies=false, args... ) # Generate binding curve for all structures - binding_curve=binding_curve_from_structure(sim, v[1], r[1]; bond_lengths=bond_lengths, height=height, surface_normal=surface_normal, atom_indices=atom_indices) - ν, J=[],[] + binding_curve = binding_curve_from_structure(sim, v[1], r[1]; bond_lengths=bond_lengths, height=height, surface_normal=surface_normal, atom_indices=atom_indices) + results = [] - # Quantise all configurations in the given vectors. + # Quantise all configurations in the given vectors. @showprogress for index in eachindex(v) try - result=quantise_diatomic(sim, v[index], r[index], binding_curve; height=height, surface_normal=surface_normal, atom_indices=atom_indices, args...) - push!(ν, result[1]) - push!(J, result[2]) + result = quantise_diatomic(sim, v[index], r[index], binding_curve; height=height, surface_normal=surface_normal, atom_indices=atom_indices, output_energies=output_energies, args...) + push!(results, result) catch e @warn "Quantisation was unsuccessful for configuration at $(index).\nThe final results vectors will show ν=-1 and J=-1 for this configuration.\n The error was: $(e)" - push!(ν, -1) - push!(J, -1) + output = output_energies ? (missing, missing, missing, missing, missing) : (missing, missing) + push!(results, output) end end - return ν, J + return results # returns (ν, J) or (ν, J, translation_energy, rotation_energy, vibration_energy) end """ quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve::BindingCurve; show_timer=false, reset_timer=false, - height=10, normal_vector=[0, 0, 1], atom_indices=[1,2], max_translation=1) + height=10, normal_vector=[0, 0, 1], atom_indices=[1,2], max_translation=1) ) Quantise the vibrational and rotational degrees of freedom for the specified -positions and velocities using the `BindingCurve` specified. -A binding curve will be automatically generated if you do not supply one. +positions and velocities using the `BindingCurve` specified. +A binding curve will be automatically generated if you do not supply one. -If the potential can be evaluated for the diatomic only, independent of position, -supplying a `Simulation` for just the diatomic will speed up evaluation. +If the potential can be evaluated for the diatomic only, independent of position, +supplying a `Simulation` for just the diatomic will speed up evaluation. When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`. If the potential is independent of centre of mass position, this has no effect. Otherwise, be sure to modify these parameters to give the intended behaviour. -If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms -will be used if positions are close to cell boundaries. -Set `max_translation` to the radius of surrounding unit cells to search. +If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms +will be used if positions are close to cell boundaries. +Set `max_translation` to the radius of surrounding unit cells to search. (e.g. 1 if positions are already wrapped around cell boundaries) """ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve::BindingCurve; - show_timer=false, reset_timer=false, - max_translation=1, height=10.0, - surface_normal=[0, 0, 1.0], atom_indices=[1, 2], + show_timer=false, + reset_timer=false, + max_translation=1, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + langer_modification=false, # The Langer modification uses L^2 = J*(J+1/2)^2ħ^2 instead of L^2 = J*(J+1)ħ^2 + output_energies=false, # Output translation, rotation and vibrational energies with the quantisation information. ) reset_timer && TimerOutputs.reset_timer!(TIMER) - if isa(sim.cell,PeriodicCell) + if isa(sim.cell, PeriodicCell) # If the simulation used a `PeriodicCell`, translate `atom_indices` so they are at their minimum distance. (This is necessary if atoms were translated back into the original unit cell) - translations=[[i,j,k] for i in -max_translation:max_translation for j in -max_translation:max_translation for k in -max_translation:max_translation] - which_translation=argmin([norm(abs.(r[:,atom_indices[2]]-r[:,atom_indices[1]]+sim.cell.vectors*operation)) for operation in translations]) - # Translate one atom for minimal distance. - if translations[which_translation]!=[0,0,0] - r[:,atom_indices[2]].=r[:,atom_indices[2]]+sim.cell.vectors*translations[which_translation] - @debug "Using a periodic copy of atom "*string(atom_indices[end])*" to bring it closer to atom "*string(atom_indices[begin]) + translations = [[i, j, k] for i in -max_translation:max_translation for j in -max_translation:max_translation for k in -max_translation:max_translation] + which_translation = argmin([norm(abs.(r[:, atom_indices[2]] - r[:, atom_indices[1]] + sim.cell.vectors * operation)) for operation in translations]) + # Translate one atom for minimal distance. + if translations[which_translation] != [0, 0, 0] + r[:, atom_indices[2]] .= r[:, atom_indices[2]] + sim.cell.vectors * translations[which_translation] + @debug "Using a periodic copy of atom " * string(atom_indices[end]) * " to bring it closer to atom " * string(atom_indices[begin]) end end @@ -405,24 +416,34 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve: v, slab_v = separate_slab_and_molecule(atom_indices, v) environment = EvaluationEnvironment(atom_indices, size(sim), slab, austrip(height), surface_normal) - @debug "After PBC check and separation from slab, diatomic positions are:\n$(r)" - r_com = subtract_centre_of_mass(r, masses(sim)[atom_indices]) v_com = subtract_centre_of_mass(v, masses(sim)[atom_indices]) p_com = v_com .* masses(sim)[atom_indices]' k = DynamicsUtils.classical_kinetic_energy(masses(sim)[atom_indices], v_com) + @debug "Diatomic Kinetic energy: $(auconvert(u"eV", k))" p = calculate_diatomic_energy(bond_length(r_com), sim.calculator.model, environment) + @debug "Diatomic Potential energy: $(auconvert(u"eV", p))" E = k + p L = total_angular_momentum(r_com, p_com) - J = (sqrt(1+4*L^2) - 1) / 2 # L^2 = J(J+1)ħ^2 + @debug "Total angular momentum: $(L)" + J = langer_modification ? L - (1 / 2) : (sqrt(1 + 4 * L^2) - 1) / 2 # L^2 = J(J+1)ħ^2 + @debug "Calculated J (before rounding): $(J)" μ = reduced_mass(masses(sim)[atom_indices]) - @debug "k=$k\np=$p\nE=$E\nL=$L\nJ=$J\n" + V = EffectivePotential(μ, J, binding_curve, langer_modification) - V = EffectivePotential(μ, J, binding_curve) + @debug begin + diag_plt = lineplot( + binding_curve.bond_lengths, + sqrt.(abs.(E .- V.(binding_curve.bond_lengths))); + title="Binding curve", xlabel="Bond length / bohr", ylabel="Energy / Hartree", + name="E - V values (should intersect 0 twice)", canvas=DotCanvas, border=:ascii + ) + show(diag_plt) + end r₁, r₂ = @timeit TIMER "Finding bounds" find_integral_bounds(E, V) @@ -437,12 +458,22 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve: TimerOutputs.complement!(TIMER) show_timer && show(TIMER) - @debug "Found ν=$ν" - return round(Int, ν), round(Int, J) + @debug "Calculated ν: $ν" + if !output_energies + return round(Int, ν), round(Int, J) + else + ν = round(Int, ν) + J = round(Int, J) + translation_energy = sum(sim.atoms.masses[atom_indices]) / 2 * norm(centre_of_mass(v[:, atom_indices], sim.atoms.masses[atom_indices]))^2 # m/2 v^2 + rotation_energy = J * (J + 1) / (2 * sum(r_com .^ 2 .* masses(sim)[atom_indices]')) # $E = \frac{J*(J+1)}{2I}$ in a.u. + vibration_energy = (ν + 0.5) * sqrt(calculate_force_constant(binding_curve) / μ) # $(ν+1/2)\sqrt{\frac{k}{μ}}$ + + return (round(Int, ν), round(Int, J), translation_energy, rotation_energy, vibration_energy) + end end function quantise_1D_vibration(model::AdiabaticModel, μ::Real, r::Real, v::Real; - bond_lengths=0.5:0.01:5.0, reset_timer=false, show_timer=false + bond_lengths=0.5:0.01:5.0, reset_timer=false, show_timer=false, langer_modification=false ) reset_timer && TimerOutputs.reset_timer!(TIMER) @@ -455,7 +486,7 @@ function quantise_1D_vibration(model::AdiabaticModel, μ::Real, r::Real, v::Real binding_curve = calculate_binding_curve(bond_lengths, model, environment) - V = EffectivePotential(μ, J, binding_curve) + V = EffectivePotential(μ, J, binding_curve, langer_modification) r₁, r₂ = @timeit TIMER "Finding bounds" find_integral_bounds(E, V) diff --git a/src/InitialConditions/QuantisedDiatomic/energy_evaluation.jl b/src/InitialConditions/QuantisedDiatomic/energy_evaluation.jl index 3ae74ec2c..6dfc7fff8 100644 --- a/src/InitialConditions/QuantisedDiatomic/energy_evaluation.jl +++ b/src/InitialConditions/QuantisedDiatomic/energy_evaluation.jl @@ -10,8 +10,10 @@ end Base.broadcastable(p::EvaluationEnvironment) = Ref(p) function EvaluationEnvironment(molecule_indices, system_size, other_atoms, height, surface_normal) + normalize!(surface_normal) # Normalise surface normal vector + # If surface atoms are explicitly included, find the height of the surface along `surface_normal` and adjust `height` accordingly. if size(other_atoms, 2) != 0 - height += maximum(other_atoms[3,:]) + height += reduce(max, mapslices(norm, other_atoms, dims=1)) end offset = surface_normal .* height orthogonal_vectors = nullspace(reshape(surface_normal, 1, :)) # Plane parallel to surface @@ -49,9 +51,10 @@ end Returns potential energy of diatomic with `bond_length` at `height` from surface. -Orients molecule parallel to the surface at the specified height, the surface is assumed -to intersect the origin. -This requires that the model implicitly provides the surface, or works fine without one. +Orients molecule parallel to the surface at the specified height within the simulation cell, +assuming the height has already been adjusted to include that of the surface. + +(this is checked in the EvaluationEnvironment constructor) """ function calculate_diatomic_energy( bond_length::Real, model::AdiabaticModel, environment::EvaluationEnvironment diff --git a/src/NQCDistributions-convenience.jl b/src/NQCDistributions-convenience.jl new file mode 100644 index 000000000..1f4423edb --- /dev/null +++ b/src/NQCDistributions-convenience.jl @@ -0,0 +1,28 @@ +""" +Add-on functions for NQCDistributions to make generating distributions more convenient. + +""" + +import Distributions + +""" + VelocityBoltzmann(temperature, sim::AbstractSimulation; center = zeros(size(sim))) + +Generates a VelocityBoltzmann distribution covering an entire Simulation. +If `NQCModels.mobileatoms` is modified, velocities for immobile atoms will always be zero. +""" +function NQCDistributions.VelocityBoltzmann(temperature, sim::AbstractSimulation; center = zeros(Float64, size(sim))) + if length(NQCModels.mobileatoms(sim))==size(sim)[2] + # Simple case: All atoms should be moving, so give them all Boltzmann distributions. + return NQCDistributions.VelocityBoltzmann(temperature, sim.atoms.masses, size(sim); center=center) + else + # Not all atoms should be moving, so assign Boltzmann distribution only to mobile atoms. + sampleable_array=convert(Matrix{Distributions.UnivariateDistribution}, hcat([[Distributions.Dirac(center[i,j]) for i in 1:size(sim)[1]] for j in 1:size(sim)[2]]...)) + for i in NQCModels.mobileatoms(sim) + for j in NQCModels.dofs(sim) + sampleable_array[j,i]=NQCDistributions.VelocityBoltzmann(temperature, sim.atoms.masses[i]; center=center[j,i]) + end + end + return NQCDistributions.UnivariateArray(sampleable_array) + end +end diff --git a/src/NQCDynamics.jl b/src/NQCDynamics.jl index c35e6a7e3..35c35dd11 100644 --- a/src/NQCDynamics.jl +++ b/src/NQCDynamics.jl @@ -17,9 +17,13 @@ include("simulations.jl") export Simulation, RingPolymerSimulation, natoms, - masses + masses, + Thermostat, + get_temperature @reexport using NQCDistributions +# Simulation-aware version of nuclear Boltzmann distribution. +include("NQCDistributions-convenience.jl") include("DynamicsUtils/DynamicsUtils.jl") @reexport using .DynamicsUtils: get_positions, get_velocities diff --git a/src/simulations.jl b/src/simulations.jl index 2b2f2014e..6eddc8d13 100644 --- a/src/simulations.jl +++ b/src/simulations.jl @@ -3,7 +3,7 @@ using Unitful: @u_str using UnitfulAtomic: austrip, auconvert using .Calculators: AbstractCalculator, Calculator -using NQCModels: Model +using NQCModels: Model, Subsystem, CompositeModel abstract type AbstractSimulation{M} end @@ -27,6 +27,18 @@ dynamics method, temperature and simulation cell. function Simulation(atoms::Atoms{T}, model::Model, method::M; temperature=0u"K", cell::AbstractCell=InfiniteCell()) where {M,T} calc = Calculator(model, length(atoms), T) + # If a thermostat is provided, check it covers the whole system. + isa(temperature, Thermostat{Vector{Int}}) ? throw(DomainError(temperature, "Thermostat must apply to all atoms.")) : nothing + # If multiple Thermostats are provided, check that each atom only has one thermostat applied to it. + if isa(temperature, Vector{<:Thermostat}) + indices = vcat([thermostat.indices for thermostat in temperature]...) + if length(unique(indices)) != length(atoms.masses) + throw(DomainError(temperature, "Every atom must have a Thermostat applied to it.")) + end + if length(indices) != length(unique(indices)) + throw(DomainError(temperature, "Atoms can only have one thermostat applied to them.")) + end + end Simulation(temperature, cell, atoms, calc, method) end @@ -88,7 +100,7 @@ get_temperature(temperature::Number, t=0) = austrip(temperature) get_temperature(temperature::Function, t=0) = austrip(temperature(t)) function get_ring_polymer_temperature(sim::RingPolymerSimulation, t::Real=0) - return get_temperature(sim, t) * nbeads(sim) + return @. get_temperature(sim, t) * nbeads(sim) end function Base.show(io::IO, sim::Simulation{M}) where {M} @@ -99,3 +111,64 @@ function Base.show(io::IO, sim::RingPolymerSimulation{M}) where {M} print(io, "RingPolymerSimulation{$M}:\n\n ", sim.atoms, "\n\n ", sim.calculator.model, "\n with ", length(sim.beads), " beads.") end + +""" +A Temperature contains both temperature information and the atom indices within a Simulation that it is applied to. + +**If you don't need to apply different temperatures to different parts of your Simulation, assign the `temperature` keyword argument of your simulation to a number or function.** + +## Parameters + +`value`: A temperature function. This can be a Number type for constant temperatures, or a function taking the time in Unitful `u"ps"` as input and giving a temperature in Unitful `u"K"` as output. + +`indices`: Indices of the atoms to which this thermostat is applied. Can be a range of indices, a single `Int`, or a `Vector{Int}`. + +""" +struct Temperature{I} + value::Function + indices::I +end + +function Base.show(io::IO, temperature::Temperature) + print(io, "Temperature:\n\tT(t=0) = $(get_temperature(temperature, 0u"fs"))\n\tApplies to atoms: $(temperature.indices)\n") +end + +function Temperature(value, indices=:) + if isa(indices, UnitRange) + indices=collect(indices) + elseif isa(indices, Int) + indices=[indices] + end + if !isa(value, Function) + temperature_function(t)=value + return Temperature(temperature_function, indices) + else + return Temperature(value, indices) + end +end + +""" + Temperature(temperature, subsystem::NQCModels.Subsystem) + +Apply a `Temperature` to all atoms in a `Subsystem`. +""" +function Temperature(temperature, subsystem::NQCModels.Subsystem) + Temperature(temperature, subsystem.indices) +end + +get_temperature(thermostat::Temperature{<:Any}, t=0u"fs") = austrip(thermostat.value(t)) + +""" + get_temperature(thermostats::Vector{<:Temperature}, t=0u"fs") + +Gets the temperature from multiple `Temperature`s and returns a vector of the temperature applied to each atom. +""" +function get_temperature(thermostats::Vector{<:Temperature}, t=0u"fs") + indices=vcat([thermostat.indices for thermostat in thermostats]...) + temperature_vector=zeros(Number, length(indices)) + for thermostat in thermostats + temperature_vector[thermostat.indices] .= austrip(thermostat.value(t)) + end + return temperature_vector +end + diff --git a/test/Core/simulations.jl b/test/Core/simulations.jl index 7f17d463c..fd759f0c5 100644 --- a/test/Core/simulations.jl +++ b/test/Core/simulations.jl @@ -1,6 +1,6 @@ using Test using NQCDynamics -using Unitful +using Unitful, UnitfulAtomic atoms = Atoms([:C, :H]) cell = InfiniteCell() @@ -22,7 +22,25 @@ model = NQCModels.Free() @test Simulation{Langevin}(atoms, model) isa Simulation{<:Langevin} @test Simulation{FSSH}(atoms, NQCModels.DoubleWell()) isa Simulation{<:FSSH} -@testset "get_temperature" begin +@testset "Thermostats" begin + # Init thermostat with Unitful quantity + thermostat1=Thermostat(10u"K", [1,2]) + @test NQCDynamics.get_temperature(thermostat1, 0) == NQCDynamics.get_temperature(thermostat1, 100) + # Init thermostat with function + thermostat2=Thermostat(x->(10+x*u"fs^-1")*u"K", [2,3]) + @test NQCDynamics.get_temperature(thermostat2) == austrip(10u"K") + @test NQCDynamics.get_temperature(thermostat2, 10u"fs") == austrip(20u"K") + # Test that thermostats with overlapping indices throw an error + @test_throws DomainError Simulation(atoms, model; temperature=[thermostat1, thermostat2]) + # Test system size / total thermostat size mismatch + @test_throws DomainError Simulation(Atoms([:N, :H, :H, :H,]), model; temperature=thermostat1) + thermostat3=thermostat2=Thermostat(x->(10+x*u"fs^-1")*u"K", [3,4]) + # Combined temperature evaluation + sim = Simulation(Atoms([:N, :H, :H, :H,]), model; temperature=[thermostat1, thermostat3]) + @test NQCDynamics.get_temperature(sim, austrip(1u"fs")) == austrip.([10u"K", 10u"K", 11u"K", 11u"K"]) +end + +@testset "get_temperature(Simulation)" begin sim = Simulation(atoms, model; temperature=10u"K") @test NQCDynamics.get_temperature(sim) isa Real sim = Simulation(atoms, model; temperature=x->10u"K")