I am trying to simulate a bead-and-spring polymer interacting with bivalent particles and free particles, following the procedure described in this article. Here’s a brief description:
Bivalent particles: Each bivalent particle consists of a core with two patches (using “bond_style hybrid fene harmonic”).
Interactions:
Excluded-volume interactions among monomers (Type 1), cores (Type 3), and free particles (Type 4).
Patches (Type 2) interact with monomers (attraction) and themselves (excluded-volume) but not with cores or free particles.
I set the LJ coefficients as:
...
#### Groups
group mon type 1
group patch type 2
group core type 3
group prot_free type 4
pair_style lj/cut 1.122462
pair_modify shift yes
pair_coeff 1 1 1 1 1.122462
pair_coeff 1 2 30 0.535 0.685 # attraction between patches and monomers
pair_coeff 1 3 1 1.122462
pair_coeff 1 4 1 1.122462
pair_coeff 2 2 1 0.178 0.2
#pair_coeff 2 3 0 0.535 0.685 # comment out: no pair interaction
#pair_coeff 2 4 0 0.535 0.685 # comment out: no pair interaction
pair_coeff 3 3 1 1 1.122462
pair_coeff 3 4 1 1 1.122462
pair_coeff 4 4 1 1 1.122462
...
Given that the “pair_coeff” documentation says, “…if coefficients for type pairs with I != J are not set explicitly by a pair_coeff command, the values are inferred from the I,I and J,J settings by mixing rules”. and the “pair style lj/cut” documentation indicates the default mix value as geometric and refers to the “pair_modify” documentation for further fine-tuning, have I set up my coefficients correctly to achieve my intended interactions?
If not, I am considering two solutions to make atom types non-interacting:
Approach one: Uncommenting the following lines above so patches interact with cores and free particles with a “zero” interaction strength:
pair_coeff 2 3 0 0.535 0.685 # pair interaction with zero strength
pair_coeff 2 4 0 0.535 0.685 # pair interaction with zero strength
Approach 2: Using the exclude setting in neigh_modify
neigh_modify exclude group patch core
neigh_modify exclude group patch prot_free
Either of your “approaches” should work. From what I can see, approach 2 may be faster if it skips computation of a large number of interactions (instead of just a few; this depends on the concentration and nature of your solution).
However, some compute calculations may not work during the simulation since they also rely on the neighbour list (for example, the patch-core RDF can’t be calculated with compute rdf during the simulation if there are no patch-core pairs in the neighbour list). You can do calculations either in post-processing or by using approach 1.
In either case, the LAMMPS log file will contain a line saying (something like) 12 of 21 pair coefficients generated by mixing rules if pair coefficients have been mixed. You can also explicitly check all pair coefficients by using write_data with the pair ij keyword-value pair.
Thanks for your response. The log file has no message saying the undetermined coefficients are set by the default “geometric” mixing rule. However, the coefficients generated by the geometric mixing rule are written in the data file. Does this mean the coefficients were generated at the beginning of the simulation and used throughout it? Or does it mean the write_data just generated them at the end of the simulation?
In any case, write_data does not recompute pair coefficients, so it should reflect the system as of the last time pair coefficients were calculated in the script.
I double-checked the log file and compared it to the page provided. My log file contains the message, “WARNING: Too many warnings : 171 vs 100. All future warnings will be suppressed.” This likely explains why I don’t see the “Mixing forced for lj coefficients” warning in my log file.
If your input creates this many warnings, there is something not quite right and those warnings should be addressed. Since you didn’t share your complete input, it is difficult to make any specific statements and give suggestions.
I employ the Kremer-Grest model (FENE + LJ) to examine a lengthy polymer. Initially, the polymer adopts a helical conformation, leading to warnings like the one below during the first several hundred integration steps of the equilibration phase:
WARNING: FENE bond too long: 74 89 90 1.426624559474392 (src/MOLECULE/bond_fene.cpp:83)
While I attempted to use a larger cutoff distance to circumvent this warning, it led to a slowdown in the simulations.
Should I increase the warning message limit to view messages generated later in the log files?