<jupyter_start><jupyter_text>Penalty FunctionThis is a simple example of how to use the flat histogram method in in Faunus<jupyter_code>%matplotlib inline
import matplotlib
import matplotlib.cm as cm
import numpy as np
import matplotlib.pyplot as plt
from math import log, fabs, pi, cos, sin<jupyter_output><empty_output><jupyter_text>Run penalty example<jupyter_code>%%bash
rm -fR penalty.dat
if [[ -z "${FAUNUS_EXECUTABLE}" ]]; then
  yason.py penalty.yml | faunus --nobar
else
  echo "Seems we're running CTest - use Faunus target from CMake"
  "${YASON_EXECUTABLE}" penalty.yml | "${FAUNUS_EXECUTABLE}" --nobar
fi<jupyter_output><empty_output><jupyter_text>Load penalty functionThis loads the penalty function from disk which is simply the negative of the final potential of mean force. Also load the final histogram, which ideally should be _flat_ over the reaction coordinates.<jupyter_code>PMF = -np.loadtxt('penalty.dat', skiprows=1).T
HIST = np.loadtxt('penalty-histogram.dat').T
HIST = np.log(HIST/HIST.min())
print('Penalty maximum barrier    =', fabs(PMF.max()-PMF.min()), 'kT')
print('Maximum simulation barrier =', HIST.max(), 'kT')
fig, (ax1, ax2) = plt.subplots(1,2, sharex=True, sharey=True)
fig.set_size_inches(10,4)
im1 = ax1.imshow(PMF, interpolation='bilinear', cmap=cm.afmhot, extent=[-2,2,-2,2],
               origin='lower', vmin=PMF.min(), vmax=PMF.max())
im2 = ax2.imshow(HIST, interpolation='bilinear', cmap=cm.afmhot, extent=[-2,2,-2,2],
               origin='lower', vmin=HIST.min(), vmax=HIST.max())
ax1.set_title('Potential of Mean Force')
ax2.set_title('Barriers during simulation')
ax1.set_xlabel(r'$x$')
ax2.set_xlabel(r'$x$')
ax1.set_ylabel(r'$y$')
ax2.set_ylabel(r'$y$')
cbar1 = fig.colorbar(im1, ax=ax1)
cbar2 = fig.colorbar(im2, ax=ax2)
cbar1.set_label('$k_BT$')
cbar2.set_label('$k_BT$')
plt.show()<jupyter_output><empty_output><jupyter_text>External PotentialNow let's plot the external potential applied during simulation. The PMF should tend towards this.<jupyter_code>def energy(x,y):
    ''' Energy function used during simulation '''
    s=1+sin(2*pi*x) + cos(2*pi*y)
    if x>=-2.00 and x<=-1.25: return 1*s
    if x>=-1.25 and x<=-0.25: return 2*s
    if x>=-0.25 and x<= 0.75: return 3*s
    if x>= 0.75 and x<= 1.75: return 4*s
    if x>= 1.75 and x<= 2.00: return 5*s
    raise Exception
x = y = np.arange(-2.0, 2.0, 0.1)
X, Y = np.meshgrid(x,y) 
Z = np.vectorize(energy)(X,Y)
Z = Z-Z.max()
fig, ax = plt.subplots()
im = ax.imshow(Z, interpolation='bilinear', cmap=cm.afmhot,
               origin='lower', extent=[-2, 2, -2, 2],
               vmax=Z.max(), vmin=Z.min())
ax.set_title('External potential')
fig.colorbar(im, ax=ax)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1618134906008109, 'fraction_numerical': 0.030224843346848506, 'mean_word_length': 4.78678038379531, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '466366', 'license': 'MIT', 'n_tokens_mistral': 1102, 'n_tokens_neox': 1055, 'n_words': 252, 'path': 'examples/penalty/penalty.ipynb', 'repo_name': 'MarcoPolimeni/faunus', 'size': '2713'}