-
Notifications
You must be signed in to change notification settings - Fork 2
/
urpm_scan.py
executable file
·130 lines (97 loc) · 5.19 KB
/
urpm_scan.py
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
#!/usr/bin/env python3
# This file is part of SunlightDPD - a home for open source software
# related to the dissipative particle dynamics (DPD) simulation
# method.
# Copyright (c) 2009-2017 Unilever UK Central Resources Ltd
# (Registered in England & Wales, Company No 29140; Registered
# Office: Unilever House, Blackfriars, London, EC4P 4BQ, UK).
# SunlightDPD 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.
# SunlightDPD 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 SunlightDPD. If not, see <http://www.gnu.org/licenses/>.
# Scan the screening behaviour of various URPM models, using command
# line parameters as a control.
# Typical use starts off with a broad (--lo, --hi) interval where the
# numbers are -log10(rho_z). Then the interval which contains
# the Kirkwood crossover between pure exponential and oscillatory
# decay can be narrowed by a factor 10 each time. For example :
# ./urpm_scan.py --lo=-2.0 --hi=-1.0
# ./urpm_scan.py --lo=-1.6 --hi=-1.5
# ./urpm_scan.py --lo=-1.53 --hi=-1.52
# Inspecting the output of the last narrows the interval to (1.524,
# 1.526) - this for the default parameters.
import argparse
import math as m
import matplotlib.pyplot as plt
from oz import wizard as w
parser = argparse.ArgumentParser(description='URPM h(r) scanner')
parser.add_argument('--ncomp', action='store', default=2, type=int, help='number of components (species) (default 2)')
parser.add_argument('--ng', action='store', default=4096, type=int, help='number of grid points (default 4096)')
parser.add_argument('--deltar', action='store', default=0.01, type=float, help='grid spacing (default 0.01)')
parser.add_argument('--alpha', action='store', default=0.2, type=float, help='Picard mixing fraction (default 0.2)')
parser.add_argument('--npic', action='store', default=6, type=int, help='number of Picard steps (default 6)')
parser.add_argument('--lb', action='store', default=1.0, type=float, help='Bjerrum length (default 1.0)')
parser.add_argument('--sigma', action='store', default=1.0, type=float, help='Gaussian charge size (default 1.0)')
parser.add_argument('--R', action='store', default=0.0, type=float, help='Grootian charge size (defaults to sigma*sqrt(15/2))')
parser.add_argument('--rc', action='store', default=1.0, type=float, help='DPD length scale (default 1.0)')
parser.add_argument('--A', action='store', default=0.0, type=float, help='DPD repulsion amplitude (default 0.0)')
parser.add_argument('--z1', action='store', default=1, type=int, help='valency of positive ions (default +1)')
parser.add_argument('--z2', action='store', default=-1, type=int, help='valency of negative ions (default -1)')
parser.add_argument('--type', action='store', default=1, type=int, help='charge type (1=Gaussian, 2=Bessel, 3=Groot, 4=Slater)')
parser.add_argument('--case', action='store', default=1, type=int, help='Slater method (1=exact, 2=good, 3=bad)')
parser.add_argument('--rho', action='store', default=3.0, type=float, help='total density if ncomp = 3 (default 3.0)')
parser.add_argument('--npt', action='store', default=11, type=int, help='number of rho points (default 11)')
parser.add_argument('--lo', action='store', default=-2.0, type=float, help='log10(rho_min) (default -2.0)')
parser.add_argument('--hi', action='store', default=-1.0, type=float, help='log10(rho_max) (default -1.0)')
args = parser.parse_args()
w.ng = args.ng
w.ncomp = args.ncomp
w.deltar = args.deltar
w.alpha = args.alpha
w.npic = args.npic
w.initialise()
w.lb = args.lb
w.sigma = args.sigma
if (args.R > 0): w.rgroot = args.R
else: w.rgroot = args.sigma * m.sqrt(7.5)
w.rc = args.rc
w.arep[:,:] = args.A
w.z[0] = args.z1
w.z[1] = args.z2
# potential type = 4 (exact), or potential type = 5 (approximate) with
# beta = 1/lambda, or beta=5/(8lambda)
if args.type < 4:
w.dpd_potential(args.type)
else:
if args.case == 1:
w.dpd_potential(w.dpd_slater_exact_charges)
else:
if args.case == 2: w.beta = 5 / (8*w.lbda)
else: w.beta = 1 / w.lbda
w.dpd_potential(w.dpd_slater_approx_charges)
n = args.npt
off = -2
eps = 1e-20
plt.figure(1)
for i in range(n):
log10rho = args.lo + (args.hi - args.lo) * i / (n - 1.0)
rhoz = 10.0**(log10rho)
w.rho[0] = rhoz / (args.z1 * (args.z1 - args.z2))
w.rho[1] = rhoz / (args.z2 * (args.z2 - args.z1))
if (w.ncomp > 2): w.rho[2] = args.rho - w.rho[0] - w.rho[1]
if (i == 0): w.write_params()
w.hnc_solve()
if w.return_code: exit()
if (w.ncomp == 2): print('%i\t%g\t%g\t%g' % (i, log10rho, rhoz, w.error))
else: print('%i\t%g\t%g\t%g\t%g' % (i, log10rho, rhoz, args.rho, w.error))
plt.plot(w.r[:], list(map(lambda x, y: i*off + m.log10(eps + m.fabs(x*y)), w.hr[:, 0, 0], w.r[:])), label="%g" % log10rho)
plt.legend(loc='upper right')
plt.xlabel('$r$')
plt.ylabel('$\log_{10} r\,h_{++}$')
plt.show()