Why am I obtaining a different polymer behavior when I turn off the WCA component of FENEWCA and use LJ instead?

[aldorvasq]

I am trying to replicate the results of a paper where a coarse-grained (CG) model of amphiphilic block copolymers in a water-oil biphasic system using Langevin Dynamics is used (https://doi.org/10.1063/5.0189156).
In that article the authors mention the use of the FENE potential for bonded interactions and LJ potential for all CG beads. To recreate that, I set all the values (epsilon, sigma and delta) that correspond to the WCA component in the hoomd.md.bond.FENEWCA() class to zero. As I understood from the documentation, neighbor lists include all beads unless explicitly told otherwise, so I would expect the 1-2 interactions to use the md.pair.LJ() interaction in addition to the FENE potential.
What I am observing in my simulations is that when I turn off the WCA values of the FENEWCA, all the polymer chains bundle up into solid-like structures, when the expected behavior should be more of a fully dissolved polymer chain. I tried turning on the WCA values to epsilon=1.0, sigma=1.0, delta=0.0 and excluding 1-2 interactions, and then I was able to obtain the expected behavior, but I was not able to reproduce the interfacial tension values the paper reported making me think that the WCA part was indeed not included in the bonded interactions.

Here is a snippet of the script that I used where I set the WCA values to zero and includes 1-2 interactions. The rest of the script is is just where I save the trajectories, loggable quantities and I set the run() line.

import hoomd
import datetime
import numpy as np

processor = hoomd.device.GPU()
sim = hoomd.Simulation(device=processor, seed=2467)
sim.create_state_from_gsd(filename='initial_configuration.gsd')
integrator = hoomd.md.Integrator(dt=0.01)
fene = hoomd.md.bond.FENEWCA()
fene.params['polymer'] = dict(k=30.0, r0=1.5, epsilon=0, sigma=0, delta=0)
integrator.forces.append(fene)
neighbor_list = hoomd.md.nlist.Cell(buffer=0.4, rebuild_check_delay=1)

epsilon = 1  # reduced units
sigma = 1  # reduced units
alpha = 0.5
beta = 0.5
lj = hoomd.md.pair.LJ(neighbor_list, default_r_cut=2.5*sigma, mode='shift')

lj.params[('A', 'A')] = dict(epsilon=epsilon, sigma=sigma)
lj.params[('B', 'B')] = dict(epsilon=epsilon, sigma=sigma)
lj.params[('S', 'S')] = dict(epsilon=epsilon, sigma=sigma)
lj.params[('W', 'W')] = dict(epsilon=epsilon, sigma=sigma)

lj.params[('A', 'B')] = dict(epsilon=beta * epsilon, sigma=sigma)
lj.params[('A', 'S')] = dict(epsilon=epsilon, sigma=sigma)
lj.params[('A', 'W')] = dict(epsilon=np.sqrt(alpha * beta) * epsilon, sigma=sigma)
lj.params[('B', 'S')] = dict(epsilon=np.sqrt(alpha * beta) * epsilon, sigma=sigma)
lj.params[('B', 'W')] = dict(epsilon=epsilon, sigma=sigma)
lj.params[('S', 'W')] = dict(epsilon=alpha * epsilon, sigma=sigma)

integrator.forces.append(lj)

ensemble = hoomd.md.methods.Langevin(filter=hoomd.filter.All(), kT=1.0, default_gamma=1.0)
integrator.methods.append(ensemble)

sim.operations.integrator = integrator

At this point I am unsure what could be the issue and I would appreciate any feedback on this problem that I am encountering.
Additionally, I used LAMMPS just as the authors did, and I was able to recreate the interfacial tension values and I did not see any solid-like structures.

I will attach pictures of the solid-like structures and the correct dissolved structures too that I obtained with LAMMPS.
I am using hoomd 4.4.1

[WeikangXian]

Hello,

It will be easier for the other to comment if the LAMMPS script is also provided. I assume that you want to reproduce the results from the paper via HOOMD.

Are you expecting the provided HOOMD snip to produce the expected results? In your snip, WCA is turned off. But LJs are included. Such a setting makes the difference in the cutoff. The cutoff for WCA is 2^(1/6) but it is 2.5, as you assigned, for LJs.

Hope it helps.

Best,
WX

[chrisjonesBSU]

Just to clarify one thing.

As I understood from the documentation, neighbor lists include all beads unless explicitly told otherwise, so I would expect the 1-2 interactions to use the md.pair.LJ() interaction in addition to the FENE potential.

The default behavior for the neighbor lists is to exclude bonded particle pairs, as you can see in the docs

For example, with the cell neighborlist:

Class hoomd.md.nlist.Cell(buffer, exclusions=('bond',), rebuild_check_delay=1, check_dist=True, deterministic=False, mesh=None, default_r_cut=0.0)

@WeikangXian makes a good point about possible inconsistencies in the r-cut value when deciding whether to use WCA or 1-2 LJ pair force with the FENE bonds. Another difference is that the WCA term in FENEWCA shifts the potential up by $\epsilon$ whereas the LJ12-6 potential doesn't, in case that makes a difference in how you are calculating surface tension.

[aldorvasq]

Thank you @WeikangXian and @chrisjonesBSU for your responses.

I was able to reproduce the expected behavior and interfacial tension values once I set exclusions=[] in the Cell() class. I thought that, since the hoomd.md.nlist.NeighborList base class mentions in the documentation that ''Neighbor lists nominally include all particles within the chosen cutoff distances.", it also meant the same for the Cell() class, but I was mistaken.

As for the WCA part, setting mode='shift' in hoomd.md.pair.LJ() should also shift the potential up by the same amount that FENEWCA does, right? I looked at the documentation for hoomd.md.pair.LJ() but I could not find where it explicitly says that mode='shift' shifts by $\epsilon$, although I would assume so.

Thank you once again,
AV

[chrisjonesBSU]

I didn't notice you were using shift in hoomd.md.pair.LJ(). I believe the value that shift uses is whatever is needed to shift the potential at r_cut to zero, so it depends on what you're using for r_cut