diff --git a/Project.toml b/Project.toml index 0ca86e1f6..dca911953 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NQCDynamics" uuid = "36248dfb-79eb-4f4d-ab9c-e29ea5f33e14" authors = ["James "] -version = "0.13.3" +version = "0.13.4" [deps] AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d" @@ -45,16 +45,16 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" [compat] -AdvancedHMC = "0.5" -AdvancedMH = "0.6, 0.7" -ComponentArrays = "0.11, 0.12, 0.13, 0.14" +AdvancedHMC = "0.5, 0.6" +AdvancedMH = "0.6, 0.7, 0.8" +ComponentArrays = "0.11, 0.12, 0.13, 0.14, 0.15" DEDataArrays = "0.2" Dictionaries = "0.3" DiffEqBase = "6" Distributions = "0.25" FastBroadcast = "0.1, 0.2" FastLapackInterface = "1, 2" -HDF5 = "0.15, 0.16" +HDF5 = "0.15, 0.16, 0.17" Interpolations = "0.13, 0.14" MuladdMacro = "0.2" NQCBase = "0.2" @@ -72,7 +72,7 @@ Reexport = "1" RingPolymerArrays = "0.1" Roots = "1, 2" Rotations = "1" -SciMLBase = "1" +SciMLBase = "1, 2" StaticArrays = "1" StatsBase = "0.33, 0.34" StochasticDiffEq = "6" @@ -92,6 +92,7 @@ DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" JuLIP = "945c410c-986d-556a-acb1-167a618e0462" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" @@ -100,4 +101,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "SafeTestsets", "CSV", "DataFrames", "DiffEqNoiseProcess", "FiniteDiff", "DiffEqDevTools", "Logging", "PyCall", "JuLIP", "Symbolics", "Statistics", "Plots"] +test = ["Test", "SafeTestsets", "CSV", "DataFrames", "DiffEqNoiseProcess", "FiniteDiff", "DiffEqDevTools", "Logging", "MKL", "PyCall", "JuLIP", "Symbolics", "Statistics", "Plots"] diff --git a/docs/Project.toml b/docs/Project.toml index 974a2dbdd..1a502da20 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -26,5 +26,5 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a" [compat] -Documenter = "0.27" +Documenter = "1" DocumenterCitations = "1" diff --git a/docs/add_profiling_info.jl b/docs/add_profiling_info.jl new file mode 100644 index 000000000..7fdeb19e3 --- /dev/null +++ b/docs/add_profiling_info.jl @@ -0,0 +1,37 @@ + +function add_text_blocks(file) + preamble = """ + ```@setup logging + @info "Expanding $file..." + start_time = time() + ``` + """ + + postamble = """ + ```@setup logging + runtime = round(time() - start_time; digits=2) + @info "...done after \$runtime s." + ``` + """ + + (tmppath, tmpio) = mktemp() + open(file) do io + write(tmpio, preamble) + for line in eachline(io, keep=true) # keep so the new line isn't chomped + write(tmpio, line) + end + write(tmpio, postamble) + end + close(tmpio) + mv(tmppath, file, force=true) +end + + +for (root, dirs, files) in walkdir("src") + for file in files + if splitext(file)[2] == ".md" + filepath = joinpath(root, file) + add_text_blocks(filepath) + end + end +end diff --git a/docs/make.jl b/docs/make.jl index ea8565109..275610ffc 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -3,10 +3,6 @@ using DocumenterCitations using NQCBase, NQCModels, NQCDistributions, NQCDynamics using CubeLDFAModel, NNInterfaces -DocMeta.setdocmeta!(NQCDynamics, :DocTestSetup, :(using NQCDynamics); recursive=true) -DocMeta.setdocmeta!(NQCModels, :DocTestSetup, :(using NQCModels, Symbolics); recursive=true) -DocMeta.setdocmeta!(NQCBase, :DocTestSetup, :(using NQCBase); recursive=true) - bib = CitationBibliography(joinpath(@__DIR__, "references.bib")) function find_all_files(directory) @@ -16,16 +12,15 @@ function find_all_files(directory) ) end -@time makedocs( - bib, +@time makedocs(; + plugins=[bib], sitename="NQCDynamics.jl", modules=[NQCDynamics, NQCDistributions, NQCModels, NQCBase, CubeLDFAModel], - strict=true, + 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"], - ansicolor=true, ), authors="James Gardner and contributors.", pages=[ diff --git a/docs/src/NQCDistributions/overview.md b/docs/src/NQCDistributions/overview.md index 5e9a32027..f115063fc 100644 --- a/docs/src/NQCDistributions/overview.md +++ b/docs/src/NQCDistributions/overview.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/NQCDistributions/overview.md..." +start_time = time() +``` # NQCDistributions.jl # Storing and sampling distributions @@ -67,7 +71,7 @@ read on to the next sections about the included sampling methods. ### VelocityBoltzmann When performing equilibrium simulations it is often desirable to initialise trajectories -when thermal velocities, e.g. in combination with positions obtained from Monte Carlo sampling. +with thermal velocities, e.g. in combination with positions obtained from Monte Carlo sampling. These can be obtained for each atom from a gaussian distribution of the appropriate width, or alternatively, using the [`VelocityBoltzmann`](@ref) distribution which simplifies the process. @@ -155,3 +159,8 @@ total_dist.electronic.state # Returns the chosen electronic state. ``` + +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/NQCModels/analyticmodels.md b/docs/src/NQCModels/analyticmodels.md index 603bed907..df60c3b8e 100644 --- a/docs/src/NQCModels/analyticmodels.md +++ b/docs/src/NQCModels/analyticmodels.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/NQCModels/analyticmodels.md..." +start_time = time() +``` # Analytic model library This page plots many of the analytic models included in `NQCDynamics`. @@ -93,3 +97,7 @@ z = range(-0.5, 4.0, length=200) contour(x, z, v1, color=:blue, levels=0:0.01:0.1, label="V11", colorbar=false) contour!(x, z, v2, color=:red, levels=0:0.01:0.1, label="V22") ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/NQCModels/ase.md b/docs/src/NQCModels/ase.md index 0f4739b96..354101ed5 100644 --- a/docs/src/NQCModels/ase.md +++ b/docs/src/NQCModels/ase.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/NQCModels/ase.md..." +start_time = time() +``` # ASE interface The easiest way to obtain potentials and forces from established codes is to @@ -23,12 +27,12 @@ associated calculator to implement the required [`potential`](@ref) and 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 PyCall +using PythonCall: pyimport, pylist +using ASEconvert -ase = pyimport("ase") emt = pyimport("ase.calculators.emt") -h2 = ase.Atoms("H2", [(0, 0, 0), (0, 0, 0.74)]) +h2 = ase.Atoms("H2", pylist([(0, 0, 0), (0, 0, 0.74)])) h2.calc = emt.EMT() nothing # hide ``` @@ -53,3 +57,7 @@ derivative(model, rand(3, 2)) [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. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/NQCModels/frictionmodels.md b/docs/src/NQCModels/frictionmodels.md index b6450267b..5b73bfd23 100644 --- a/docs/src/NQCModels/frictionmodels.md +++ b/docs/src/NQCModels/frictionmodels.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/NQCModels/frictionmodels.md..." +start_time = time() +``` # [Electronic friction models](@id models-friction) To perform [molecular dynamics with electronic friction (MDEF)](@ref mdef-dynamics) @@ -69,3 +73,7 @@ to obtain the time-dependent perturbation theory friction from the atomic positi As with LDFA, one of these models is used in the [reactive scattering example](@ref example-h2scattering). +```@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 06552d504..e8011c54a 100644 --- a/docs/src/NQCModels/neuralnetworkmodels.md +++ b/docs/src/NQCModels/neuralnetworkmodels.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/NQCModels/neuralnetworkmodels.md..." +start_time = time() +``` # Neural network models Using the [ASE interface](@ref) we can directly use models trained using @@ -59,3 +63,7 @@ r = [0 0; 0 0; 0 ustrip(auconvert(0.74u"Å"))] potential(model, r) derivative(model, r) ``` +```@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 0a7f9cd69..8e4567cc8 100644 --- a/docs/src/NQCModels/overview.md +++ b/docs/src/NQCModels/overview.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/NQCModels/overview.md..." +start_time = time() +``` # NQCModels.jl To perform nonadiabatic molecular dynamics simulations, it is necessary to define @@ -102,3 +106,7 @@ AbstractTrees.print_tree(Model) # hide 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/NQCBase/nqcbase.md b/docs/src/api/NQCBase/nqcbase.md index e8f1dbabb..92dc31056 100644 --- a/docs/src/api/NQCBase/nqcbase.md +++ b/docs/src/api/NQCBase/nqcbase.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCBase/nqcbase.md..." +start_time = time() +``` # NQCBase ```@autodocs Modules=[NQCBase] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDistributions/nqcdistributions.md b/docs/src/api/NQCDistributions/nqcdistributions.md index 78298807a..6ad334e7b 100644 --- a/docs/src/api/NQCDistributions/nqcdistributions.md +++ b/docs/src/api/NQCDistributions/nqcdistributions.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCDistributions/nqcdistributions.md..." +start_time = time() +``` # NQCDistributions ```@autodocs Modules=[NQCDistributions] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/calculators.md b/docs/src/api/NQCDynamics/calculators.md index de5822883..934009651 100644 --- a/docs/src/api/NQCDynamics/calculators.md +++ b/docs/src/api/NQCDynamics/calculators.md @@ -1,5 +1,13 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/calculators.md..." +start_time = time() +``` # Calculators ```@autodocs Modules=[NQCDynamics.Calculators] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/dynamicsmethods.md b/docs/src/api/NQCDynamics/dynamicsmethods.md index 3fca2317d..b1f2e39fc 100644 --- a/docs/src/api/NQCDynamics/dynamicsmethods.md +++ b/docs/src/api/NQCDynamics/dynamicsmethods.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/dynamicsmethods.md..." +start_time = time() +``` # DynamicsMethods @@ -34,3 +38,7 @@ Modules=[NQCDynamics.DynamicsMethods.EhrenfestMethods] ```@autodocs Modules=[NQCDynamics.DynamicsMethods.IntegrationAlgorithms] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/dynamicsoutputs.md b/docs/src/api/NQCDynamics/dynamicsoutputs.md index d998122f9..8abea7684 100644 --- a/docs/src/api/NQCDynamics/dynamicsoutputs.md +++ b/docs/src/api/NQCDynamics/dynamicsoutputs.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/dynamicsoutputs.md..." +start_time = time() +``` # DynamicsOutputs @@ -15,3 +19,7 @@ Private=false Modules=[NQCDynamics.DynamicsOutputs] Public=false ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/dynamicsutils.md b/docs/src/api/NQCDynamics/dynamicsutils.md index 9bb644d68..f23ff3011 100644 --- a/docs/src/api/NQCDynamics/dynamicsutils.md +++ b/docs/src/api/NQCDynamics/dynamicsutils.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/dynamicsutils.md..." +start_time = time() +``` # DynamicsUtils ```@autodocs Modules=[NQCDynamics.DynamicsUtils] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/ensembles.md b/docs/src/api/NQCDynamics/ensembles.md index 4f94345d7..36a1f7e05 100644 --- a/docs/src/api/NQCDynamics/ensembles.md +++ b/docs/src/api/NQCDynamics/ensembles.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/ensembles.md..." +start_time = time() +``` # Ensembles ```@autodocs Modules=[NQCDynamics.Ensembles] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/estimators.md b/docs/src/api/NQCDynamics/estimators.md index fb8c7232a..d7fec7a37 100644 --- a/docs/src/api/NQCDynamics/estimators.md +++ b/docs/src/api/NQCDynamics/estimators.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/estimators.md..." +start_time = time() +``` # Estimators ```@autodocs Modules=[NQCDynamics.Estimators] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/initialconditions.md b/docs/src/api/NQCDynamics/initialconditions.md index bafb46be1..59d86aabf 100644 --- a/docs/src/api/NQCDynamics/initialconditions.md +++ b/docs/src/api/NQCDynamics/initialconditions.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/initialconditions.md..." +start_time = time() +``` # InitialConditions @@ -22,3 +26,7 @@ Modules=[NQCDynamics.InitialConditions.QuantisedDiatomic] ```@autodocs Modules=[NQCDynamics.InitialConditions.MetropolisHastings] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/nonadiabaticmoleculardynamics.md b/docs/src/api/NQCDynamics/nonadiabaticmoleculardynamics.md index 27e6affb5..f9605474f 100644 --- a/docs/src/api/NQCDynamics/nonadiabaticmoleculardynamics.md +++ b/docs/src/api/NQCDynamics/nonadiabaticmoleculardynamics.md @@ -1,5 +1,13 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/nonadiabaticmoleculardynamics.md..." +start_time = time() +``` # NQCDynamics ```@autodocs Modules=[NQCDynamics] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/numericutils.md b/docs/src/api/NQCDynamics/numericutils.md index 1ebbbd03a..e182434e0 100644 --- a/docs/src/api/NQCDynamics/numericutils.md +++ b/docs/src/api/NQCDynamics/numericutils.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/numericutils.md..." +start_time = time() +``` # Numerical utilities ```@autodocs Modules=[NQCDynamics.FastDeterminant] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/ringpolymers.md b/docs/src/api/NQCDynamics/ringpolymers.md index e9ac27ae7..441797c42 100644 --- a/docs/src/api/NQCDynamics/ringpolymers.md +++ b/docs/src/api/NQCDynamics/ringpolymers.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/ringpolymers.md..." +start_time = time() +``` # RingPolymers ```@autodocs Modules=[NQCDynamics.RingPolymers] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCDynamics/timecorrelationfunctions.md b/docs/src/api/NQCDynamics/timecorrelationfunctions.md index 12f030d45..ae72e6678 100644 --- a/docs/src/api/NQCDynamics/timecorrelationfunctions.md +++ b/docs/src/api/NQCDynamics/timecorrelationfunctions.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCDynamics/timecorrelationfunctions.md..." +start_time = time() +``` # TimeCorrelationFunctions ```@autodocs Modules=[NQCDynamics.TimeCorrelationFunctions] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCModels/adiabaticmodels.md b/docs/src/api/NQCModels/adiabaticmodels.md index 3b19ec29f..0e723a2ff 100644 --- a/docs/src/api/NQCModels/adiabaticmodels.md +++ b/docs/src/api/NQCModels/adiabaticmodels.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCModels/adiabaticmodels.md..." +start_time = time() +``` # AdiabaticModels ```@autodocs Modules=[NQCModels.AdiabaticModels] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCModels/cubeldfamodel.md b/docs/src/api/NQCModels/cubeldfamodel.md index 7e7f4ede8..4277b1fad 100644 --- a/docs/src/api/NQCModels/cubeldfamodel.md +++ b/docs/src/api/NQCModels/cubeldfamodel.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCModels/cubeldfamodel.md..." +start_time = time() +``` # CubeLDFAModel ```@autodocs Modules=[CubeLDFAModel] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCModels/diabaticmodels.md b/docs/src/api/NQCModels/diabaticmodels.md index 7603ca18d..b31316696 100644 --- a/docs/src/api/NQCModels/diabaticmodels.md +++ b/docs/src/api/NQCModels/diabaticmodels.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCModels/diabaticmodels.md..." +start_time = time() +``` # DiabaticModels ```@autodocs Modules=[NQCModels.DiabaticModels] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCModels/frictionmodels.md b/docs/src/api/NQCModels/frictionmodels.md index 31349d849..edc8151b7 100644 --- a/docs/src/api/NQCModels/frictionmodels.md +++ b/docs/src/api/NQCModels/frictionmodels.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCModels/frictionmodels.md..." +start_time = time() +``` # FrictionModels ```@autodocs Modules=[NQCModels.FrictionModels] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCModels/nninterfaces.md b/docs/src/api/NQCModels/nninterfaces.md index 02aaec52c..9cb5a6c87 100644 --- a/docs/src/api/NQCModels/nninterfaces.md +++ b/docs/src/api/NQCModels/nninterfaces.md @@ -1,6 +1,14 @@ +```@setup logging +@info "Expanding src/api/NQCModels/nninterfaces.md..." +start_time = time() +``` # NNInterfaces ```@autodocs Modules=[NNInterfaces] ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/api/NQCModels/nonadiabaticmodels.md b/docs/src/api/NQCModels/nonadiabaticmodels.md index 9139d38dc..69698e4a0 100644 --- a/docs/src/api/NQCModels/nonadiabaticmodels.md +++ b/docs/src/api/NQCModels/nonadiabaticmodels.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/api/NQCModels/nonadiabaticmodels.md..." +start_time = time() +``` # NQCModels @@ -7,3 +11,7 @@ Modules=[NQCModels] +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/atoms.md b/docs/src/atoms.md index 1e404f343..6c43c01a7 100644 --- a/docs/src/atoms.md +++ b/docs/src/atoms.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/atoms.md..." +start_time = time() +``` # [Handling Atoms](@id atoms) This package makes the choice to separate the atomic parameters from their positions and @@ -82,3 +86,7 @@ AtomsIO.save_trajectory("Si.xyz", trajectory) # Save a trajectory AtomsIO also provides `load_system` and `load_trajectory` which can be converted to the NQCDynamics format as above to initialise simulations. Refer to [AtomsIO](https://mfherbst.github.io/AtomsIO.jl/stable/) for more information. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/devdocs/diffeq.md b/docs/src/devdocs/diffeq.md index 470e5a7c3..febd73d8a 100644 --- a/docs/src/devdocs/diffeq.md +++ b/docs/src/devdocs/diffeq.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/devdocs/diffeq.md..." +start_time = time() +``` # DifferentialEquations.jl integration @@ -56,3 +60,7 @@ the simulation cell, and the termination caused the simulation to exit early. The callback setup we're using is exactly that provided by DifferentialEquations.jl, if you want more details on callbacks, please refer to their [documentation](https://diffeq.sciml.ai/dev/features/callback_functions/#callbacks). +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/devdocs/models.md b/docs/src/devdocs/models.md index 9df43a58c..cda85c399 100644 --- a/docs/src/devdocs/models.md +++ b/docs/src/devdocs/models.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/devdocs/models.md..." +start_time = time() +``` # [Implementing a new model](@id devdocs-model) [NQCModels.jl](@ref) aims to provide a unified interface for defining model @@ -49,3 +53,7 @@ For further examples, you can also take a look into the source code of the `NQCM package to see how the analytic models have been implemented. If you have any issues or questions about implementing a new model, open up an issue on Github and we can work together to resolve the problem. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/devdocs/new_methods.md b/docs/src/devdocs/new_methods.md index e8432045c..9a35503e6 100644 --- a/docs/src/devdocs/new_methods.md +++ b/docs/src/devdocs/new_methods.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/devdocs/new_methods.md..." +start_time = time() +``` # Contributing a new method @@ -162,3 +166,7 @@ The `DifferentialEquations.jl` framework provides a simple interface for adding algorithms, check out the [developer documentation](https://devdocs.sciml.ai/dev/) to learn how it works. You can also find some examples of custom algorithms in the `DynamicsMethods.IntegrationAlgorithms` module. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/classical.md b/docs/src/dynamicssimulations/dynamicsmethods/classical.md index f583e0737..54d37abf2 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/classical.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/classical.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/classical.md..." +start_time = time() +``` # [Classical molecular dynamics](@id classical-dynamics) Classical (molecular) dynamics proceeds @@ -33,3 +37,7 @@ traj = run_dynamics(sim, (0.0, 10.0), u0; dt=0.1, output=OutputPosition) plot(traj, :OutputPosition) ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/ehrenfest.md b/docs/src/dynamicssimulations/dynamicsmethods/ehrenfest.md index 5dd5033e4..f8a504016 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/ehrenfest.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/ehrenfest.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/ehrenfest.md..." +start_time = time() +``` # [Ehrenfest molecular dynamics](@id ehrenfest-dynamics) The Ehrenfest method is a mixed quantum-classical dynamics method in which the total wavefunction is factorized into slow (nuclear) variables, which are treated classically, and fast ones (electrons) which remain quantum-mechanical. In the Ehrenfest method, nuclei move according to classical mechanics on a potential energy surface given by the expectation value of the electronic Hamiltonian. @@ -45,3 +49,7 @@ using Plots histogram(momenta) xlims!(-20,20) ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/fssh.md b/docs/src/dynamicssimulations/dynamicsmethods/fssh.md index 61a37259e..f71744c34 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/fssh.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/fssh.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/fssh.md..." +start_time = time() +``` # [Fewest-switches surface hopping (FSSH)](@id fssh-dynamics) Tully's FSSH [Tully1990](@cite) is one of the most popular methods for nonadiabatic @@ -120,3 +124,7 @@ plot(traj, :OutputDiabaticPopulation) [Another example is available](@ref examples-tully-model-two) where we use FSSH and other methods to reproduce some of the results from [Tully1990](@cite). +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/langevin.md b/docs/src/dynamicssimulations/dynamicsmethods/langevin.md index ab8d9b11e..58ce6195b 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/langevin.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/langevin.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/langevin.md..." +start_time = time() +``` # [Classical Langevin dynamics](@id langevin-dynamics) Langevin dynamics can be used to sample the canonical ensemble for a classical system. @@ -122,3 +126,7 @@ Estimators.@estimate kinetic_energy(sim, traj[:OutputVelocity]) ``` Again, this takes a similar value since the total energy is evenly split between the kinetic and potential for a classical harmonic system. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/mdef.md b/docs/src/dynamicssimulations/dynamicsmethods/mdef.md index a86e3f1db..9b4719e56 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/mdef.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/mdef.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/mdef.md..." +start_time = time() +``` # [Molecular dynamics with electronic friction (MDEF)](@id mdef-dynamics) ## Introduction @@ -117,3 +121,7 @@ View the [friction models page](@ref models-friction) to learn about how this ca If you would like to see an example using both LDFA and TDPT during full dimensional dynamics, refer to the [reactive scattering example](@ref example-h2scattering). +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/nrpmd.md b/docs/src/dynamicssimulations/dynamicsmethods/nrpmd.md index eb7c33489..bd2c178bb 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/nrpmd.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/nrpmd.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/nrpmd.md..." +start_time = time() +``` # [Nonadiabatic ring polymer molecular dynamics (NRPMD)](@id nrpmd-dynamics) ## Theory @@ -179,3 +183,7 @@ plot(0:0.1:30, [p[1,1]-p[2,1] for p in ensemble[:PopulationCorrelationFunction]] xlabel!("Time") ylabel!("Population difference") ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/rpmd.md b/docs/src/dynamicssimulations/dynamicsmethods/rpmd.md index 0662509b7..b61f24d3e 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/rpmd.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/rpmd.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/rpmd.md..." +start_time = time() +``` # [Ring polymer molecular dynamics (RPMD)](@id rpmd-dynamics) Ring polymer molecular dynamics is a quantum dynamics methods that attempts @@ -111,3 +115,7 @@ end Since this package is focused on nonadiabatic dynamics, you won't see much adiabatic RPMD elsewhere in the documentation, but it's useful to understand how the original adiabatic version works before moving on to the nonadiabatic extensions. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicsmethods/rpsh.md b/docs/src/dynamicssimulations/dynamicsmethods/rpsh.md index c79d1f3c7..991ba53d3 100644 --- a/docs/src/dynamicssimulations/dynamicsmethods/rpsh.md +++ b/docs/src/dynamicssimulations/dynamicsmethods/rpsh.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicsmethods/rpsh.md..." +start_time = time() +``` # [Ring polymer surface hopping (RPSH)](@id rpsh-dynamics) Ring polymer surface hopping was one of the early attempts to extend @@ -90,3 +94,7 @@ and it is possible that better results could be achieved using a free ring polym [Welsch2016](@cite) provides a theoretical description of how nonequilibrium simulations using RPMD should be performed. This techniques here should likely be applied to RPSH too. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/dynamicssimulations/dynamicssimulations.md b/docs/src/dynamicssimulations/dynamicssimulations.md index 836e4a9e3..be5d4d953 100644 --- a/docs/src/dynamicssimulations/dynamicssimulations.md +++ b/docs/src/dynamicssimulations/dynamicssimulations.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/dynamicssimulations/dynamicssimulations.md..." +start_time = time() +``` # [Introduction](@id dynamicssimulations) Performing dynamics simulations is at the core of this package's functionality @@ -92,3 +96,7 @@ plot(out, :OutputAdiabaticPopulation) All of the dynamics methods work in a similar way. For details on a specific method along with examples, please see the method specific page in the sidebar under `Dynamics methods`. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/ensemble_simulations.md b/docs/src/ensemble_simulations.md index 91223db88..7ffc0bd0a 100644 --- a/docs/src/ensemble_simulations.md +++ b/docs/src/ensemble_simulations.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/ensemble_simulations.md..." +start_time = time() +``` # [Ensemble simulations](@id ensembles) Typically we'll be interested in computing observables based upon the statistics @@ -123,3 +127,7 @@ Throughout the documentation, ensemble simulations like this one are used to dem many of the dynamics methods. Now that you have understood the contents of this page, all of the ensemble simulations will appear familiar. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/examples/reactive_scattering.md b/docs/src/examples/reactive_scattering.md index a6ece42b6..d91f48a27 100644 --- a/docs/src/examples/reactive_scattering.md +++ b/docs/src/examples/reactive_scattering.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/examples/reactive_scattering.md..." +start_time = time() +``` # [Reactive scattering from a metal surface](@id example-h2scattering) Our implementation allows us to simulate vibrational de-excitation probability during reactive scattering events at metal surfaces for any diatomic molecule @@ -153,3 +157,7 @@ individual trajectory. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/examples/spinboson.md b/docs/src/examples/spinboson.md index d3db5a3a1..e04902837 100644 --- a/docs/src/examples/spinboson.md +++ b/docs/src/examples/spinboson.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/examples/spinboson.md..." +start_time = time() +``` # Ohmic spin-boson nonequilibrium population dynamics The spin-boson model is widely used as a model for condensed phase quantum dynamics. @@ -77,3 +81,7 @@ We can see that even with just 100 trajectories, our Ehrenfest result closely ma The FSSH is quite clearly underconverged with only 100 trajectories due to the discontinuous nature of the individual trajectories. Feel free to try this for yourself and see what the converged FSSH result looks like! +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/examples/threestatemorse.md b/docs/src/examples/threestatemorse.md index 7e09b1f3f..d13f219cd 100644 --- a/docs/src/examples/threestatemorse.md +++ b/docs/src/examples/threestatemorse.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/examples/threestatemorse.md..." +start_time = time() +``` # Time-dependent populations with the ThreeStateMorse model In this example we investigate the time-dependent populations of the three state @@ -103,3 +107,7 @@ To reduce the build time for the documentation the results here are underconverg already it is clear that both of these methods come close to the exact result shown by [Coronado2001](@cite). After performing enough trajectories to converge the population dynamics, we would be better able to judge the effectiveness of FSSH and Ehrenfest at reproducing the exact quantum dynamics for this model. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/examples/tully_scattering.md b/docs/src/examples/tully_scattering.md index b91b2f57e..0a344f21c 100644 --- a/docs/src/examples/tully_scattering.md +++ b/docs/src/examples/tully_scattering.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/examples/tully_scattering.md..." +start_time = time() +``` # [Scattering probabilities for TullyModelTwo](@id examples-tully-model-two) In this section we aim to reproduce the results of Fig. 5 from [Tully1990](@cite). @@ -108,3 +112,7 @@ state 1 and transmission on state 2 respectively. For this model, surface hopping is successful in closely approximating the exact quantum result, especially at higher momentum values. Refer to [Tully1990](@cite) and [Shakib2017](@cite) for a detailed discussion of the results. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md index b40520b03..8ec6d1768 100644 --- a/docs/src/getting_started.md +++ b/docs/src/getting_started.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/getting_started.md..." +start_time = time() +``` # Getting started To get started with the package we can identify the necessary ingredients to @@ -278,3 +282,7 @@ Now that we've covered the basics of classical dynamics, we're ready to explore world of nonadiabatic dynamics. All the dynamics methods follow these patterns and anything you find elsewhere in the documentation should now seem relatively familiar. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/index.md b/docs/src/index.md index 16cd7844a..5f0d40033 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/index.md..." +start_time = time() +``` # Introduction Welcome to the documentation for NQCDynamics, @@ -115,3 +119,7 @@ The first page to read is the [Getting started](@ref) section which walks throug needed to perform a conventional classical molecular dynamics simulation. After this, the reader is free to explore at their leisure since everything else builds directly upon sections from the [Getting started](@ref) page. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/initialconditions/ebk.md b/docs/src/initialconditions/ebk.md index 757726eb0..b98610f5c 100644 --- a/docs/src/initialconditions/ebk.md +++ b/docs/src/initialconditions/ebk.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/initialconditions/ebk.md..." +start_time = time() +``` # [Semiclassical EBK quantisation](@id ebk-sampling) In surface science, it is often of interest to investigate how collisions with surfaces @@ -94,3 +98,7 @@ procedure, but this technique can use any arbitrary potential. In the [hydrogen scattering example](@ref example-h2scattering) we build on this example and use the sample procedure to perform scattering simulations starting from this distribution. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/initialconditions/hamiltonian.md b/docs/src/initialconditions/hamiltonian.md index e305bb0e6..fbef92bb4 100644 --- a/docs/src/initialconditions/hamiltonian.md +++ b/docs/src/initialconditions/hamiltonian.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/initialconditions/hamiltonian.md..." +start_time = time() +``` # [Thermal Hamiltonian Monte Carlo](@id hmc-sampling) Our implementation of Hamiltonian Monte Carlo (HMC) is a light wrapper around the @@ -39,3 +43,7 @@ with the equipartition theorem: Estimators.@estimate potential_energy(sim, chain) austrip(sim.temperature) * 3 * 4 / 2 ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/initialconditions/langevin.md b/docs/src/initialconditions/langevin.md index 991cb5f55..d080f522b 100644 --- a/docs/src/initialconditions/langevin.md +++ b/docs/src/initialconditions/langevin.md @@ -1,4 +1,12 @@ +```@setup logging +@info "Expanding src/initialconditions/langevin.md..." +start_time = time() +``` # [Thermal Langevin dynamics](@id langevin-sampling) Thermal initial conditions can be obtained directly from a Langevin dynamics simulations. See the [Langevin dynamics page](@ref langevin-dynamics) for more info. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/initialconditions/metropolishastings.md b/docs/src/initialconditions/metropolishastings.md index fdc748a23..e188b4295 100644 --- a/docs/src/initialconditions/metropolishastings.md +++ b/docs/src/initialconditions/metropolishastings.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/initialconditions/metropolishastings.md..." +start_time = time() +``` # [Thermal Metropolis-Hastings Monte Carlo](@id mhmc-sampling) Metropolis-Hastings Monte Carlo is a popular method for sampling the canonical @@ -113,3 +117,7 @@ presented in the previous section and `output.R` can be readily passed to provid the position distribution. The Monte Carlo sampling does not include velocities but these can be readily obtained from the Maxwell-Boltzmann distribution. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/integration_algorithms.md b/docs/src/integration_algorithms.md index 294ebd321..06686a7fe 100644 --- a/docs/src/integration_algorithms.md +++ b/docs/src/integration_algorithms.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/integration_algorithms.md..." +start_time = time() +``` # Integration algorithms At the core of NQCDynamics.jl is the DifferentialEquations.jl package that performs all of the dynamics simulations. @@ -116,3 +120,7 @@ Using this approach, NQCDynamics.jl implements the BAOAB algorithm for tensorial The `DynamicalSDEProblem` was originally implemented for performing Langevin thermostatted dynamics simulations using the BAOAB algorithm. At the time of writing, BAOAB is the only algorithm implemented in StochasticDiffEq.jl for these problems. In future it would be useful to implement further algorithms and allow for more general noise profiles. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/references.md b/docs/src/references.md index 78f29bc41..ca5e6f7c4 100644 --- a/docs/src/references.md +++ b/docs/src/references.md @@ -1,4 +1,12 @@ +```@setup logging +@info "Expanding src/references.md..." +start_time = time() +``` # References ```@bibliography ``` +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/docs/src/saving_loading.md b/docs/src/saving_loading.md index 23f1a560e..f04699859 100644 --- a/docs/src/saving_loading.md +++ b/docs/src/saving_loading.md @@ -1,3 +1,7 @@ +```@setup logging +@info "Expanding src/saving_loading.md..." +start_time = time() +``` # [Saving and loading](@id saving-and-loading) If you would like to split your workflow into multiple scripts @@ -49,3 +53,7 @@ model = load("parameters.jld2", "model") [JLD2](https://github.com/JuliaIO/JLD2.jl) is compatible with any Julia type so it widely usable for most of the types you encounter is NQCDynamics.jl and across all Julia packages. +```@setup logging +runtime = round(time() - start_time; digits=2) +@info "...done after $runtime s." +``` diff --git a/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl b/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl index f4128850a..e3a1575ad 100644 --- a/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl +++ b/src/DynamicsMethods/IntegrationAlgorithms/mdef_baoab.jl @@ -9,6 +9,18 @@ using NQCDynamics: get_temperature StochasticDiffEq.alg_compatible(::DiffEqBase.AbstractSDEProblem,::MDEF_BAOAB) = true +""" + MDEF_BAOAB(): + We define two types of variables structure in the cache: + - `MDEF_BAOABConstantCache` is used for constant variables, which are not changed during the integration. + + - `MDEF_BAOABCache` is used for variables that are changed during the integration. + + Generally, the mutableChahe is faster than the constant cache. + + Mathmatically speaking, they are doing the same things. MutableCache is more well-optimize for the Julia compiler. +""" + struct MDEF_BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEq.StochasticDiffEqConstantCache k::uType half::uEltypeNoUnits @@ -28,6 +40,9 @@ struct MDEF_BAOABCache{uType,rType,vType,uTypeFlat,uEltypeNoUnits,rateNoiseType, c2::Matrix{uEltypeNoUnits} end +""" + Insecting the inputs into a Cache structure. +""" function StochasticDiffEq.alg_cache(::MDEF_BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) k = zero(rate_prototype.x[1]) MDEF_BAOABConstantCache(k, uEltypeNoUnits(1//2)) @@ -52,6 +67,15 @@ function StochasticDiffEq.alg_cache(::MDEF_BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype MDEF_BAOABCache(tmp, utmp, dutmp, k, flatdutmp, tmp1, tmp2, gtmp, noise, half, c1, c2) end +""" + Initialize the the velocities du and positions u. + + k : gradient of the potential energy + + integrator.f.f1 : the gradient operator ---> acceleration!() in diabtic_mdef.jl + +""" + function StochasticDiffEq.initialize!(integrator, cache::MDEF_BAOABConstantCache) @unpack t,uprev,p = integrator du1 = uprev.x[1] @@ -68,6 +92,29 @@ function StochasticDiffEq.initialize!(integrator, cache::MDEF_BAOABCache) integrator.f.f1(cache.k,du1,u1,p,t) end +""" + perform_step! perform a single step of BAOAB algorithm. + + uprev: previous configureation(velocities and positions) + dt : length of time step + t : current time + k : gradient of the potential energy in terms of velocities + p : parameters like temperature, mass, etc. + + The B and A steps are the same as Leimkuhler and Matthews (2013) and the O step is the same as Bussi and Parrinello (2007). + + B and A step: https://academic.oup.com/amrx/article/2013/1/34/166771 + + O step : general part https://academic.oup.com/amrx/article/2013/1/34/166771 + + friction part https://journals.aps.org/pre/abstract/10.1103/PhysRevE.75.056707 + and https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.118.256001 (supplementary material) + + friction tensor Λ: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.118.256001 (supplementary material) + + +""" + @muladd function StochasticDiffEq.perform_step!(integrator,cache::MDEF_BAOABConstantCache,f=integrator.f) @unpack t,dt,sqdt,uprev,p,W = integrator @unpack k, half = cache @@ -81,14 +128,16 @@ end u2 = u1 + half*dt*du2 # O - Λ = integrator.g(u2,p,t+dt*half) + Λ = 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))) + # eigen decomposition of Λ γ, c = LAPACK.syev!('V', 'U', Λ) # symmetric eigen clamp!(γ, 0, Inf) c1 = diagm(exp.(-γ.*dt)) c2 = diagm(sqrt.(1 .- diag(c1).^2)) - noise = σ.*W.dW[:] / sqdt - du3 = c*c1*c'*du2[:] + c*c2*c'*noise + noise = σ.*W.dW[:] / sqdt # W.dW[:] noise rescaled by square root of time step to make it normal distributed + du3 = c*c1*c'*du2[:] + c*c2*c'*noise # O step from Leimkuhler and Matthews (2013) du3 = reshape(du3, size(du2)) # A @@ -101,6 +150,12 @@ end integrator.u = ArrayPartition((du, u)) end +""" + The details of the following perform_step! using mutableChahe can be found in steps.jl in intergrationAlgorithm. + + Before the checking, you are recommonded to read the function with ConstantCache input first. +""" + @muladd function StochasticDiffEq.perform_step!(integrator,cache::MDEF_BAOABCache,f=integrator.f) @unpack t,dt,sqdt,uprev,u,p,W = integrator @unpack utmp, dutmp, k, half, gtmp = cache diff --git a/test/Core/calculators.jl b/test/Core/calculators.jl index badc2436b..875d09f34 100644 --- a/test/Core/calculators.jl +++ b/test/Core/calculators.jl @@ -4,6 +4,10 @@ using NQCDynamics.Calculators using LinearAlgebra: tr, Diagonal, eigvecs, eigvals using RingPolymerArrays: RingPolymerArrays +# For the allocation tests check against both backends as we can't be sure which is in use +const MKL_EIGEN_ALLOCATIONS = 54912 +const OPENBLAS_EIGEN_ALLOCATIONS = 56896 + @testset "General constructors" begin model = NQCModels.DoubleWell() @test Calculators.Calculator(model, 1, Float64) isa Calculators.DiabaticCalculator @@ -265,7 +269,8 @@ end @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0 @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0 - @test @allocated(Calculators.evaluate_eigen!(calc, r)) == 56896 # nonzero due to eigenroutines + eigen_allocations = @allocated(Calculators.evaluate_eigen!(calc, r)) + @test (eigen_allocations == MKL_EIGEN_ALLOCATIONS) || (eigen_allocations == OPENBLAS_EIGEN_ALLOCATIONS) @test @allocated(Calculators.evaluate_adiabatic_derivative!(calc, r)) == 0 @test @allocated(Calculators.evaluate_nonadiabatic_coupling!(calc, r)) == 0 end @@ -287,7 +292,8 @@ end @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0 @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0 - @test @allocated(Calculators.evaluate_eigen!(calc, r)) == 568960 # nonzero due to eigenroutines + eigen_allocations = @allocated(Calculators.evaluate_eigen!(calc, r)) / 10 + @test (eigen_allocations == MKL_EIGEN_ALLOCATIONS) || (eigen_allocations == OPENBLAS_EIGEN_ALLOCATIONS) @test @allocated(Calculators.evaluate_adiabatic_derivative!(calc, r)) == 0 @test @allocated(Calculators.evaluate_nonadiabatic_coupling!(calc, r)) == 0 end diff --git a/test/Dynamics/mdef_baoab.jl b/test/Dynamics/mdef_baoab.jl index 6c411f387..990f306a7 100644 --- a/test/Dynamics/mdef_baoab.jl +++ b/test/Dynamics/mdef_baoab.jl @@ -6,6 +6,7 @@ using DiffEqNoiseProcess using DiffEqDevTools using Random using Logging +# Set up the seed Random.seed!(100) u0 = zeros(2,2) @@ -17,8 +18,7 @@ g(u,p,t) = diagm(ones(length(u))) p = Simulation(Atoms([:H,:H]), NQCModels.Free(2); temperature=100) -ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g) -prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,0.5),p) +prob1 = DynamicalSDEProblem(f1_harmonic, f2_harmonic ,g,v0,u0,(0.0,0.5),p) sol1 = solve(prob1,DynamicsMethods.IntegrationAlgorithms.MDEF_BAOAB();dt=1/10,save_noise=true)