One of the typical applications of ESPResSo is the simulation of polymer chains with a bead-spring-model. For this we need a repulsive interaction between all beads, for which one usually takes a shifted and truncated Lennard-Jones (so-called Weeks–Chandler–Andersen or WCA) interaction, and additionally a bonded interaction between adjacent beads to hold the polymer together. You have already learned that the command
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=1.0, sigma=1.0,
shift=0.25, cutof=1.226)
creates a Lennard-Jones interaction with $\varepsilon=1.$, $\sigma=1.$, $r_\text{cut} = 1.125$ and $\varepsilon_\text{shift}=0.25$ between particles of type 0, which is the desired repulsive interaction. The command
from espressomd import interactions
fene = interactions.FeneBond(k=7, d_r_max=2)
creates a FeneBond object (see ESPResSo manual for the details). Still ESPResSo does not know between which beads this interaction should be applied. This can be either be specified explicitly or done with the polymer module. This creates a given number of beads, links them with the given bonded interaction and places them following a certain algorithm. We will use the pruned self-avoiding walk to place the monomers with a fixed distance between adjacent bead positions. The syntax is:
from espressomd import polymer
# mpc: monomers per chain
mpc = 30
poly = polymer.Polymer(N_P=1, MPC=mpc, bond=fene, bond_length=1)
Compared to a random walk, the pruned self-avoiding walk limits the risk of the polymer crossing itself (or closely approaching itself). Since the LJ potential is very steep, this would raise the potential energy enormously and would make the monomers shoot through the simulation box. The pruned self-avoiding walk should prevent that, but to be sure we perform some MD steps with a capped LJ potential, i.e. forces above a certain threshold will be set to the threshold to prevent the system from exploding. The script below contains quite a long warmup command to allow longer polymers.
This allows to quickly change the number of monomers without editing the script. For the warmup a Langevin thermostat is used to keep the temperature constant. Furthermore we want to compute the diffusion constant of the polymer for different numbers of monomers. For this purpose we can again use the multiple tau correlator. Have a look at the ESPResSo-script for the single particle diffusion and add the adapted commands for the polymer. Find out how many integration steps are necessary to capture the long-time diffusion regime of the polymer. The script already computes the time averaged hydrodynamic radius and stores it in a file rh_out.dat whose first column denotes the number of monomers.
In [ ]:
import espressomd
espressomd.assert_features(['CUDA', 'LENNARD_JONES'])
from espressomd import System, interactions, lb, polymer
from espressomd.observables import ComPosition
from espressomd.accumulators import Correlator
from numpy import savetxt, zeros
import numpy as np
# Setup constant
time_step = 0.01
loops = 100
step_per_loop = 100
# System setup
system = System(box_l=[32.0, 32.0, 32.0])
system.set_random_state_PRNG()
np.random.seed(seed=system.seed)
system.cell_system.skin = 0.4
mpc = 20 # The number of monomers has been set to be 20 as default
# Change this value for further simulations
# Lennard-Jones interaction
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=1.0, sigma=1.0,
shift="auto", cutoff=2.0**(1.0 / 6.0))
# Fene interaction
fene = interactions.FeneBond(k=7, d_r_max=2)
system.bonded_inter.add(fene)
# Setup polymer of part_id 0 with fene bond
positions = polymer.positions(n_polymers=1,
beads_per_chain=mpc,
bond_length=1, seed=5642,
min_distance=0.9)
for i, pos in enumerate(positions[0]):
id = len(system.part)
system.part.add(id=id, pos=pos)
if i > 0:
system.part[id].add_bond((fene, id - 1))
print("Warming up the polymer chain.")
# For longer chains (>100) an extensive
# warmup is necessary ...
system.time_step = 0.002
system.thermostat.set_langevin(kT=1.0, gamma=10, seed=42)
for i in range(100):
system.force_cap = float(i) + 1
system.integrator.run(1000)
print("Warmup finished.")
system.force_cap = 0
system.integrator.run(10000)
system.time_step = time_step
system.integrator.run(50000)
system.thermostat.turn_off()
system.part[:].v = [0, 0, 0]
lbf = lb.LBFluidGPU(kT=1, seed=123, agrid=1, dens=1, visc=5, tau=time_step)
system.actors.add(lbf)
system.thermostat.set_lb(LB_fluid=lbf, seed=142, gamma=5)
print("Warming up the system with LB fluid.")
system.integrator.run(1000)
print("LB fluid warming finished.")
# configure correlators
com_pos = ComPosition(ids=(0,))
c = Correlator(obs1=com_pos, tau_lin=16, tau_max=loops * step_per_loop, delta_N=1,
corr_operation="square_distance_componentwise", compress1="discard1")
system.auto_update_accumulators.add(c)
print("Sampling started.")
for i in range(loops):
system.integrator.run(step_per_loop)
system.analysis.append()
print("\rSampling: {:.0f}%".format(i * 100. / loops), end="")
print("\rSampling finished.")
c.finalize()
corrdata = c.result()
corr = zeros((corrdata.shape[0], 2))
corr[:, 0] = corrdata[:, 0]
corr[:, 1] = (corrdata[:, 2] + corrdata[:, 3] + corrdata[:, 4]) / 3
savetxt("./msd_nom" + str(mpc) + ".dat", corr)
with open("./rh_out.dat", "a") as datafile:
rh = system.analysis.calc_rh(chain_start=0, number_of_chains=1, chain_length=mpc - 1)
datafile.write(str(mpc) + " " + str(rh[0]) + "\n")
Run the script for different numbers of monomers and determine the evolution of the diffusion coefficient as a function of the chain length. Compare the results of your ESPResSo simulations with the given Kirkwood-Zimm formula.
[1] S Succi. The lattice Boltzmann equation for fluid dynamics and beyond. Clarendon Press, Oxford, 2001.
[2] B. Dünweg and A. J. C. Ladd. Advanced Computer Simulation Approaches for Soft Matter Sciences III, chapter II, pages 89–166. Springer, 2009.
[3] B. Dünweg, U. Schiller, and A.J.C. Ladd. Statistical mechanics of the fluctuating lattice-boltzmann equation. Phys. Rev. E, 76:36704, 2007.
[4] P. G. de Gennes. Scaling Concepts in Polymer Physics. Cornell University Press, Ithaca, NY, 1979.
[5] M. Doi. Introduction to Polymer Physics. Clarendon Press, Oxford, 1996.
[6] Michael Rubinstein and Ralph H. Colby. Polymer Physics. Oxford University Press, Oxford, UK, 2003.
[7] Daan Frenkel and Berend Smit. Understanding Molecular Simulation. Academic Press, San Diego, second edition, 2002.