Metadensity Functional Learning for Classical Fluids: Regularizing with Pair Correlations.

1. # %% [markdown]
2. # # Neural functional theory for inhomogeneous fluids – Tutorial
4. # %% [markdown]
5. # This material is available at https://github.com/sfalmo/NeuralDFT-Tutorial and it supplements the following publication:
6. #
7. # **Why neural functionals suit statistical mechanics**
8. # *Florian Sammüller, Sophie Hermann, and Matthias Schmidt, [J. Phys.: Condens. Matter **36**, 243002](https://doi.org/10.1088/1361-648X/ad326f) (2024); [arXiv:2312.04681](https://arxiv.org/abs/2312.04681).*
9. #
10. # We show how the physics of simple fluids can be described with **many-body simulations** and **classical density functional theory** (DFT) and how the two approaches can be combined with the help of machine learning to a **neural functional theory**.
11. # The connection of the different methods is exemplified for the case of the one-dimensional hard rod fluid with hands-on code examples and exercises in the programming language Julia.
12. # Besides the pedagogical overview given here, a more in-depth account of the presented methodology with applications to the hard sphere and Lennard-Jones fluids in three dimensions can be found in
13. #
14. # **Neural functional theory of inhomogeneous fluids: Fundamentals and applications**
15. # *Florian Sammüller, Sophie Hermann, Daniel de las Heras, and Matthias Schmidt, [Proc. Natl. Acad. Sci. **120**, e2312484120](https://doi.org/10.1073/pnas.2312484120) (2023); [arXiv:2307.04539](https://arxiv.org/abs/2307.04539).*
17. # %% [markdown]
18. # ## Technical remarks
20. # %% [markdown]
21. # Running some of the tasks in this notebook requires a considerable amount of compute power, e.g. for the generation of simulation data and for the training of the neural network.
22. # Keep this in mind if you are using an online service with limited resources.
23. # For convenience, we also provide ready-to-use data sets and models, which can be downloaded and applied in the relevant parts.
24. #
25. # Particularly for the machine learning tasks, we recommend using a GPU, as this speeds up training and inference considerably.
26. # Using a recent Nvidia GPU should work out of the box.
27. # For other GPU vendors, you might need to install the corresponding Julia packages, see https://juliagpu.org/ for further details.
29. # %% [markdown]
30. # ### Google Colab setup
32. # %% [markdown]
33. # 1. Connect to a GPU runtime:
34. # - Click the arrow (▼) next to *Connect* and select *Change runtime type*
35. # - Under *Runtime type*, select "julia"
36. # - Under *Hardware accelerator*, select some option with "GPU" in its name
37. # - Click *Connect*
38. # 2. Run the following cell to copy remaining files into your workspace and to install the required packages:
40. # %%
41. run(`git clone https://github.com/sfalmo/NeuralDFT-Tutorial.git`)
42. run(`cp NeuralDFT-Tutorial/simulation.jl NeuralDFT-Tutorial/dft.jl NeuralDFT-Tutorial/neural.jl .`)
44. using Pkg
45. Pkg.add(["BSON", "CUDA", "DelimitedFiles", "FFTW", "Flux", "IJulia", "Plots", "cuDNN", "Dates", "Downloads", "Printf"])
46. Pkg.instantiate()
48. # %% [markdown]
49. # ## Part 1: Many-body simulations
51. # %% [markdown]
52. # Many-body simulations have long become a standard tool for the investigation of classical fluids.
53. # Conceptually, they offer a rather straightforward way to predict the behavior of a fluid from a microscopic description, i.e. by specifying the interactions of its constituent particles.
54. # However, simulations often come with a significant demand of computational resources.
55. #
56. # In the following, we focus on Monte Carlo methods for the description of equilibrium systems.
57. # Specifically, we choose the grand ensemble where the temperature $T$ and the chemical potential $\mu$ are kept fixed (we also specify the length $L$ of the one-dimensional domain).
58. # The statistical mechanics of such a system is then determined by the equilibrium distribution function $\psi(x^{(N)}) \sim \exp(-\beta (U(x^{(N)}) - \mu N))$ where $\beta = 1 / (k_B T)$ with the Boltzmann constant $k_B$.
59. # The potential energy $U(x^{(N)}) = u(x^{(N)}) + \sum_{i=1}^N V_\mathrm{ext}(x_i)$ of a given microstate $x^{(N)} = (x_1, x_2, \dots, x_N)$ of $N$ particles consists of a contribution due to an external potential $V_\mathrm{ext}(x)$ and of the internal energy $u(x^{(N)})$.
60. # If the particles in the fluid possess pairwise interactions, $u(x^{(N)}) = \sum_{i=1}^N \sum_{j>i}^N \phi(|x_j - x_i|)$, where $\phi(r)$ is the interaction potential for a given distance $r$ of two particles.
61. #
62. # Monte Carlo is based on the generation of microstates according to their known equilibrium distribution.
63. # This is done iteratively by mutating an initial state $A$ into a new state $B$ with a probability such that a given distribution $P$ of states is kept intact.
64. # The mutation happens in two stages: a trial transition selects a new state $B$, and a criterion $\mathrm{acc}(A \rightarrow B)$ determines if the new state $B$ shall be accepted or if the system shall be reset to the previous state $A$.
65. # Specifically, a valid choice for this acceptance probability is the Metropolis criterion $\mathrm{acc}(A \rightarrow B) = \min(1, P(B) / P(A))$.
66. #
67. # We now apply this scheme to the grand ensemble, which yields the standard grand canonical Monte Carlo (GCMC) method.
68. # There are three possible trial transitions: i) a particle is moved to a new position, ii) a particle is inserted, iii) a particle is removed.
69. # The goal distribution of states $P$ is the grand canonical equilibrium distribution $\psi$ (see above).
70. # When considering particle displacements, the Metropolis rate becomes
71. # $$
72. # \mathrm{acc}(x^{(N)} \rightarrow \tilde{x}^{(N)}) = \min\left(1, \exp[-\beta(U(\tilde{x}^{(N)}) - U(x^{(N)}))]\right)
73. # $$
74. # for moving a particle and changing the initial microstate $x^{(N)}$ to the new configuration $\tilde{x}^{(N)}$.
75. #
76. # Let us now illustrate this trial move with code.
77. # For simplification, we provide some utilities for common tasks in [simulation.jl](simulation.jl), e.g. `calc_energy` for calculating the potential energy of one particle, `pbc!` for applying periodic boundary conditions, and definitions of containers (structs) which hold the state of the `System` and the `Histograms`.
79. # %%
80. include("simulation.jl");
82. function trial_move(system::System; Δxmax=0.1)
83. if isempty(system.particles) # If no particles are in the system, do nothing
84. return
85. end
86. i = rand(1:length(system.particles)) # Select random particle with index i
87. xbefore = system.particles[i] # Save its initial position
88. Ebefore = calc_energy(system, i) # Calculate the initial potential energy of particle i
89. system.particles[i] += Δxmax * (2 * rand() - 1) # Move particle to a new position
90. pbc!(system, i) # Apply periodic boundary conditions (this places the particle back in the box if it has moved outside of the valid range)
91. Eafter = calc_energy(system, i) # Calculate the potential energy of particle i after it has been moved
92. if rand() > exp(-system.β * (Eafter - Ebefore))
93. system.particles[i] = xbefore # Trial move rejected. Reset particle to previous state
94. end
95. end
97. # %% [markdown]
98. # As the particle number can fluctuate in the grand ensemble, there are additional transitions which add and remove particles at random throughout the simulation.
99. # The acceptance probabilities of these transitions can be derived by considering particle exchanges of the system with a virtual reservoir.
100. # We spare this derivation here and only give the results
101. # $$
102. # \mathrm{acc}(x^{(N)} \rightarrow x^{(N+1)}) = \min\left(1, \frac{L}{N+1} \exp[\beta(\mu - U(x^{(N+1)}) + U(x^{(N)}))]\right),
103. # $$
104. # $$
105. # \mathrm{acc}(x^{(N)} \rightarrow x^{(N-1)}) = \min\left(1, \frac{N}{L} \exp[-\beta(\mu + U(x^{(N-1)}) - U(x^{(N)}))]\right).
106. # $$
107. #
108. # These trial transitions are already provided as `trial_insert` and `trial_delete` in [simulation.jl](simulation.jl), feel free to take a look for the implementation details.
109. #
110. # With all transitions being specified, we can write a simulation loop which consists of an equilibration stage and a stage in which measurements take place.
111. # For each simulation step, we perform a sweep over a fixed number of trial transitions which are chosen at random without introducing bias in the insertions and deletions.
113. # %%
114. using Dates
116. function sweep(system::System; transitions=10, insert_delete_probability=0.2)
117. for _ in 1:transitions
118. if rand() < insert_delete_probability # Randomly select trial transition: either move or insert/delete
119. rand() < 0.5 ? trial_insert(system) : trial_delete(system) # Do trial insertions and removals with equal probability
120. else
121. trial_move(system)
122. end
123. end
124. end
126. function simulate(L::Number, μ::Number, T::Number, Vext::Function, ϕ::Function; equilibration_time=Dates.Second(1), production_time=Dates.Second(2), sweep_transitions=10)
127. system = System(L, μ, T, Vext, ϕ) # The state of the system is encapsulated in this struct
128. histograms = Histograms(system)
129. equilibration_start = now()
130. while now() - equilibration_start < equilibration_time # Equilibration stage, no sampling
131. sweep(system; transitions=sweep_transitions)
132. end
133. production_start = now()
134. while now() - production_start < production_time
135. sweep(system; transitions=sweep_transitions)
136. sample(system, histograms)
137. end
138. get_results(system, histograms) # Normalizes histograms and returns (xs, ρ)
139. end
141. # %% [markdown]
142. # Equilibrium averages such as the one-body density profile $\rho(x) = \langle \sum_{i=1}^N \delta(x - x_i) \rangle$ are obtained by sampling.
143. # For this, the particle configuration is recorded in a position-resolved histogram, which yields the desired average after normalization:
145. # %%
146. function sample(system::System, histograms::Histograms)
147. for x in system.particles
148. bin = ceil(Int, x / L * histograms.bins) # Calculate the bin index from the given particle position
149. histograms.ρ[bin] += 1
150. end
151. histograms.count += 1 # Needed later for normalization
152. end
154. # %% [markdown]
155. # Let us now do some simulations to illustrate the usage of the code.
156. # We just pass the length $L$ of our system, the thermodynamic statepoint $\mu$ and $T$, the external potential $V_\mathrm{ext}(x)$ and the pair interaction potential $\phi(r)$ to the `simulate` function.
157. # In this tutorial, we focus on the hard rod fluid, which is characterized by forbidding the overlap of particles as illustrated below.
158. #
159. # 
160. #
161. # This behavior is reflected by the pair interaction potential
162. # $$
163. # \phi(r) = \begin{cases}
164. # \infty, &r < \sigma,\\
165. # 0, &r \geq \sigma,
166. # \end{cases}
167. # $$
168. # where $\sigma$ specifies the particle size (we set $\sigma = 1$ in the code).
169. # As a simple test case, we choose confinement between hard walls by setting $V_\mathrm{ext}(x) = \infty$ at the boundaries of the domain.
170. # Change the form of the external potential as well as the rest of the system parameters to explore the behavior of the hard rod fluid.
172. # %%
173. using Plots
175. L = 10.0
176. μ, T = 2.0, 1.0
177. Vext(x) = x < 0.5 || x > L - 0.5 ? Inf : 0 # Hard walls at the boundaries
178. ϕ(r) = r < 1.0 ? Inf : 0 # Hard core repulsion
180. xs, ρ = simulate(L, μ, T, Vext, ϕ; equilibration_time=Dates.Second(1), production_time=Dates.Second(2), sweep_transitions=10)
182. plot(xs, ρ, label="ρ", xlims=(0, L), xlabel="x/σ")
184. # %% [markdown]
185. # ## Part 2: Classical density functional theory
187. # %% [markdown]
188. # Classical DFT is founded on a minimization principle of the grand potential $\Omega[\rho]$, which can be expressed as a functional of the one-body density profile $\rho(x)$.
189. # Its strength is therefore the reduction of a many-body problem (see Part 1) to a formally exact description on the level of one-body quantities.
190. # By writing out entropic, external and internal contributions of the grand potential, carrying out the functional derivative $\delta \Omega[\rho] / \delta \rho(x)$ and demanding that it vanishes at the equilibrium density, one arrives at the Euler-Lagrange equation of DFT,
191. # $$
192. # c_1(x) = \ln\rho(x) + \beta (V_\mathrm{ext}(x) - \mu),
193. # $$
194. # where $c_1(x)$ is the one-body direct correlation function.
195. # This rather abstract object attains a fundamental meaning in classical DFT as it captures the nontrivial effects of the internal interactions within the fluid on the one-body level.
196. # Formally, it arises from a functional derivative,
197. # $$
198. # c_1(x; [\rho]) = - \frac{\delta \beta F_\mathrm{exc}[\rho]}{\delta \rho(x)},
199. # $$
200. # where $F_\mathrm{exc}[\rho]$ is the excess free energy, i.e. the part of $\Omega[\rho]$ which accounts for the beyond-ideal-gas behavior of the interacting fluid.
201. # As an immediate consequence of $F_\mathrm{exc}[\rho]$ being a density functional, $c_1(x; [\rho])$ also attains a functional dependence on $\rho(x)$, which is made explicit by the bracket notation.
202. # In principle, $F_\mathrm{exc}[\rho]$ and $c_1(x; [\rho])$ also depend parametrically on the temperature $T$, but this dependence becomes trivial for the hard rod fluid.
203. #
204. # By rearranging the above equation to
205. # $$
206. # \rho(x) = \exp(-\beta(V_\mathrm{ext}(x) - \mu) + c_1(x; [\rho])),
207. # $$
208. # we can reveal its use in actual DFT applications.
209. # Given a suitable expression for $c_1(x; [\rho])$ (which one has to obtain somehow for a given type of model fluid), a self-consistent iteration scheme can be used to solve for $\rho(x)$.
210. # One such method is the Picard iteration with mixing, in which the iteration is performed as follows:
211. # $$
212. # \rho(x) \leftarrow (1 - \alpha) \rho(x) + \alpha \rho_\mathrm{EL}(x).
213. # $$
214. # Here, $\alpha$ is a mixing parameter and $\rho_\mathrm{EL}(x) = \exp(-\beta(V_\mathrm{ext}(x) - \mu) + c_1(x; [\rho]))$ is the right hand side of the rearranged Euler-Lagrange equation.
215. #
216. # A simple DFT program proceeds as follows:
217. # The system parameters $\mu$ and $T$, an external potential $V_\mathrm{ext}(x)$ and a functional form of $c_1(x; [\rho])$ are given, and the density profile is initialized (e.g. with a constant value) on a discrete numerical grid that spans the domain of length $L$.
218. # The iteration is then started and Picard steps are performed to update $\rho(x)$ as described above.
219. # If the change of $\rho(x)$ between iteration steps falls below a predefined numerical tolerance, the iteration is stopped and the converged self-consistent density profile is obtained as a result.
220. # We give an example of such a program in the following:
222. # %%
223. function minimize(L::Number, μ::Number, T::Number, Vext::Function, get_c1::Function; α::Number=0.03, maxiter::Int=10000, dx::Number=0.01, floattype::Type=Float32, tol::Number=max(eps(floattype(1e3)), 1e-8))
224. L, μ, T = floattype.((L, μ, T)) # Technical detail: we will use Float32 in the machine learning part as neural networks usually operate on single precision floats
225. xs = collect(floattype, dx/2:dx:L) # Construct the numerical grid
226. Vext = Vext.(xs) # Evaluate the external potential on the grid
227. infiniteVext = isinf.(Vext) # Check where Vext is infinite to set ρ = 0 there
228. ρ, ρEL = zero(xs), zero(xs) # Preallocate the density profile and an intermediate buffer for iteration
229. fill!(ρ, 0.5) # Start with a bulk density of 0.5
230. c1 = get_c1(xs) # Obtain the c1 functional for the given numerical grid
231. i = 0
232. while true
233. ρEL .= exp.((μ .- Vext) ./ T .+ c1(ρ)) # Evaluate the RHS of the Euler-Lagrange equation
234. ρ .= (1 - α) .* ρ .+ α .* ρEL # Do a Picard iteration step to update ρ
235. ρ[infiniteVext] .= 0 # Set ρ to 0 where Vext = ∞
236. clamp!(ρ, 0, Inf) # Make sure that ρ does not become negative
237. Δρmax = maximum(abs.(ρ - ρEL)[.!infiniteVext]) # Calculate the remaining discrepancy to check convergence
238. i += 1
239. if Δρmax < tol
240. println("Converged (step: $(i), ‖Δρ‖ = $(Δρmax) < $(tol) = tolerance)")
241. break # The remaining discrepancy is below the tolerance: break out of the loop and return the result
242. end
243. if !isfinite(Δρmax) || i >= maxiter
244. println("Did not converge (step: $(i) of $(maxiter), ‖Δρ‖: $(Δρmax), tolerance: $(tol))")
245. return nothing # The iteration did not converge, there is no valid result
246. end
247. end
248. xs, ρ
249. end
251. # %% [markdown]
252. # The way in which $c_1(x; [\rho])$ appears both in the theory as well as in the code example above seems innocuous at first.
253. # However, the crux of DFT is finding a suitable functional expression for $c_1(x; [\rho])$ (or equivalently for $F_\mathrm{exc}[\rho]$) for a given fluid model, and much of the ongoing research deals exactly with this problem.
254. # In Part 3, we will show how to efficiently use neural networks to capture such a nontrivial functional mapping from simulation data.
255. #
256. # In the following, we proceed analytically as the focus lies on the hard rod system, which is at present one of very few model fluids where the exact excess free energy functional could be found.
257. # This success can be traced back to the purely geometrical nature of this particular system.
258. # For hard spheres in 3D, fundamental measure theory provides a similar geometric approach to the construction of excess functionals, although the results remain approximate in this more general geometry.
259. #
260. # For the exact hard rod density functional, which is attributed to Jerry Percus, the excess free energy can be expressed as
261. # $$
262. # \beta F_\mathrm{exc}[\rho] = \int \mathrm{d}x \Phi(n_0(x; [\rho]), n_1(x; [\rho]))
263. # $$
264. # where the free energy density has the specific form $\Phi(n_0, n_1) = - n_0 \ln(1 - n_1)$.
265. # The functions $n_0(x; [\rho])$ and $n_1(x; [\rho])$ are *weighted* densities which arise from convolutions of the density profile with the weight functions $\omega_0(x) = (\delta(x-R) + \delta(x+R)) / 2$ and $\omega_1(x) = \Theta(R - |x|)$, i.e. $n_\alpha(x; [\rho]) = (\omega_\alpha \star \rho)(x)$.
266. # From this excess free energy functional, one can easily obtain by functional differentiation the result
267. # $$
268. # c_1(x; [\rho]) = - \sum_{\alpha=0,1} \left(\omega_\alpha \star \frac{\partial \Phi}{\partial n_\alpha}\right)(x)
269. # $$
270. # for the one-body direct correlation function.
271. #
272. # We implement a method that constructs the Percus $c_1(x; [\rho])$ for a given numerical grid in the following.
273. # As assistance, we provide `conv_fft` to evaluate convolutions efficiently in Fourier space and `get_weights_Percus` to obtain $\omega_\alpha(x)$ on the numerical grid in the file [dft.jl](dft.jl).
275. # %%
276. include("dft.jl")
278. function get_c1_Percus(xs)
279. ω0, ω1 = get_weights_Percus(xs)
280. conv(f, g) = conv_fft(f, g; dx=xs[2]-xs[1])
281. function (ρ)
282. n0, n1 = conv(ω0, ρ), conv(ω1, ρ) # Calculate the weighted densities
283. ∂ϕ_∂n0 = -log.(1 .- n1) # Evaluate the partial derivatives of the excess free energy density
284. ∂ϕ_∂n1 = n0 ./ (1 .- n1)
285. -(conv(ω0, ∂ϕ_∂n0) .+ conv(ω1, ∂ϕ_∂n1)) # Return c1
286. end
287. end
289. # %% [markdown]
290. # Now we can perform some DFT minimizations to see if the results match the simulations.
291. # Note that in general, you might have to adjust the mixing parameter $\alpha$ to ensure convergence of the Picard iteration (see the keyword arguments of `minimize`).
293. # %%
294. L = 10.0
295. μ, T = 2.0, 1.0
296. Vext(x) = x < 0.5 || x > L - 0.5 ? Inf : 0
298. xs, ρ = minimize(L, μ, T, Vext, get_c1_Percus)
300. plot(xs, ρ, label="ρ", xlims=(0, L), xlabel="x/σ")
302. # %% [markdown]
303. # By comparing simulation and DFT results for the same choices of $L$, $\mu$, $T$ and $V_\mathrm{ext}(x)$, we can cross-check the validity of the density profiles of the hard rod fluid.
304. # However, it is obvious that both methods have substantial (and quite contrary) restrictions:
305. # - The simulation data is noisy, which can only be improved with longer simulation runs. This quickly becomes prohibitively expensive, in particular if one wishes to perform many individual simulations to explore the behavior for different system parameters. However, one is free to change the type of considered fluid by simply modifying the form of the internal interactions.
306. # - The DFT calculation is fast and does not suffer from noisy results, thus enabling vast and efficient parameter studies. The results are exact for the case of hard rods, but as illustrated above, one had to find and implement a suitable density functional in order to capture the nontrivial intrinsic correlations. For more complex fluids, there is little hope in deriving an exact functional analytically, and the construction of good approximations often proves to be very difficult.
307. #
308. # In the final part of this tutorial, we will combine the advantages of both methods with the help of machine learning.
309. # We will demonstrate how to acquire an accurate and flexible representation of $c_1(x; [\rho])$ by training a neural network with simulation data.
310. # The resulting *neural functional* can be used in the minimization scheme as seen above, but it can also reveal more fundamental information about the statistical mechanics of the considered fluid.
312. # %% [markdown]
313. # ## Part 3: Neural functional theory
315. # %% [markdown]
316. # 
317. #
318. # The survey in Part 2 indicates that classical DFT is a powerful concept for the determination of fluid equilibra, but that it is limited by the difficulty of finding analytic expressions for $c_1(x; [\rho])$.
319. # In the following, we will train a neural network with simulation data (see Part 1) to obtain a representation of $c_1(x; [\rho])$.
320. # We will then show applications of this neural functional, which include the prediction of density profiles and the investigation of neural functional calculus.
321. #
322. # To generate a suitable data set, simulations of systems with randomized inhomogeneous external potentials $V_\mathrm{ext}(x)$ and random values of the chemical potential $\mu$ are performed.
323. # For the construction of an appropriate form of $V_\mathrm{ext}(x)$, we combine Fourier modes, linear segments and hard walls, which are defined in the following.
324. # We also write a method which generates a combination of these elementary functions with randomly chosen parameters.
326. # %%
327. Vext_sin(x; n::Int, A::Number, φ::Number, L::Number) = A * sin(2π * x * n / L + φ)
329. Vext_lin(x; x1::Number, x2::Number, E1::Number, E2::Number) = x > x1 && x < x2 ? E1 + (x - x1) * (E2 - E1) / (x2 - x1) : 0
331. Vext_wall(x; xw::Number, L::Number) = x < xw || x > L - xw ? Inf : 0
333. function generate_Vext(L::Number; num_sin=4, num_lin=rand(1:5), wall=true)
334. Avar = 1.0
335. sin_parameters = []
336. for n in 1:num_sin # Generate random parameters for periodic sine functions with increasing frequency
337. push!(sin_parameters, (n = n, A = randn() * Avar, φ = rand() * 2π, L = L))
338. end
339. Evar = 1.0
340. lin_parameters = []
341. for _ in 1:num_lin # Generate random parameters for discontinuous linear segments
342. push!(lin_parameters, (x1 = round(rand() * L, digits=2), x2 = round(rand() * L, digits=2), E1 = randn() * Evar, E2 = randn() * Evar))
343. end
344. xwmax = 1.0
345. wall_params = (xw = round(rand() * xwmax, digits=2), L = L) # Set a random wall width
346. function (x) # Return a method which evaluates a combination of all functions with the chosen parameters above
347. result = 0.0
348. for sin_params in sin_parameters
349. result += Vext_sin(x; sin_params...)
350. end
351. for lin_params in lin_parameters
352. result += Vext_lin(x; lin_params...)
353. end
354. if wall
355. result += Vext_wall(x; wall_params...)
356. end
357. result
358. end
359. end
361. # %% [markdown]
362. # To get an idea of the randomized inhomogeneous environments, let us generate and plot some $V_\mathrm{ext}(x)$ profiles.
363. # Run the following code cell a few times:
365. # %%
366. L = 10.0
367. dx = 0.01
368. xs = dx/2:dx:L
369. Vext_generated = generate_Vext(L)
370. plot(xs, Vext_generated.(xs), label="Vext", xlims=(0, L), xlabel="x/σ")
372. # %% [markdown]
373. # We can now proceed to generate reference data by running simulations with these randomized external potentials and random values of the chemical potential as input.
374. # Of course, this requires spending some computational resources.
375. # You can run the simulations yourself, but be aware that this does take a considerable amount of compute time (~hours) if you want to get good data.
376. # For this, change `use_prepared_simulations` to `false` in the following code cell and set an appropriate number of simulations and the time to spend for the equilibration and production stage.
377. # Alternatively, you could "fake" the simulations by using the exact Percus DFT as illustrated in Part 2.
378. #
379. # For readers without the required computational resources and/or patience, we have prepared a ready-to-use data set for the following machine learning tasks, which is downloaded by default as follows:
381. # %%
382. using DelimitedFiles, Downloads
384. use_prepared_simulations = true # Choose if you want to use pregenerated simulations (true) or if you want to generate a simulation data set from scratch (false)
386. if use_prepared_simulations
387. println("Downloading pregenerated simulation data set")
389. Downloads.download("https://www.staff.uni-bayreuth.de/~bt306964/neuraldft-tutorial/data.tar", "data.tar")
390. run(`tar xf data.tar`) # You should now have a directory "data" with a bunch of simulation files
391. datadir = "data"
392. else
393. num_sim = 512
394. equilibration_time = Dates.Second(10)
395. production_time = Dates.Second(200)
396. nthreads = Threads.nthreads()
397. println("Generating simulation reference data from scratch. This will take approximately $(canonicalize(num_sim * (equilibration_time + production_time) / nthreads)).")
399. L = 10
400. ϕ(r) = r < 1.0 ? Inf : 0 # Hard core repulsion
401. μmin, μmax = -7.0, 5.0
402. T = 1.0
404. datadir = mkdir("data_$(now())")
405. println("Saving results to $(datadir)")
406. Threads.@threads for i in 1:num_sim
407. μ = μmin + rand() * (μmax - μmin)
408. Vext_generated = generate_Vext(L)
409. println("Simulation $(i) running...")
410. xs, ρ = simulate(L, μ, T, Vext_generated, ϕ; equilibration_time, production_time)
411. βμloc = (μ .- Vext_generated.(xs)) / T
412. writedlm("$(datadir)/$(i).dat", [xs βμloc ρ])
413. println("Simulation $(i) done")
414. end
415. end
417. println("There are $(length(readdir(datadir))) result files in the directory $(datadir)")
419. # %% [markdown]
420. # If you have generated your own training set, the results are now located in the directory `data_<timestamp>` (use the file explorer of JupyterLab on the left side).
421. # Otherwise, proceed with the pregenerated simulations which are located in `data` after successful download and unpacking.
422. #
423. # Let us first plot some of the data files to check if they are reasonable:
425. # %%
426. datadir = "data" # ... or use your own simulation results from above
428. sim = rand(readdir(datadir, join=true)) # Select a random simulation within the data directory
429. xs, βμloc, ρ = eachcol(readdlm(sim)) # Read the result file and split columns
430. display(plot(xs, βμloc, label="β μloc", xlims=(0, L), xlabel="x/σ"))
431. display(plot(xs, ρ, label="ρ", xlims=(0, L), xlabel="x/σ"))
433. # %% [markdown]
434. # Note that besides $\rho(x)$, we have saved the local chemical potential $\mu_\mathrm{loc}(x) = \mu - V_\mathrm{ext}(x)$ as a further one-body field.
435. # This additional information suffices to calculate the one-body direct correlation function $c_1(x) = \ln \rho(x) - \beta \mu_\mathrm{loc}(x)$ for each $x$ where $\rho(x) \neq 0$, which is the target of our following investigation.
436. #
437. # In Part 2, we have shown that the one-body direct correlation function cannot only be obtained pointwise by the above relation, but that it also constitutes a universal functional $c_1(x; [\rho])$ of the density profile, without having to invoke the local chemical potential $\mu_\mathrm{loc}(x)$ explicitly.
438. # In the following, we attempt to machine-learn this functional relationship with a neural network.
439. # For constructing the specific input-output mapping, we appeal to some physical background.
440. # Formally, the functional dependence on $\rho(x)$ is given with respect to the profile of the whole system.
441. # However, for short-ranged pair interactions, the influence of the surrounding density profile on the value of $c_1(x; [\rho])$ at a certain position $x$ also remains very short-ranged.
442. # Therefore, we can adopt a *local* learning strategy: the density profile is taken only within a narrow window around a given position as input to the neural network, which is trained to yield the *scalar value* of $c_1(x; [\rho])$ at that position.
443. # The whole one-body direct correlation profile can be recovered by evaluating the neural network at different positions across the system domain.
444. # The following schematic illustrates this setup.
445. #
446. # 
447. #
448. # Due to this specific construction, the data has to be prepared accordingly prior to the training, i.e. the input-output pairs of $\rho$-windows and $c_1$-values have to be generated from the full profiles.
449. # For this, some utility functions are given in the file [neural.jl](neural.jl).
450. # Note that we choose a window width of $1 \sigma$ to each side.
451. # For the hard rod fluid, we know that this choice suffices, as the exact analytic solution for $c_1(x; [\rho])$ also has a convolutional range of $1 \sigma$.
452. # For other fluid models, one would have to determine the window width heuristically, for which the neural functional calculus (see below) could provide some guidance.
454. # %%
455. include("neural.jl")
457. ρ_profiles, c1_profiles = read_sim_data(datadir)
458. ρ_windows, c1_values = generate_inout(ρ_profiles, c1_profiles; window_width=1.0, dx=xs[2]-xs[1]) # window width = 1σ suffices for hard rods
460. size(ρ_windows), size(c1_values)
462. # %% [markdown]
463. # The data is now ready to be given to a neural network for training.
464. # We choose a simple architecture with fully-connected hidden layers and continuous activation functions (e.g. `softplus`).
465. # For the construction of the neural network and the subsequent machine learning routines, the framework [Flux.jl](https://fluxml.ai/) is used.
466. # We also specify a standard optimizer and give the generated input-output pairs to a data loader, which automates shuffling and batching during training.
467. # As this is a common regression task, we select the mean squared error as the loss function and the mean absolute error as the metric.
468. # The learning rate is decreased exponentially after the first few epochs, which improves the final training result.
469. #
470. # You can choose in the following code cell via `use_pretrained_model` whether you want to download a ready-to-use model or whether you want to do the training yourself.
471. # Running the training on a GPU is recommended.
472. # Otherwise, the code runs as-is on the CPU, but training might be slow.
474. # %%
475. using BSON, Dates, Downloads, Flux, Printf
476. using CUDA # For Nvidia GPU support, defaults to CPU if no CUDA device is available. Switch out this package if you have a different GPU manufacturer.
478. use_pretrained_model = true # Choose if you want to use the pretrained model (true) or if you want to do the training yourself (false)
480. if use_pretrained_model
481. println("Downloading pretrained model")
483. Downloads.download("https://www.staff.uni-bayreuth.de/~bt306964/neuraldft-tutorial/model.bson", "model.bson")
484. BSON.@load "model.bson" model
485. else
486. println("Training model from scratch")
487. ρ_windows, c1_values = (ρ_windows, c1_values) |> gpu
489. model = Chain(
490. Dense(size(ρ_windows)[1] => 128, softplus),
491. Dense(128 => 64, softplus),
492. Dense(64 => 32, softplus),
493. Dense(32 => 1)
494. ) |> gpu
496. display(model) # Show a summary of the model with the number of fittable parameters
498. opt = Flux.setup(Adam(), model) # Set up a standard Adam optimizer
500. loader = Flux.DataLoader((ρ_windows, c1_values), batchsize=256, shuffle=true) # Initialize the DataLoader to yield shuffled batches
502. loss(m, x, y) = Flux.mse(m(x), y) # Use mean squared error as loss
503. metric(m, x, y) = Flux.mae(m(x), y) # Use mean absolute error as metric
505. get_learning_rate(epoch; initial=0.0001, rate=0.03, wait=5) = epoch < wait ? initial : initial * (1 - rate)^(epoch - wait)
507. model_savefile = "model_$(now()).bson"
508. println("Saving model to $(model_savefile)")
509. for epoch in 1:250 # Do the training loop
510. learning_rate = get_learning_rate(epoch)
511. Flux.adjust!(opt, learning_rate)
512. @printf "Epoch: %3i (learning_rate: %.2e)..." epoch learning_rate; flush(stdout)
513. Flux.train!(loss, model, loader, opt)
514. @printf " loss: %.5f, metric: %.5f\n" loss(model, ρ_windows, c1_values) metric(model, ρ_windows, c1_values); flush(stdout)
515. model = model |> cpu
516. BSON.@save model_savefile model
517. model = model |> gpu
518. end
519. end
521. model
523. # %% [markdown]
524. # The model is now ready to be used for inference.
525. # For convenience, we write a method that lets us evaluate the whole $c_1$-profile via the trained model for a given density profile.
526. # This is necessary due to the local nature of the mapping, which prevents $\rho(x)$ from being used directly as input.
527. # Instead, as shown above, the density profile must be restructured into appropriate windows which are passed to the model.
529. # %%
530. function get_c1_neural(model, xs)
531. window_bins = length(model.layers[1].weight[1,:]) # Get the number of input bins from the shape of the first layer
532. model = model |> gpu
533. function (ρ)
534. ρ_windows = generate_windows(ρ; window_bins) |> gpu # The helper function generate_windows is defined in neural.jl
535. model(ρ_windows) |> cpu |> vec # Evaluate the model, make sure the result gets back to the CPU, and transpose it to a vector
536. end
537. end
539. # %% [markdown]
540. # The central use of the neural functional is revealed by applying it for the self-consistent calculation of density profiles.
541. # Just like we have used the analytic functional due to Percus in Part 2, we can now employ the neural functional to obtain results for $c_1(x; [\rho])$ within a DFT minimization.
542. # Importantly, the specific minimization scheme to update and iterate $\rho(x)$ until convergence remains unchanged upon switching out the analytic one-body direct correlation functional with its neural counterpart.
543. # Hence, this *neural* DFT is an effective means to circumvent the cumbersome search for analytic functionals while retaining the merits of the DFT minimization procedure.
544. #
545. # We show in the following that neural DFT indeed yields accurate results for the considered hard rod fluid.
546. # Exact reference data for comparison can be obtained immediately by the Percus theory (one would have to resort to simulation results for other fluid models).
548. # %%
549. L = 10.0
550. μ, T = 2.0, 1.0
551. Vext(x) = x < 0.5 || x > L - 0.5 ? Inf : 0
553. xs, ρ_neural = minimize(L, μ, T, Vext, xs -> get_c1_neural(model, xs))
554. xs, ρ_Percus = minimize(L, μ, T, Vext, get_c1_Percus)
556. plot(xs, [ρ_neural ρ_Percus], label=["ρ (neural)" "ρ (Percus)"], linestyle=[:solid :dash], xlims=(0, L), xlabel="x/σ")
558. # %% [markdown]
559. # Although the box size has been kept constant during training and in the above example, we can immediately apply the neural functional for different choices of $L$ due to its local nature.
560. # This also gives us the possibility to run some rather demanding quality checks.
561. # For narrow confinement, which has clearly not been included in our training data, the hard rod fluid yields very interesting density profiles.
562. # The following test reveals that the neural functional is able to extrapolate to these unseen cases and that the results from neural DFT match the analytic theory with surprising accuracy.
564. # %%
565. L = 3.0 # Try other values between 2.0 and 4.0, the density profiles will be interesting
566. μ, T = 2.0, 1.0
567. Vext(x) = x < 0.5 || x > L - 0.5 ? Inf : 0
569. xs, ρ_neural = minimize(L, μ, T, Vext, xs -> get_c1_neural(model, xs))
570. xs, ρ_Percus = minimize(L, μ, T, Vext, get_c1_Percus)
572. plot(xs, [ρ_neural ρ_Percus], label=["ρ (neural)" "ρ (Percus)"], linestyle=[:solid :dash], xlims=(0, L), xlabel="x/σ")
574. # %% [markdown]
575. # We can also increase $L$ significantly compared to the previous cases.
576. # The ability of the neural functional to be used straightforwardly in this manner creates the opportunity of efficient multiscale investigations.
577. # In the following, this is illustrated by considering the sedimentation of the hard rod fluid in a column with hard walls at the top and at the bottom.
578. # The external potential is chosen to slowly increase linearly with height $x$ in order to model the effect of gravity.
579. # Within the sedimentation column, the density closely follows the equation of state of the hard rod fluid.
581. # %%
582. L = 100.0 # You can try to increase L even further
583. μ, T = 2.0, 1.0
584. Vext(x) = x < 0.5 || x > L - 0.5 ? Inf : 0.05 * x # Hard walls at the boundaries and linearly increasing within the sedimentation column
586. xs, ρ_neural = minimize(L, μ, T, Vext, xs -> get_c1_neural(model, xs))
587. xs, ρ_Percus = minimize(L, μ, T, Vext, get_c1_Percus)
589. plot(xs, [ρ_neural ρ_Percus], label=["ρ (neural)" "ρ (Percus)"], linestyle=[:solid :dash], xlims=(0, L), xlabel="x/σ")
591. # %% [markdown]
592. # Besides using the trained model directly in DFT minimizations, we show in the following that one can extract much more fundamental physical information from the neural functional.
593. # For this, we aim at implementing functional calculus on the basis of the neural network.
594. # An important prerequisite for this is to understand the concept of automatic differentiation.
595. # In fact, we have used autodifferentiation already above to realize the adaptation of the neural network parameters during its training.
596. # This has been hidden within the `Flux.train!` method, which evaluates for the given training data the output of the neural network as well as gradients of the loss function with respect to all trainable parameters.
597. # These gradients are used to trace the error of the neural network predictions back to the weights in order to nudge them in the "right direction", i.e. to decrease the loss.
598. #
599. # Reverse mode automatic differentiation computes derivatives by recording individual subexpressions of a given composite function, e.g. of the neural network, during evaluation.
600. # The chain rule is then applied programmatically to yield the derivative of the function in a second computational pass.
601. # As this mechanism is central to machine learning, many frameworks come with ready-to-use implementations.
602. # For instance, Flux.jl provides the `gradient` method, which can be used as follows for scalar-valued functions:
604. # %%
605. f(x, y) = x^2 + sin(x * y) # Some arbitrary test function
606. df(x, y) = (2*x + y * cos(x * y), x * cos(x * y)) # Analytic derivative to test autodiff results
608. x, y = π, 2.0 # Values for which the derivative will be computed
610. println("analytic: ", df(x, y))
611. println("autodiff: ", Flux.gradient(f, x, y))
613. # %% [markdown]
614. # Similarly, we can derive the output of the neural functional with respect to the input density via autodifferentiation.
615. # But what do we get out of this?
616. #
617. # In Part 2, we have seen that the one-body direct correlation functional is obtained from the functional derivative of the excess free energy, $c_1(x; [\rho]) = -\delta \beta F_\mathrm{exc}[\rho] / \delta \rho(x)$, and that it is a central object in DFT calculations.
618. # Applying an additional functional derivative yields another important quantity, the two-body direct correlation function
619. # $$
620. # c_2(x, x'; [\rho]) = \frac{\delta c_1(x; [\rho])}{\delta \rho(x')}.
621. # $$
622. # In fact, a complete hierarchy of direct correlation functions $c_n(x, x', ...; [\rho])$ is defined by iterating subsequent functional derivatives.
623. #
624. # From fundamental physical considerations, further relations can be obtained for these direct correlation functions.
625. # In particular, by applying Noether's theorem to statistical mechanical systems, one arrives at so-called sum rules, which constrain the interrelation of different quantities of interest.
626. # One such sum rule which connects $c_1(x)$, $c_2(x, x')$ and $\rho(x)$ reads as follows:
627. # $$
628. # \nabla c_1(x) = \int \mathrm{d}x' c_2(x, x') \nabla' \rho(x').
629. # $$
630. # As these identities emerge from fundamental invariances of equilibrium many-body systems, they can be used as further valuable tests to check the consistency of the neural functional.
631. #
632. # We implement in the following a method to calculate $c_2(x, x'; [\rho])$ via autodifferentiation of a given function for $c_1(x; [\rho])$.
633. # Note that autodifferentiation can be used very generically, which gives us a way to obtain $c_2(x, x'; [\rho])$ also from the Percus $c_1(x; [\rho])$ without having to derive and implement the analytic expression by hand.
635. # %%
636. function get_c2_autodiff(c1_function, xs)
637. dx = xs[2] - xs[1]
638. function (ρ)
639. Flux.jacobian(c1_function, ρ)[1] / dx
640. end
641. end
643. # %% [markdown]
644. # First, we run some quick tests to illustrate the usage of this method.
645. # As test density input, we take profiles that follow from a randomized continuous external potential.
646. # Run the cell a few times to see different results for $c_2(x, x'; [\rho])$.
647. # Also try to switch out the neural correlation functional for the Percus expression, the use of autodifferentiation remains straightforward.
649. # %%
650. L = 10.0
651. μ, T = 0.0, 1.0
653. xs, ρ = minimize(L, μ, T, generate_Vext(L; num_sin=5, num_lin=0, wall=false), get_c1_Percus)
655. c1_func = get_c1_neural(model, xs)
656. # c1_func = get_c1_Percus(xs) # The Percus c1 can also be used in the following autodifferentiation, no need for a pen-and-paper derivation :)
658. c2_func = get_c2_autodiff(c1_func, xs)
660. c1 = c1_func(ρ)
661. c2 = c2_func(ρ)
663. display(plot(xs, ρ, label="ρ", xlims=(0, L), xlabel="x/σ"))
664. display(plot(xs, c1, label="c1", xlims=(0, L), xlabel="x/σ"))
665. heatmap(xs, xs, c2, colorbar_title="c2", aspect_ratio=1, xlims=(0, L), ylims=(0, L), xlabel="x/σ", ylabel="x'/σ")
667. # %% [markdown]
668. # Now we can check the above Noether sum rule.
669. # The gradients are evaluated numerically on the spatial grid via finite differences, for which we have supplied the function `finite_diff` in [neural.jl](neural.jl).
670. # The following plot shows the right and left hand side of the identity, which coincide when evaluating the quantities with a neural functional that has been trained sufficiently well.
671. # Recall that we have neither specified $c_2(x, x'; [\rho])$ nor the Noether identity during training.
672. # The successful reproduction of such fundamental results implies that the neural functional does not merely interpolate values of $c_1(x; [\rho])$, but that it instead also captures the essential physics of the underlying many-body problem.
674. # %%
675. dx = xs[2] - xs[1]
676. grad_c1 = finite_diff(c1; dx)
677. grad_ρ = finite_diff(ρ; dx)
678. lhs = grad_c1 # Left hand side of the Noether sum rule
679. rhs = c2 * grad_ρ * dx # Right hand side of the Noether sum rule
680. plot(xs, [lhs rhs], label=["∇c1(x)" "∫dx' c2(x,x') ∇'ρ(x')"], linestyle=[:solid :dash], xlims=(0, L), xlabel="x/σ")
682. # %% [markdown]
683. # To complete the functional calculus, we now consider the *inverse* operation of functional differentiation.
684. # This gives us direct access to the excess free energy $F_\mathrm{exc}[\rho]$, as can be deduced from the definition of $c_1(x; [\rho])$ via its functional derivative.
685. # Formally, functional line integration serves as such an inverse operation, and it amounts to performing a line integral in function space.
686. # By choosing a linear parameterization $\rho_a(x) = a \rho(x)$, one can derive the explicit result
687. # $$
688. # \beta F_\mathrm{exc}[\rho] = - \int \mathrm{d}x \rho(x) \int_0^1 \mathrm{d}a c_1(x; [\rho_a]).
689. # $$
690. #
691. # On the level of the neural functional (or, in fact, any representation of the one-body direct correlation functional), this expression is easy to evaluate numerically by discretizing both integrals.
692. # We write a method for this in the following.
694. # %%
695. function get_Fexc_funcintegral(c1_function, xs; num_a=100)
696. dx = xs[2] - xs[1]
697. da = 1 / num_a
698. as = da/2:da:1
699. function (ρ)
700. aintegral = zero(ρ)
701. for a in as
702. aintegral .+= c1_function(a .* ρ)
703. end
704. -sum(ρ .* aintegral) * dx * da
705. end
706. end
708. # %% [markdown]
709. # In the following code cell, randomized inhomogeneous systems are considered as above to test the excess free energy calculation via functional line integration.
710. # We additionally implement the exact analytic Percus expression for $F_\mathrm{exc}[\rho]$ (see Part 2) in order to judge the results.
711. # The values of $F_\mathrm{exc}[\rho]$ as obtained by functional line integration of the neural correlation functional show almost no discrepancy to the exact analytic theory.
713. # %%
714. L = 10.0
715. μ, T = 0.0, 1.0
717. xs, ρ = minimize(L, μ, T, generate_Vext(L), get_c1_Percus)
719. c1_func = get_c1_neural(model, xs)
720. # c1_func = get_c1_Percus(xs)
722. Fexc_func = get_Fexc_funcintegral(c1_func, xs)
724. function get_Fexc_Percus(xs)
725. dx = xs[2] - xs[1]
726. ω0, ω1 = get_weights_Percus(xs)
727. conv(f, g) = conv_fft(f, g; dx)
728. function (ρ)
729. n0, n1 = conv(ω0, ρ), conv(ω1, ρ)
730. ϕ = n0 .* log.(1 .- n1)
731. -sum(ϕ) * dx
732. end
733. end
735. Fexc_Percus = get_Fexc_Percus(xs)
737. println("Fexc from analytic expression (Percus): ", Fexc_Percus(ρ))
738. println("Fexc from functional line integral (neural): ", Fexc_func(ρ))
739. plot(xs, ρ, label="ρ", xlims=(0, L), xlabel="x/σ")
741. # %% [markdown]
742. # ## Conclusions
744. # %% [markdown]
745. # In this tutorial, we have investigated various methods for the description of classical statistical many-body systems.
746. # In particular, we have focused on developing a machine learning framework based on classical DFT and many-body simulations, which combines the strengths of both methods and which facilitates to overcome their individual limitations.
747. #
748. # As preparation, the fundamentals of GCMC simulations were laid out in Part 1.
749. # We have shown that many-body simulations are suitable to obtain equilibrium averages such as the density profile from the microscopic description of a system.
750. # This was exemplified specifically for the hard rod fluid in inhomogeneous environments.
751. # However, the computational demands of simulations often turn out to be too restrictive in concrete applications.
752. #
753. # In Part 2, classical DFT was presented as a coarse-grained method that operates directly on the level of the one-body density rather than on the full many-body picture, which greatly reduces the computational complexity of the problem.
754. # For the case of the hard rod fluid, we have reproduced the exact analytic density functional due to Percus and have shown its application in a standard DFT minimization for the self-consistent calculation of density profiles.
755. # While DFT is formally exact for arbitrary fluid models, we have emphasized the difficulty of finding analytic functionals and the exceptionality of the hard rod result.
756. # These considerations led us to the application of machine learning with the goal of capturing functional maps from simulation data.
757. #
758. # Part 3 gives a detailed account of our machine learning procedure, in which a neural network serves as the central object for the representation of a functional relationship.
759. # Training data has been acquired with GCMC simulations, and we have described in this regard how one can generate suitable external potentials randomly to cover diverse inhomogeneous environments.
760. # The actual training procedure was based on the functional mapping from the density profile $\rho(x)$ to the one-body direct correlation functional $c_1(x; [\rho])$, which were hence taken respectively as input and target output of the neural network.
761. # Crucially, however, we have adopted a *local* learning strategy.
762. # The neural network was constructed to yield only the scalar value of $c_1(x; [\rho])$ at a certain location $x$ when provided with the surrounding density profile within a narrow window around $x$ as input.
763. # This setup imposes no substantial restriction, as the whole profile of $c_1(x; [\rho])$ can be recovered by evaluating the neural network (in parallel) at multiple positions within the system domain.
764. # However, there are numerous benefits to considering the functional relationship locally, such as an improvement of training statistics, the prescription of the short-ranged nature of $c_1(x; [\rho])$ and the applicability of the neural functional to multiscale problems.
765. #
766. # We have then investigated applications of the trained neural correlation functional and its use in DFT minimizations.
767. # The results of this *neural* DFT were highly accurate even for unseen external environments, which could be demonstrated for the hard rod fluid by comparison to the exact Percus solution.
768. # Further, we have shown the implementation of functional calculus on the basis of the neural network, which enabled us to obtain related quantities such as the two-body direct correlation functional $c_2(x, x'; [\rho])$ or the excess free energy $F_\mathrm{exc}[\rho]$.
769. # In this regard, automatic differentiation was utilized for the evaluation of functional derivatives, and a functional line integral served as the explicit inverse operation.
770. # The results of the neural functional calculus could be verified by the exact Percus theory.
771. # A special case, which has not been covered in this tutorial, is the application of the neural functional calculus to systems with spatially constant density, where important bulk properties such as the equation of state and the bulk pair structure can be determined.
772. # It is a valuable excercise to use the above tools for the investigation of these bulk results.
773. # Further details are given in our research paper ([doi:10.1073/pnas.2312484120](https://doi.org/10.1073/pnas.2312484120), [arXiv:2307.04539](https://arxiv.org/abs/2307.04539)).
774. #
775. # As a further consistency check of the neural functional, we have shown the validity of a sum rule which arises from the application of Noether's theorem to statistical many-body systems, and which can hence be traced back to a fundamental invariance property of phase space.
776. # It is therefore nontrivial that the neural network reproduces the identity given that no such information was provided during training.
777. # There are many more sum rules which result from thermal Noether invariance and which are described in more detail in the accompanying manuscript ([arXiv:2312.04681](https://arxiv.org/abs/2312.04681)).
778. # We encourage the reader to implement and verify them; the above code provides the required utilities and may serve as guidance.
779. #
780. # Even though we have focused on the hard rod fluid, which is owed to the pedagogical purpose of this tutorial and to the availability of exact theoretical results, the presented machine learning framework is generically applicable to other fluid models with short-ranged interactions.
781. # You can try this yourself by modifying the above code examples.
782. # Simulation data can be generated for a different type of fluid by simply changing the form of the pair potential $\phi(r)$ and by selecting appropriate parameters for the random generation of chemical and external potentials.
783. # The training routine of the neural functional then proceeds identically, provided that the input range of the density is chosen sufficiently large for the specific interaction type (this could be checked e.g. by evaluating $c_2(x, x'; [\rho])$ via autodifferentiation).
784. # The trained neural functional can be applied as before in DFT minimizations and in the functional calculus, which is particularly relevant for fluid models where no satisfactory analytic functionals exist.
785. # One could even go beyond an isothermal setting and include variations of the temperature, which then serves as a further input parameter to the neural correlation functional.
786. # In total, the neural functional theory leaves much room for further research and we have demonstrated here that it serves as a promising hybrid approach to the investigation of inhomogeneous fluids.

Tutorial.ipynb at commit 357f85d, under MIT · at the source