-
Notifications
You must be signed in to change notification settings - Fork 2
Expand file tree
/
Copy pathdpd_demo.py
More file actions
executable file
·216 lines (164 loc) · 8.09 KB
/
dpd_demo.py
File metadata and controls
executable file
·216 lines (164 loc) · 8.09 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# This program is part of pyHNC, copyright (c) 2023 Patrick B Warren (STFC).
# Email: patrick.warren{at}stfc.ac.uk.
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with this program. If not, see
# <http://www.gnu.org/licenses/>.
# Demonstrate the capabilities of the HNC package for solving DPD
# potentials, comparing with SunlightHNC if requested, and plotting
# the pair distribution function and the structure factor too. For
# details here see also the SunlightHNC documentation.
# For standard DPD at A = 25 and ρ = 3, we have the following table
# ∆t = 0.02 ∆t = 0.01 Monte-Carlo HNC deviation
# pressure 23.73±0.02 23.69±0.02 23.65±0.02 23.564 (0.4%)
# energy 13.66±0.02 13.64±0.02 13.63±0.02 13.762 (1.0%)
# mu^ex 12.14±0.02 12.16±0.02 12.25±0.10 12.170 (0.7%)
# The first two columns are from dynamic simulations. The excess
# chemical potential (final row) is measured by Widom insertion. The
# HNC results from the present code are in agreement with those from
# SunlightHNC to at least the indicated number of decimals. The
# deviation is between HNC and simulation results.
# Data is from a forthcoming publication on osmotic pressure in DPD.
import os
import pyHNC
import argparse
import numpy as np
from numpy import pi as π
from pyHNC import Grid, PicardHNC, truncate_to_zero
parser = argparse.ArgumentParser(description='DPD HNC calculator')
pyHNC.add_grid_args(parser)
pyHNC.add_solver_args(parser)
parser.add_argument('-v', '--verbose', action='count', help='more details (repeat as required)')
parser.add_argument('-A', '--A', default=25.0, type=float, help='repulsion amplitude, default 25.0')
parser.add_argument('-r', '--rho', default=3.0, type=float, help='density, default 3.0')
parser.add_argument('--dlambda', default=0.05, type=float, help='for coupling constant integration, default 0.05')
parser.add_argument('--rcut', default=3.0, type=float, help='maximum in r for plotting, default 3.0')
parser.add_argument('--qcut', default=25.0, type=float, help='maximum in q for plotting, default 25.0')
parser.add_argument('--sunlight', action='store_true', help='compare to SunlightHNC')
parser.add_argument('--rpa', action='store_true', help='solve RPA instead of HNC')
parser.add_argument('--exp', action='store_true', help='solve EXP instead of HNC')
parser.add_argument('-s', '--show', action='store_true', help='show results')
parser.add_argument('-o', '--output', help='write pair function to a file')
args = parser.parse_args()
args.script = os.path.basename(__file__)
A, ρ = args.A, args.rho
grid = Grid(**pyHNC.grid_args(args)) # make the initial working grid
r, Δr, q = grid.r, grid.deltar, grid.q # extract the co-ordinate arrays for use below
if args.verbose:
print(f'{args.script}: {grid.details}')
# Define the DPD potential, and its derivative, then solve the HNC
# problem. The arrays here are all size ng-1, same as r[:]
φ = truncate_to_zero(A/2*(1-r)**2, r, 1) # the DPD potential
f = truncate_to_zero(A*(1-r), r, 1) # the force f = -dφ/dr
solver = PicardHNC(grid, **pyHNC.solver_args(args))
if args.verbose:
print(f'{args.script}: {solver.details}')
soln = solver.solve(φ, ρ, monitor=args.verbose) # solve for the DPD potential
h, c, hq = soln.hr, soln.cr, soln.hq # extract for use in a moment
# For the integrals here, see Eqs. (2.5.20) and (2.5.22) in Hansen &
# McDonald, "Theory of Simple Liquids" (3rd edition): for the (excess)
# energy density, e = 2πρ² ∫_0^∞ dr r² φ(r) g(r) and virial pressure,
# p = ρ + 2πρ²/3 ∫_0^∞ dr r³ f(r) g(r) where f(r) = −dφ/dr is the
# force. An integration by parts shows that the mean-field
# contributions, being these with g(r) = 1, are the same.
# Here specifically the mean-field contributions are
# 2πρ²/3 ∫_0^∞ dr r³ f(r) = A ∫_0^1 dr r³ (1−r) = πAρ²/30 .
e_mf = p_mf = π*A*ρ**2/30
e_xc = 2*π*ρ**2 * np.trapz(r**2*φ*h, dx=Δr)
e_ex = e_mf + e_xc
e = 3*ρ/2 + e_ex
p_xc = 2*π*ρ**2/3 * np.trapz(r**3*f*h, dx=Δr)
p_ex = p_mf + p_xc
p = ρ + p_ex
μ_ex = 4*π*ρ * np.trapz(r**2*(h*(h-c)/2 - c), dx=Δr)
μ = np.log(ρ) + μ_ex
# Coupling constant integration for the free energy
# Function to calculate the excess non-mean-field energy with coupling
# constant λ. Uses a bunch of main script variables that are 'in
# scope' here.
def excess(λ):
'''Return the excess correlation energy with coupling λ'''
h = 0 if λ == 0 else solver.solve(λ*φ, ρ).hr # presumed will converge !
e_xc = 2*π*ρ**2 * np.trapz(r**2*φ*h, dx=Δr) # the integral above
if args.verbose and args.verbose > 1:
print(f'{args.script}: excess: λ = {λ:0.3f}, e_xc = {e_xc:f}')
return e_xc
λ_arr = np.linspace(0, 1, 1+round(1/args.dlambda))
dλ = pyHNC.grid_spacing(λ_arr)
e_xc_arr = np.array([excess(λ) for λ in np.flip(λ_arr)]) # descend, to assure convergence
f_xc = np.trapz(e_xc_arr, dx=dλ) # the coupling constant integral
f_ex = e_mf + f_xc # f_mf = e_mf in this case
closure = 'EXP' if args.exp else 'RPA' if args.rpa else 'HNC'
print(f'{args.script}: model: standard DPD with A = {A:g}, ρ = {ρ:g}, {closure} closure')
print(f'{args.script}: Monte-Carlo (A,ρ = 25,3): energy, virial pressure =',
'\t13.63±0.02\t\t\t23.65±0.02')
print(f'{args.script}: pyHNC v{pyHNC.version}: energy, free energy, virial pressure =',
'\t%0.5f\t%0.5f\t%0.5f' % (e_ex, f_ex, p))
print(f'{args.script}: pyHNC v{pyHNC.version}: chemical potential, free energy =',
'\t%0.5f\t%0.5f' % (μ, ρ*μ_ex - p_ex))
if args.rpa or args.exp:
c = -φ # the RPA
cq = grid.fourier_bessel_forward(c) # forward transform to reciprocal space
hq = cq / (1 - ρ*cq) # solve the OZ relation
h = grid.fourier_bessel_backward(hq) # back transform to real space
if args.exp: # implement EXP if requested
h = (np.exp(h) - 1)
if args.sunlight:
from oz import wizard as w
w.ng = grid.ng
w.deltar = grid.deltar
w.initialise()
w.arep[0,0] = A
w.dpd_potential()
w.rho[0] = ρ
w.rpa_solve() if args.rpa else w.hnc_solve()
sunlight_version = str(w.version, 'utf-8').strip()
print(f'{args.script}: SunlightHNC v{sunlight_version}: energy, free energy, virial pressure =',
'\t%0.5f\t%0.5f\t%0.5f' % (w.uex, w.aex, w.press))
print(f'{args.script}: SunlightHNC v{sunlight_version}: chemical potential =',
'\t%0.5f' % (np.log(ρ)+w.muex[0]))
if args.show:
import matplotlib.pyplot as plt
gr = 1.0 + h # the pair function
sq = 1.0 + ρ*hq # the structure factor
plt.figure(1)
cut = r < args.rcut
if args.sunlight:
imax = int(args.rcut / w.deltar)
plt.plot(w.r[0:imax], 1.0+w.hr[0:imax,0,0], '.')
plt.plot(r[cut], gr[cut], '--')
else:
plt.plot(r[cut], gr[cut])
plt.xlabel('$r$')
plt.ylabel('$g(r)$')
plt.figure(2)
cut = q < args.qcut
if args.sunlight:
jmax = int(args.qcut / w.deltak)
plt.plot(w.k[0:jmax], w.sk[0:jmax,0,0]/ρ, '.')
plt.plot(q[cut], sq[cut], '--')
else:
plt.plot(q[cut], sq[cut])
plt.xlabel('$k$')
plt.ylabel('$S(k)$')
plt.show()
if args.output:
import pandas as pd
g = 1.0 + h # the pair function
cut = r < args.rcut
df = pd.DataFrame({'r': r[cut], 'g': g[cut]})
df_agr = pyHNC.df_to_agr(df[['r', 'g']])
with open(args.output, 'w') as f:
f.write(f'# DPD with A = {A:g}, ρ = {ρ:g}, {closure} closure\n')
f.write(df_agr) # use a utility here to convert to xmgrace format
f.write('\n')
print(f'{args.script}: written (r, g) to {args.output}')