Monte Carlo simulation of a system of Lennard-Jones + FENE potential

Question

I want to generate two linear chains of 20 monomers each at some distance to each other. The following code generates a single chain. Could someone help me with how to generate the second chain? The two chains are fixed to a surface i.e the first monomer of the chain is fixed and the rest of the monomers move freely in x-y-z directions but the z component of the monomers should be positive.

Something like this:

import numpy as np
import numba as nb
#import pandas as pd

@nb.jit()
def gen_chain(N):

    x = np.zeros(N)
    y = np.zeros(N)
    z = np.linspace(0, (N)*0.9, num=N)

    return np.column_stack((x, y, z)), np.column_stack((x1, y1, z1))

    #coordinates = np.loadtxt('2GN_50_T_10.txt', skiprows=199950)
    #return coordinates

@nb.jit()
def lj(rij2):

    sig_by_r6 = np.power(sigma**2 / rij2, 3)
    sig_by_r12 = np.power(sigma**2 / rij2, 6)
    lje = 4 * epsilon * (sig_by_r12 - sig_by_r6)
    return lje

@nb.jit()
def fene(rij2):

    return (-0.5 * K * np.power(R, 2) * np.log(1 - ((np.sqrt(rij2) - r0) / R)**2))

@nb.jit()
def total_energy(coord):
    # Non-bonded energy.
    e_nb = 0.0
    for i in range(N):
        for j in range(i - 1):
            ri = coord[i]
            rj = coord[j]
            rij = ri - rj
            rij2 = np.dot(rij, rij)
            if (rij2 < rcutoff_sq):
                e_nb += lj(rij2)

    # Bonded FENE potential energy.
    e_bond = 0.0
    for i in range(1, N):
        ri = coord[i]
        rj = coord[i - 1] # Can be [i+1] ??
        rij = ri - rj
        rij2 = np.dot(rij, rij)
        e_bond += fene(rij2)
    return e_nb + e_bond

@nb.jit()
def move(coord):

    trial = np.ndarray.copy(coord)
    for i in range(1, N):
        while True:
            delta = (2 * np.random.rand(3) - 1) * max_delta

            trial[i] += delta
        #while True:
            if trial[i,2] > 0.0:
                break
            trial[i] -= delta



    return trial

@nb.jit()
def accept(delta_e):

    beta = 1.0 / T
    if delta_e < 0.0:
        return True
    random_number = np.random.rand(1)
    p_acc = np.exp(-beta * delta_e)
    if random_number < p_acc:
        return True
    return False

if __name__ == "__main__":

    # FENE potential parameters.
    K = 40.0
    R = 0.3
    r0 = 0.7
    # L-J potential parameters
    sigma = 0.5716
    epsilon = 1.0
    # MC parameters
    N = 20   # Numbers of monomers
    rcutoff = 2.5 * sigma
    rcutoff_sq = rcutoff * rcutoff
    max_delta = 0.01
    n_steps = 100000
    T = 10

    # MAIN PART OF THE CODE
    coord = gen_chain(N)
    energy_current = total_energy(coord)

    traj = open('2GN_20_T_10.xyz', 'w')
    traj_txt = open('2GN_20_T_10.txt', 'w')

    for step in range(n_steps):
        if step % 1000 == 0:
            traj.write(str(N) + '\n\n')
            for i in range(N):
                traj.write("C %10.5f %10.5f %10.5f\n" % (coord[i][0], coord[i][1], coord[i][2]))
                traj_txt.write("%10.5f %10.5f %10.5f\n" % (coord[i][0], coord[i][1], coord[i][2]))
            print(step, energy_current)
        coord_trial = move(coord)
        energy_trial = total_energy(coord_trial)
        delta_e = energy_trial - energy_current
        if accept(delta_e):
            coord = coord_trial
            energy_current = energy_trial


    traj.close()

I except the chain of particles to collapse into a globule.

Answer 1

Below you can find a modified version of your code. To make things simpler, in this version there is no box and no boundary conditions, just a chain in free space. The chain is initialized as a linear sequence of particles each distant 80% of R0 from the next, since R0 is the maximum length of the FENE bond. The code considers that particle i is bonded with i+1 and the bond is not broken. This code is just a proof of concept.

#!/usr/bin/python

import numpy as np

def gen_chain(N, R):
    x = np.linspace(0, (N-1)*R*0.8, num=N)
    y = np.zeros(N)
    z = np.zeros(N)
    return np.column_stack((x, y, z))

def lj(rij2):
    sig_by_r6 = np.power(sigma/rij2, 3)
    sig_by_r12 = np.power(sig_by_r6, 2)
    lje = 4.0 * epsilon * (sig_by_r12 - sig_by_r6)
    return lje

def fene(rij2):
    return (-0.5 * K * R0**2 * np.log(1-(rij2/R0**2)))

def total_energy(coord):
    # Non-bonded
    e_nb = 0
    for i in range(N):
        for j in range(i-1):
            ri = coord[i]
            rj = coord[j]
            rij = ri - rj
            rij2 = np.dot(rij, rij)
            if (rij2 < rcutoff):
                e_nb += lj(rij2)
    # Bonded
    e_bond = 0
    for i in range(1, N):
        ri = coord[i]
        rj = coord[i-1]
        rij = ri - rj
        rij2 = np.dot(rij, rij)
        e_bond += fene(rij2)
    return e_nb + e_bond

def move(coord):
    trial = np.ndarray.copy(coord)
    for i in range(N):
        delta = (2.0 * np.random.rand(3) - 1) * max_delta
        trial[i] += delta
    return trial

def accept(delta_e):
    beta = 1.0/T
    if delta_e <= 0.0:
        return True
    random_number = np.random.rand(1)
    p_acc = np.exp(-beta*delta_e)
    if random_number < p_acc:
        return True
    return False


if __name__ == "__main__":

    # FENE parameters
    K = 40
    R0 = 1.5

    # LJ parameters
    sigma = 1.0
    epsilon = 1.0

    # MC parameters
    N = 50 # number of particles
    rcutoff = 3.5
    max_delta = 0.01
    n_steps = 10000000
    T = 1.5

    coord = gen_chain(N, R0)
    energy_current = total_energy(coord)

    traj = open('traj.xyz', 'w') 

    for step in range(n_steps):
        if step % 1000 == 0:
            traj.write(str(N) + '\n\n')
            for i in range(N):
                traj.write("C %10.5f %10.5f %10.5f\n" % (coord[i][0], coord[i][1], coord[i][2]))
            print(step, energy_current)
        coord_trial = move(coord)
        energy_trial = total_energy(coord_trial)
        delta_e =  energy_trial - energy_current
        if accept(delta_e):
            coord = coord_trial
            energy_current = energy_trial

    traj.close()

The code prints the current configuration at each step, you can just load it up on VMD and see how it behaves. The bonds will not show correctly at first on VMD, you must use a bead representation for the particles and define the bonds manually or with a script within VMD. In any case, you don't need to see the bonds to notice that the chain does not collapse.

Please bear in mind that if you want to simulate a chain at a certain density, you need to be careful to generate the correct topology. I recommend the EMC package to efficiently generate polymers at the desired thermodynamic conditions. It is by no means a trivial problem, especially for larger chains.

By the way, your code had an error in the FENE energy evaluation. rij2 is already squared, you squared it again.

Below you can see how the total energy as a function of the number of steps behaves for T = 1.0, N = 20, rcutoff = 3.5, and also the last current configuration after 10 thousand steps.

And below for N = 50, T = 1.5, max_delta = 0.01, K = 40, R = 1.5, rcutoff = 3.5, and 10 million steps. This is the last current configuration.

The full "trajectory", which isn't really a trajectory since this is MC, you can find here (it's under 6 MB).

Answer 2

The FENE potential models bond interactions. What your code is saying is that all particles within the cutoff are bonded by FENE springs, and that the bonds are not fixed but rather defined by the cutoff. With a r_cutoff = 3.0, larger than equilibrium distance of the LJ well, you are essentially considering that each particle is bonded to potentially many others. You are treating the FENE potential as a non-bonded one.

For the bond interactions you should ignore the cutoff and only evaluate the energy for the actual pairs that are bonded according to your topology, which means that first you need to define a topology. I suggest generating a linear molecule of N atoms in a box big enough to contain the whole stretched molecule, and consider the i-th atom as bonded to the (i-1)-th atom, with i = 2, ..., N. In this way the topology is well defined and persistent. Then consider both interactions separately, non-bonded and bond, and add them at the end.

Something like this, in pseudo-code:

e_nb = 0
for particle i = 1 to N:
    for particle j = 1 to i-1:
        if (dist(i, j) < rcutoff):
            e_nb += lj(i, j)

e_bond = 0
for particle i = 2 to N:
    e_bond += fene(i, i-1)

e_tot = e_nb + e_bond

Answer 3

There is some problem with the logic of the MC you are implementing. To perform a MC you need to ATTEMPT a move, evaluate the energy of the new state and then accept/reject according to a random number. In your code there is not the slightest sign of the attempt to move a particle. You need to move one (or more of them), evaluate the energy, and then update your coordinates. By the way, I suppose this is not your entire code. There are many parameters that are not defined like the "k" and the "R0" in your fene potential