-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathhs_diameter.py
More file actions
126 lines (107 loc) · 4.34 KB
/
hs_diameter.py
File metadata and controls
126 lines (107 loc) · 4.34 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
"""Module for effective hard-sphere diameters"""
import numpy as np
import scipy.optimize as SciOpt
import scipy.integrate as SciInt
def calc_dhs_bh(pot, T):
'''Calculate Barker-Henderson hard-sphere diameter [m].'''
if pot.hardcore:
return pot.sigma
sigma = pot.sigmaeff
def integrand(w):
Phi = pot.calc_pot_divk(sigma*w)
return (1-np.exp(-Phi/T))
# Integration using the substitution w = r/sigma -> dr = sigma*dw
integral, error = SciInt.quad(func=integrand, a=0, b=1.0,
epsabs=sigma*1e-10, limit=10000)
assert error/integral<1e-4, error
dhs = integral*sigma
return dhs
def calc_dhs_rep(pot, T):
'''A hard-sphere diameter calculated from the repulsive region, Eq. 4 in Noro-Frenkel paper (m)'''
if pot.hardcore:
return pot.sigma
sigma = pot.sigmaeff
epsdivk = pot.epsdivkeff
def integrand(w):
Phi = pot.calc_pot_divk(sigma*w) + epsdivk
return (1-np.exp(-Phi/T))
# Integration using the substitution w = r/sigma -> dr = sigma*dw
integral, error = SciInt.quad(func=integrand, a=0, b=pot.rmin/sigma,
epsabs=sigma*1e-10, limit=10000)
assert error/integral<1e-4, error
dhs = integral*sigma
return dhs
def calc_dhs_wca(pot, T, rho):
'''Calculate Weeks-Chandler-Andersen hard-sphere diameter [m].'''
if pot.hardcore:
return pot.sigma
sigma = pot.sigmaeff
epsdivk = pot.epsdivkeff
rmin = pot.rmin
def resid(d):
if d>rmin:
return np.inf
if d<=0:
return -np.inf
def integrand_left(r):
Phi = np.exp(-(pot.calc_pot_divk(r)+epsdivk)/T)
return Phi*y_hs_desousa_amotz(r, rho, d)
def integrand_right(r):
Phi_m1 = np.exp(-(pot.calc_pot_divk(r)+epsdivk)/T) - 1
return Phi_m1*y_hs_desousa_amotz(r, rho, d)
integral_left, error = SciInt.quad(func=integrand_left, a=1e-10*sigma, b=d,
epsabs=sigma*1e-10, limit=10000)
integral_right, error = SciInt.quad(func=integrand_right, a=d, b=rmin,
epsabs=sigma*1e-10, limit=10000)
return integral_left + integral_right
sol = SciOpt.root(resid, x0=sigma)
return sol.x[0]
def calc_dhs_wca_lado(pot, T, rho):
'''Calculate Weeks-Chandler-Andersen hard-sphere diameter following Lado [m].'''
if pot.hardcore:
return pot.sigma
sigma = pot.sigmaeff
epsdivk = pot.epsdivkeff
rmin = pot.rmin
def resid(d):
if d<0.1*sigma:
return np.inf
def integrand_left(r):
Phi = np.exp(-(pot.calc_pot_divk(r)+epsdivk)/T)
return np.exp(-(pot.calc_pot_divk(r)+epsdivk)/T)*dy_hs_desousa_amotz_ddhs(r, rho, d)
def integrand_right(r):
Phi_m1 = np.exp(-(pot.calc_pot_divk(r)+epsdivk)/T) - 1
return Phi_m1*dy_hs_desousa_amotz_ddhs(r, rho, d)
integral_left, error = SciInt.quad(func=integrand_left, a=1e-10*sigma, b=d,
epsabs=sigma*1e-10, limit=10000)
integral_right, error = SciInt.quad(func=integrand_right, a=d, b=rmin,
epsabs=sigma*1e-10, limit=10000)
return integral_left + integral_right
sol = SciOpt.root(resid, x0=sigma)
return sol.x[0]
def calc_dhs_BFC(pot, T):
'''Calculate hard-sphere diameter according to Boltzmann-Factor-Criterion [m].'''
if pot.hardcore:
return pot.sigma
def resid(r):
return pot.calc_pot_divk(r)/T-1
sol = SciOpt.root(resid, x0=pot.sigma)
return sol.x[0]
def y_hs_desousa_amotz(r, rho, dhs):
'''Analytic approximation of hard sphere cavity correlation function y
at position r, valid for pure hard-spheres. The approximation is
due to de Souza and Ben-Amotz (10.1080/00268979300100131)
'''
# if r>1.5*dhs:
# return 1.0
eta = np.pi/6 * rho * dhs**3
denom = (1.0-eta)**3
A = (3.0-eta)/denom - 3.0
B = -3*eta*(2.0-eta)/denom
C = np.log((2.0-eta)*(2*denom)) - eta*(2.0-6*eta + 3*eta**2)/denom
y = np.exp(A + B*(r/dhs) + C*(r/dhs)**3)
return y
def dy_hs_desousa_amotz_ddhs(r, rho, dhs):
delta = dhs*1e-5
numdiff = (y_hs_desousa_amotz(r, rho, dhs+delta) - y_hs_desousa_amotz(r, rho, dhs-delta))/(2*delta)
return numdiff