Bond atoms 13 14 missing on proc 0 with FENE potential in a polymer

Hi LAMPPS users,

I am simulating a polymer consisting of two types of beads with FENE potential as the bonded interaction and I am facing the error

    3561    22.705409     3.937074    26.642483    1.025e-05    2.6571199 -0.00011116411    11.367596 
    3562    22.700677    3.9103955    26.611072    1.025e-05    2.6391147 -0.000111255    11.367358 
    3563    22.696394    3.9114014    26.607796    1.025e-05    2.6397935 -0.00011112602    11.367103 
    3564    22.691122    3.9632069    26.654329    1.025e-05    2.6747569 -0.00011063481    11.366862 
    3565    22.685491     3.979644    26.665135    1.025e-05    2.6858503 -0.00011037473    11.366632 
    3566    22.680275    3.9782669    26.658542    1.025e-05    2.6849209 -0.00011025417    11.366414 
    3567    22.675294    3.9765629    26.651857    1.025e-05    2.6837708 -0.00011011049    11.366219 
    3568     22.67018    3.9600274    26.630207    1.025e-05    2.6726111 -0.00011007658    11.366011 
    3569    22.664926    3.9324554    26.597381    1.025e-05    2.6540028 -0.00011011125    11.365821 
    3570    22.659025    3.9038209    26.562846    1.025e-05    2.6346775 -0.00011016212    11.365639 
    3571    22.652774    3.9083975    26.561171    1.025e-05    2.6377662 -0.00010998461    11.365455 
    3572    22.646045    3.9102789    26.556324    1.025e-05     2.639036 -0.00010986949    11.365272 
    3573    22.639761    3.9093228    26.549084    1.025e-05    2.6383907 -0.00010976298    11.365089 
    3574    22.633287    3.9228475    26.556134    1.025e-05    2.6475185 -0.00010952026    11.364912 
    3575    22.627362    3.9335308    26.560893    1.025e-05    2.6547286 -0.00010930775    11.364731 
    3576    22.621503    3.9240301    26.545533    1.025e-05    2.6483166 -0.00010925498    11.364553 
    3577     22.61526    3.9226568    26.537917    1.025e-05    2.6473898 -0.00010912143    11.364365 
    3578 1.9171013e+32 1.8759915e+69 1.8759915e+69    1.025e-05 1.2661013e+69 1.2819275e+64    11.364181 
ERROR on proc 0: Bond atoms 13 14 missing on proc 0 at step 3579 (src/ntopo_bond_all.cpp:63)
Last command: run 10000000

However, for a input file with different number of beads like 61, I don’t get any error. I can see that the potential energy is blowing up but I dont know why is that happening.

Since I am a new user, I cant attach my input files but I am willing to share my lammps script, input file and the files for the non-bonded potentials.

Please help me resolve this issue.

Merry Christmas to all!!!

Thank you.

FENE bonds are guaranteed to blow up if they are somehow stretched past their constraint (this is a feature, not a bug), so I would immediately suspect that has happened. You may have a poor starting configuration that needs to be relaxed using minimisation or fix viscous. You may also need shorter integration timesteps. And you need to check if your force field is appropriate for your simulation.

Thank you for your response. I am running energy minimization. What is concerning is that it only happens when the chain length is short. Now that I can attach files, I am attaching the files here
dimer_41_2d_2LJ.input (5.5 KB)
LJ_sigma_1_epsilon_1.dat (888.7 KB)
simulate.lammps (1.5 KB)
SS_sigma_1_epsilon_1.dat (888.7 KB)

Also, could you give me some tips or point me to a resource where I can learn how to check if my force field is appropriate or not? I keep reading this everywhere but have failed to find out how to do the checks.

Thank you so much for taking time to answer.

A good force field correctly predicts enough things that you do know that you trust it when it predicts things that you don’t know.

For example, if you are studying polymers in solvents (I don’t study them), you would know what power laws are obeyed by their radius of gyration in good and poor solvents based on chain length (I don’t know those laws), and you will be able to test how well a force field reproduces those power laws (but I won’t). You can then argue based on those replications that the fine details of your simulations are also believable.

Notice I say “you” a lot there, and I hope the reason for that is clear from what I said. I don’t know what theory you study, what phenomena are more or less important to you, and therefore what force fields will be better or worse for your research. A theorist who studies polypeptides as theoretical polymers and a biochemist who studies polypeptides as drug candidates will have very different needs and use very different force fields even though they are studying “the exact same thing” (at least to a grant committee). For the same reason, even if I do look through your input and think it looks okay (and I might not have time to!), that does not guarantee that it looks okay. At the end of the day you will be writing your own paper and you will need to defend your own choices to the reviewers.

I have just briefly read through your attachments. You mentioned “obstacles” in one of your comments – this would have been very useful information to begin with.

I can see that you are trying to use energy minimisation to get you out of a possible unphysical initial condition – that may not work well with FENE bonds, especially if you are also using WCA beads*, which exhibit sharp repulsive forces.

If you know for sure that your initial configuration really is “close enough”, it might work for you to run minimisation with bond_style harmonic with the same equilibrium length, double-check that all bonds are shorter than the FENE maximum, save the configuration with write_data, and then run your simulation from that starting configuration with FENE bonds. But it may be simpler for you to just stick to generating initial conditions without overlaps, perhaps using a tool like Packmol.

*Why not use pair_style lj/cut? That would be far more readable both to helpers and to future you.

1 Like

Thank you so much for the advice and taking a look at the input. Since I am generating the initial conditions by self, I am running energy minimization. My initial condition in a straight polymer which surely is not close enough to the equilibrium structure and hence the EM.

Also the line about the obstacle, is an artifact from one my of previous scripts. I am not using any obstacles here. It is just a straight polymer in vacuum.

When I was first trying to use pair_style lj/cut the last pair_style was overwriting all the previous ones. However, just yesterday I figured out how to use two WCA and one LJ potential. This is what I am using now:

pair_style hybrid lj/cut 2.5 lj/cut 2.5 lj/cut 3.0
pair_coeff 1 1 lj/cut 1 1.0 1.0 1.122
pair_modify shift yes
pair_coeff 1 2 lj/cut 2 1.0 1.0 1.122
pair_modify shift yes
pair_coeff 2 2 lj/cut 3 ${eps} 1.0 3.0

That must be a consequence of a mistake on your side.

This is very much the same as:

pair_style lj/cut 3.0
pair_coeff 1 * 1.0 1.0 1.122
pair_coeff 2 2 ${eps} 1.0
pair_modify shift yes

The equilibrium ensemble of a WCA-FENE polymer (without external forces) will be overwhelmingly dictated by entropic, not energetic, considerations. This should tell you that energy minimisation should not make much difference to your run.

1 Like

A common approach for unoverlapping and equilibrating atoms in such a system would be the use of the soft pair style with fix adapt. Please see the micelle example in the LAMMPS examples folder for a demonstration of that strategy.

@akohlmey and @srtee, thank you so much for your advice and sorry for the late reply.

I will surely take a look at the micelle example.

Once again thank you so much.

Wishing you a happy new year in advance!