Fortran Simulation Features

Semiflexable polymers

The feature that wlcsim was originally designed around and for which it is named is it’s ability to quickly and accurately simulation polymers with stiffness ranging from quite rigid to very flexible.

For very stiff polymers, the usual wormlike chain is simulated.

For relatively more flexible polymers, the “stretchable, shearable” chain is used.

For VERY stretchable polymers, a purely Gaussian chain is used.

Of these, the middle strategy is where this code base stands out. The “stretchable, shearable wormlike chain” is described in “Discretizing elastic chains for coarse-grained polymer models” by E. F. Koslover and A. J. Spakowitz (2013). This model allows for quick but still accurate simulation of polymer statistics of semiflexible polymers of a wide range of stiffnesses (i.e. persistence lengths).

The effective potential for the stretchable, searable wormlike chain is given in terms of constants. The bending modulus, the stretch modulus, the shear modulus, the bend-shear coupling, and the ground state segment length. These constants are tabulated for different descritizations in input/dssWLCparams. The module MC_wlc handles interpolating from this table and calculating the potential between two beads. See MC_wlc.

Field theoretic interactions

Monte-Carlo simulations of solutions and melts of polymers are greatly accelerated by the use of a field-theoretic, binned-density approach for calculating interbead interactions. This approach is described in “Theoretically informed coarse grain simulations of polymeric systems.” by Pike et. al. (2009). With a further description of our implementation given in “Field-theoretic simulations of random copolymers with structural rigidity” by Mao et. al. (2017).

Under this approach the local volume fraction of each constituent is calculated for each of a grid of bins. The code for calculating the volume fractions and updating them after a Monte-Carlo move are

Subroutines and functions

subroutine mc_int_intitialize(wlc_p)

This subroutine calculates the field energy from scratch based on the positions and types of the beads as described in the hamiltonian. This subroutine also sets the value fractions PhiA and PhiB to 0 and sets x_chi. Volume fractions are linearly interpolated between bins.

Parameters:wlc_p [wlcsim_params,in]
Use :params (dp(), wlcsim_params()), energies (energyof(), maiersaupe_())
Called from:calculate_energies_from_scratch()
Call to:y2_calc(), interp(), hamiltonian(), calc_dphi()
subroutine mc_int_update(wlc_p)

This subroutine calculates the change in the field energy after a Monte Carlo move. Calculation only performed for bins containing moved (or changed) beads.

Parameters:wlc_p [wlcsim_params,in]
Use :params (wlcsim_params())
Called from:mcsim()
Call to:calc_dphi(), hamiltonian()
subroutine calc_dphi(wlc_p, ib)

This subroutine calculated the change in volume fraction in various bins assocated with the motion of bead IB. The old positon is taken from wlc_R while the new position is taken from wlc_RP.

Developer Note: This subroutine assumes you have set wlc_bin_in_list=FALSE some time before the start of the move.

Parameters:
  • wlc_p [wlcsim_params,in] :: data structure
  • ib [integer,in] :: Index of moved bead
Use :

params (dp(), wlcsim_params()), energies (energyof(), maiersaupe_())
Called from:

mc_int_intitialize(), mc_int_update()
Call to:

interp(), y2_calc()

The process of calculating the change in volume fractions after a move is often the computational bottle-neck in the simulation. Routines for doing this that are optimized for different move types are documented here Functions for updating densities.

The density is interpolated linearly between bins as described by Pike. This interpolation is performed in

Subroutines and functions

subroutine interp(wlc_p, rbin, ix, iy, iz, wx, wy, wz)

This program linearly interpolates a bead at RBin into 8 bins indexed by IX, IY, IZ with weights WX, WY, WZ. Interpolation is based on Pike (2009), “Theoretically informed coarse grain simulation of polymeric systems”.

For exampe a if a bead is closer to the center of IX(1) then IX(2) then WX(1) will be proportionally higher then WX(2). The total weight in any of the eight bins is given by WX*WY*WZ.

Parameters:
  • wlc_p [wlcsim_params,in] :: system varibles
  • rbin (3) [real,inout] :: position of bead
  • ix (2) [integer,out] :: Range of bins in x direction
  • iy (2) [integer,out] :: Range of bins in y direction
  • iz (2) [integer,out] :: Range of bins in z direction
  • wx (2) [real,out] :: WX(1) = (IX(2)*WLC_P__DBIN-RBin(1))/(IX(2)*WLC_P__DBIN-IX(1)*WLC_P__DBIN)
  • wy (2) [real,out] :: WY(1) = (IY(2)*WLC_P__DBIN-RBin(2))/(IY(2)*WLC_P__DBIN-IY(1)*WLC_P__DBIN)
  • wz (2) [real,out] :: Fractional contributions in z
Use :

params (dp(), wlcsim_params())
Called from:

mc_int_intitialize(), calc_dphi(), mc_int_chem(), mc_int_rep(), mc_int_super_rep(), mc_int_swap()

To run a simulation with the field theoretic interaction turned on set WLC_P__FIELD_INT_ON to true. You will also need to specify the relevant interaction parameters to your problem, for example WLC_P__CHI. The details of the interaction depend on problem (i.e. WLC_P__FIELDINTERACTIONTYPE). What each of these interactions are (or to add your own) see src/mc/mc_hamiltonain.f03.

Subroutines and functions

subroutine hamiltonian(wlc_p, initialize)

This subroutine calculates the field energy (Hamiltonian) based on volume fractions phi_A and phi_B. The expresion used to calculate the energy is different for different systems. The setting WLC_P__FIELDINTERACTIONTYPE determines which expression to use.

If initialize is true then the energy for all bins is calculated. Otherwise calculate only for bins specified by wlc_inDPHI.

Parameters:
  • wlc_p [wlcsim_params,inout] :: data
  • initialize [logical,in] :: Need to do all beads?
Use :

params (dp(), wlcsim_params()), energies (energyof(), chi_(), couple_(), kap_(), field_(), maiersaupe_())
Called from:

mc_int_intitialize(), mc_int_update(), mc_int_chem(), mc_int_rep(), mc_int_super_rep(), mc_int_swap()

If you edit this file be sure to edit both the initial calculation of volume fraction and change in volume fraction options.

As a final note, for liquid crystal systems the simulation has additional volume fraction fields that are turned on by WLC_P__CHI_L2_ABLE. There are five such fields and they correspond to the m values for sphericals harmonics of l=2 between -2 and 2.

Nucleosome geometry

For simulations of chromatin that require the geometry of nucleosomes to determine the entry exit angles of the DNA linkers the module nucleosome has the details.

Nuncleosome Module

Streamlined Energy Components

The are many different energy contributions. To keep track of these (or to add more) you should look to the energies module.

Energies Module