Decomposition of electrostatic interaction energy between two multipolar molecules

This is a Jupyter notebook (http://jupyter.org) for studying the interaction between multipolar particles.

Load required packages



In [ ]:

    
%matplotlib inline
import matplotlib
import matplotlib.cm as cm
import unittest
import numpy as np
import matplotlib.pyplot as plt
import sys
plt.rcParams.update({'font.size': 18, 'figure.figsize': [8.0, 6.0]})

Run simulation



In [ ]:

    
%%bash
if [[ -z "${FAUNUS_EXECUTABLE}" ]]; then
  yason.py multipole.yml | faunus --nobar -s multipole.state.ubj
else
  echo "Seems we're running CTest - use Faunus target from CMake"
  "${YASON_EXECUTABLE}" multipole.yml | "${FAUNUS_EXECUTABLE}" --nobar -s multipole.state.ubj
fi

Plot multipolar energies as a function of separation



In [ ]:

    
R, exact, total, ionion, iondip, dipdip, ionquad, mucorr = np.loadtxt('multipole.dat', unpack=True, skiprows=2)
lw=4
plt.plot(R, ionion, label='ion-ion', lw=lw)
plt.plot(R, iondip, label='ion-dipole', lw=lw)
plt.plot(R, dipdip, label='dipole-dipole', lw=lw)
plt.plot(R, ionquad, label='ion-quadrupole', lw=lw)
plt.plot(R, total, label='sum of multipoles', lw=lw)
plt.plot(R, exact, 'ko', label='exact')
plt.xlabel('separation (Å)')
plt.ylabel('energy ($kT/\lambda_B$)')
plt.legend(loc=0, frameon=False)
plt.savefig('multipole.png', bbox_inches='tight')

Unittests

Compare distributions with previously saved results



In [ ]:

    
class TestMultipole(unittest.TestCase):
    def test_Exact(self):
        self.assertAlmostEqual(exact.mean(), -0.12266326530612245, places=3)

    def test_IonIon(self):
        self.assertAlmostEqual(ionion.mean(), -0.11624285714285715, places=2)

    def test_IonDipole(self):
        self.assertAlmostEqual(iondip.mean(), -0.006777551020408164, places=2)
        
    def test_DipoleDipole(self):
        self.assertAlmostEqual(dipdip.mean(), -0.0008632653061224489, places=2)
        
    def test_IonQuadrupole(self):
        self.assertAlmostEqual(ionquad.mean(), 0.0005040816326530612, places=3)

test = TestMultipole()
suite = unittest.TestLoader().loadTestsFromModule(test)
ret = unittest.TextTestRunner(verbosity=2).run(suite)
if (not ret.wasSuccessful()):
    sys.exit(ret)