wsg_fig5_compare.py

#!/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.

# Based on an original code copyright (c) 2007 Lucian Anton.
# Modifications copyright (c) 2008, 2009 Andrey Vlasov.  Additional
# modifications 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).
# Additional modifications copyright (c) 2020-2021 Patrick B Warren
# <patrick.warren@stfc.ac.uk> and STFC.

# 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/>.

# The results can be directly compared with Fig 5 of Wijmans et al
# [J. Chem. Phys. v114, 7644 (2001)].  The chemical potential
# difference muex_1 - muex_0 is exactly equal to the free energy
# change to convert a particle of species 0 into a particle of species
# 1.  Note that the calculation is done as a two component problem
# with a mole fraction of the second species set to zero.

# Here is the data from Fig 5 of Wijmans et al.  

#          dA dF

data1 = [[ 0, 0],
         [ 2, 0.699745],
         [ 4, 1.370293],
         [ 6, 2.0262583],
         [ 8, 2.6384122],
         [10, 3.2359834],
         [12, 3.7897591],
         [14, 4.3435348],
         [16, 4.853515],
         [18, 5.3488966],
         [20, 5.8296797],
         [22, 6.2812658],
         [24, 6.7182534],
         [26, 7.1406424],
         [28, 7.5338344],
         [30, 7.9270106]]

wsgx = list(data1[i][0] for i in range(len(data1)))
wsgy = list(data1[i][1] for i in range(len(data1)))

# Here is another set of data generated by calculations of the
# chemical potential difference, using Widom insertion.

#    dA   d(mu)   std-err

data2 = [[0, 0.0007, 0.00299772],
         [1, 0.3571, 0.00355517],
         [2, 0.7052, 0.004213],
         [3, 1.0451, 0.00498353],
         [4, 1.3771, 0.00587967],
         [5, 1.7014, 0.00690962],
         [6, 2.0181, 0.00808128],
         [7, 2.3276, 0.00940074],
         [8, 2.6299, 0.0108687],
         [9, 2.9254, 0.0124841],
         [10, 3.2143, 0.014242],
         [11, 3.4969, 0.0161377],
         [12, 3.7734, 0.0181624],
         [13, 4.0441, 0.0203056],
         [14, 4.3093, 0.0225574],
         [15, 4.5692, 0.0249076],
         [16, 4.824, 0.0273447],
         [17, 5.0741, 0.0298606],
         [18, 5.3197, 0.0324429],
         [19, 5.5609, 0.0350847],
         [20, 5.7982, 0.0377775],
         [21, 6.0317, 0.0405129],
         [22, 6.2615, 0.0432836],
         [23, 6.488, 0.0460852],
         [24, 6.7113, 0.0489117],
         [25, 6.9317, 0.0517575],
         [26, 7.1492, 0.0546188],
         [27, 7.3641, 0.0574919],
         [28, 7.5765, 0.0603722],
         [29, 7.7866, 0.0632587],
         [30, 7.9945, 0.0661484]]

mcx = list(data2[i][0] for i in range(len(data2)))
mcy = list(data2[i][1] for i in range(len(data2)))
mce = list(data2[i][2] for i in range(len(data2)))

from oz import wizard as w

w.ncomp = 2
w.initialise()

rho = 3.0
A = 25.0
Amax = 55.0

w.arep[0,0] = w.arep[1,1] = A
w.rho[0] = rho
w.rho[1] = 0.0

npt = 41

x = []
y = []

print('#dA\tmuex0\tmuex1\tdmuex\terror')

for i in range(npt):
    dA = (Amax - A) * i / (npt - 1.0)
    w.arep[0,1] = A + dA
    w.dpd_potential()
    w.hnc_solve()
    x.append(dA)
    y.append(w.muex[1]-w.muex[0])
    print("%f\t%f\t%f\t%f\t%g" % (dA, w.muex[0], w.muex[1], w.muex[1]-w.muex[0], w.error))

import matplotlib.pyplot as plt

# plt.errorbar(mcx, mcy, yerr=mce, fmt='co', label='Monte-Carlo')
# error bars are smaller than symbols so just plot data

plt.plot(mcx, mcy, 'co', label='Widom insertion')
plt.plot(wsgx, wsgy, 'ro', label='Wijmans et al. (2001)')
plt.plot(x, y, label='HNC')
plt.xlabel('$\\Delta A$')
plt.ylabel('$\\Delta F$')
plt.legend(loc='lower right')

plt.show()