|
| 1 | +# # ACEpotentials.jl + AtomsBase.jl Tutorial |
| 2 | + |
| 3 | +# ## Introduction |
| 4 | +# |
| 5 | +# This tutorial demonstrates how the ACEpotentials.jl package |
| 6 | +# interoperates with the AtomsBase.jl ecosystem. |
| 7 | + |
| 8 | +# ## Installation |
| 9 | +# |
| 10 | +# add and load general packages used in this notebook. |
| 11 | + |
| 12 | +using Pkg |
| 13 | +## Uncomment the next line if installing Julia for the first time |
| 14 | +## Pkg.Registry.add("General") |
| 15 | +## Pkg.Registry.add |
| 16 | +Pkg.add(["ExtXYZ", "Unitful", "Distributed", "AtomsCalculators", |
| 17 | + "Molly", "AtomsCalculatorsUtilities", "AtomsBuilder", |
| 18 | + ]) |
| 19 | + |
| 20 | +## ACEpotentials installation: |
| 21 | +## If ACEpotentials has not been installed yet, uncomment the following lines |
| 22 | +## Add the ACE registry, which stores the ACEpotential package information |
| 23 | +Pkg.Registry.add(RegistrySpec(url="https://github.com/ACEsuit/ACEregistry")) |
| 24 | +Pkg.add(["GeomOpt", ]) |
| 25 | +## Pkg.add(["ACEpotentials",]) |
| 26 | + |
| 27 | +# We can check the status of the installed packages. |
| 28 | + |
| 29 | +using Pkg; Pkg.activate(".") |
| 30 | +Pkg.status() |
| 31 | + |
| 32 | +# Import all the packages that we will be using, create some processes |
| 33 | +# for parallel model fitting |
| 34 | + |
| 35 | +using ExtXYZ, Unitful, AtomsCalculators, Distributed, ACEpotentials, |
| 36 | + AtomsCalculatorsUtilities |
| 37 | +using AtomsCalculatorsUtilities.SitePotentials: cutoff_radius |
| 38 | +addprocs(10, exeflags="--project=$(Base.active_project())") |
| 39 | +@everywhere using ACEpotentials |
| 40 | + |
| 41 | +# ## Fit a potential for Cu |
| 42 | +# |
| 43 | +# The tutorial can be adapted trivially to use datasets for Ni, Cu, Li, Mo, Si, Ge. |
| 44 | +# |
| 45 | +# We generate a smallish model (about 300 basis functions) for Cu, using |
| 46 | +# correlation-order 3 (body-order 4), and default for rcut. Then we estimate |
| 47 | +# the model parameters using the standard BLR solver. |
| 48 | + |
| 49 | +## generate a model for Cu |
| 50 | +sym = :Cu |
| 51 | +model = ace1_model(elements = [sym,], order = 3, totaldegree = [ 20, 16, 12 ]) |
| 52 | +@show length_basis(model) |
| 53 | +@show cutoff_radius(model) |
| 54 | +## estimate parameters |
| 55 | +train, test, _ = ACEpotentials.example_dataset("Zuo20_$sym") |
| 56 | +solver = ACEfit.BLR(; factorization = :svd) |
| 57 | +acefit!(train, model; solver=solver); GC.gc() |
| 58 | +## quickly check test errors => 0.5 meV/atom and 0.014 eV/A are ok |
| 59 | +ACEpotentials.compute_errors(test, model); |
| 60 | + |
| 61 | +# ## Geometry Optimization with GeomOpt |
| 62 | +# |
| 63 | +# ( Note: we should use GeometryOptimization.jl, but this is not yet updated to |
| 64 | +# AtomsBase.jl version 0.4. ) |
| 65 | +# |
| 66 | +# First import some stuff + a little hack to make GeomOpt play nice with |
| 67 | +# the latest AtomsBase. This is a shortcoming of DecoratedParticles.jl |
| 68 | +# and requires some updates to fully implement the AtomsBase interface. |
| 69 | + |
| 70 | +using AtomsBuilder, GeomOpt, AtomsCalculators, AtomsBase |
| 71 | +using AtomsBase: FlexibleSystem, FastSystem |
| 72 | +using AtomsCalculators: potential_energy |
| 73 | +function _flexiblesystem(sys) |
| 74 | + c3ll = cell(sys) |
| 75 | + particles = [ AtomsBase.Atom(species(sys, i), position(sys, i)) |
| 76 | + for i = 1:length(sys) ] |
| 77 | + return FlexibleSystem(particles, c3ll) |
| 78 | +end; |
| 79 | + |
| 80 | +# We generate a cubic Cu unit cell, but our potential might not have the same |
| 81 | +# equilibrium bond distance as the default in AtomsBuilder, so we optimize |
| 82 | +# the unit cell. |
| 83 | + |
| 84 | +ucell = bulk(sym, cubic=true) |
| 85 | +ucell, _ = GeomOpt.minimise(ucell, model; variablecell=true) |
| 86 | + |
| 87 | +# We keep the energy of the equilibrated unit cell to later compute the |
| 88 | +# defect formation energy. |
| 89 | + |
| 90 | +Eperat = potential_energy(ucell, model) / length(ucell) |
| 91 | +@show Eperat; |
| 92 | + |
| 93 | +# Now that we have an equilibrated unit cell we enlarge it, and then delete |
| 94 | +# an atom to generate a vacancy defect. |
| 95 | + |
| 96 | +sys = _flexiblesystem(ucell) * (2,2,2) |
| 97 | +deleteat!(sys, 1) |
| 98 | +sys |
| 99 | + |
| 100 | +# Now we do another geometry optimization to get the equilibrium geometry. |
| 101 | + |
| 102 | +vacancy_equil, result = GeomOpt.minimise(sys, model; variablecell = false) |
| 103 | +@show result.g_residual; |
| 104 | + |
| 105 | +# We get an estimate of the formation energy. Note this is likely a poor |
| 106 | +# estimate since we didn't train the model on vacancy configurations. |
| 107 | + |
| 108 | +E_def = potential_energy(vacancy_equil, model) - length(sys) * Eperat |
| 109 | +@show E_def; |
| 110 | + |
| 111 | +# ## Molecular Dynamics with Molly |
| 112 | +# |
| 113 | +# First import some stuff + a little hack to make GeomOpt play nice with |
| 114 | +# the latest AtomsBase. This is a shortcoming of DecoratedParticles.jl |
| 115 | +# and requires some updates to fully implement the AtomsBase interface. |
| 116 | + |
| 117 | +import Molly |
| 118 | +sys = rattle!(bulk(sym, cubic=true) * (2,2,2), 0.03) |
| 119 | +sys_md = Molly.System(sys; force_units=u"eV/Å", energy_units=u"eV") |
| 120 | +temp = 298.0u"K" |
| 121 | +sys_md = Molly.System( |
| 122 | + sys_md; |
| 123 | + general_inters = (model,), |
| 124 | + velocities = Molly.random_velocities(sys_md, temp), |
| 125 | + loggers=(temp=Molly.TemperatureLogger(100),) ) |
| 126 | +## energy = Molly.PotentialEnergyLogger(100),), ) |
| 127 | +## can't add an energy logger because Molly internal energies are per mol |
| 128 | +simulator = Molly.VelocityVerlet( |
| 129 | + dt = 1.0u"fs", |
| 130 | + coupling = Molly.AndersenThermostat(temp, 1.0u"ps"), ) |
| 131 | + |
| 132 | +Molly.simulate!(sys_md, simulator, 1000) |
| 133 | + |
| 134 | +## the temperature seems to fluctuate a bit, but at least it looks stable? |
| 135 | +@info("Temperature history:", sys_md.loggers.temp.history) |
| 136 | + |
0 commit comments