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) 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: The monomers are set according to a pruned self-avoiding walk (in 3D) 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)Using a random walk to create a polymer causes trouble: The random walk may cross itself (or closely approach itself) and 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, this means forces above a certain threshold will be set to the threshold in order to prevent the system from exploding. To see how this is done, look at the script polymer_diffusion.py (see the code below). It contains a quite long warmup command so that also longer polymers are possible. You can probably make it shorter.
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 [ ]:
from __future__ import print_function, division
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
import sys
# 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
polymer.create_polymer(N_P=1, MPC=mpc, bond=fene, bond_length=1,
start_pos = [16.0, 16.0, 16.0])
print("Warming up the polymer chain.")
## For longer chains (>100) an extensive
## warmup is neccessary ...
system.time_step = 0.002
system.thermostat.set_langevin(kT=1.0, gamma=10)
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(agrid=1, dens=1, visc=5, tau=time_step, fric=5)
system.actors.add(lbf)
system.thermostat.set_lb(kT=1)
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()
sys.stdout.write("\rSampling: %05i"%i)
sys.stdout.flush()
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 do 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.
In [ ]: