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.
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.
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.
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.
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:
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.
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.