I am writing a code to simulate a very simple coarse-grained polymer with fene bond assuming there are no angles or dihedrals and finally need to calculate the MSD of the polymer. I tried the code below but the beads become detached and segregated during the simulation (maybe due to the dominance of repulsive forces).
Can anyone help me to understand why this happens?
Here are my code and initial configuration file of a polymer of 50 beads.
Thanks,
Mobi
###############################################
# LAMMPS script for a polymer
###############################################
###
# Box and units (use LJ units and periodic boundaries)
###
units lj # use lennard-jones (i.e. dimensionless) units
atom_style bond # atoms with bonds
boundary p p p # all boundaries are periodic
###
# Pair interactions require lists of neighbours to be calculated
###
neighbor 1.9 bin
neigh_modify every 1 delay 1 check yes
###
# READ "start" data file
###
read_data initial_config50.txt
###
# Reset timestep
###
reset_timestep 0
###
# Define groups
###
group all type 1 #(atom type 1 is group 'all')
variable seed equal 54654652 # a seed for the thermostat
# Equilibration (Langevin dynamics at 300 K)
velocity all create 1.0 69878
fix 1 all nve/limit 0.05
fix 2 all langevin 1.0 1.0 1.0 ${seed}
thermo_style custom step temp epair emol press vol etotal
thermo 10000
timestep 1
dump 1 all xyz 1000 Dump50.xyz
###
# Potentials
## Between bonded atoms
bond_style fene
special_bonds fene #<=== I M P O R T A N T prevents LJ from being counted twice
# For style FENE, specify:
# * bond type
# * K (energy/distance^2)
# * R0 (distance)
# * epsilon
# * sigma
bond_coeff 1 30.0 1.5 1.0 1.0
## Between non-bonded atoms
pair_style lj/cut 1.12246152962189
pair_modify shift yes # option to ensure energy is calculated corectly
# pair_coeff for LJ, specify 4:
# * atom type interacting with
# * atom type
# * energy
# * mean diameter of the two atom types
# * cutoff
pair_coeff 1 1 1.0 1.0 1.12246152962189
run 50000
#### write a final restart file
write_restart final.restart
Please have a look at the warnings that LAMMPS is printing.
Your first “run” command is before you define your interactions with pair/bond style and coeffs. So you are effectively simulating something that is more like an ideal gas and not a polymer.
This is the “LAMMPS Beginners” category, so there is more leeway to ask basic and simple questions than elsewhere. However, you should keep in mind that you should by preference ask questions that are related to technical aspects and LAMMPS specific, for more general questions, the first source of advice should always be your adviser and/or your colleagues.
Those are OK, if they only happen at the beginning and and not later in the simulation. They are usually a consequence of your initial geometry not being consistent with the force field settings and then - while the system is equilibrating - you may have (temporarily) a lot of kinetic energy locally that will eventually be dissipated.
Have a look at the evolution of the potential and kinetic energy in the thermo output.
There are multiple ways to mitigate this:
create a “better” initial geometry that is closer to equilibrium
start the simulation with a minimization, or - in more severe cases - do multiple minimizations (sometimes in this context referred to as a “quench”) interleaved with shorter MD runs until the potential energy of the system has stabilized
use (temporarily) a shorter time constant for the langevin thermostat (so it dissipates more energy due to larger friction and stronger randomization) and fix nve/limit instead of fix nve to reduce the maximum displacement of atoms
Please note that your system is very small, so it will have larger fluctuations than a larger system. But it is quite common to start this way and then enlarge the system with the replicate command after initial equilibration and save some time.
Those are followed by an error “ERROR on proc 0: Bad FENE bond (…/bond_fene.cpp:90)” and the code is not running now.
Thanks for your advice, but I am self-learning!
And this is nice that in this digital world, one can find a colleague or advisor in a far distance, isn’t it?
Sorry, but I disagree with that assessment. It won’t work for MD so well for multiple reasons:
answering questions in public forums related to software is usually left to the developers of that software and now there is the problem that beginners mix up questions about the methodology with questions about the software and thus people like me are expected to become adviser to lots of people. That cannot work unless I would stop developing and even then the numbers would be stacked against me.
Since advising is time consuming and research is extremely competitive, it leads to very few people that should be sharing their knowledge (e.g. all of those that had been helped by a developer in a similar situation) to do so. On one hand this is understandable behavior, but on the other hand, this is rather unsocial, and also a missed opportunity: there is no better way to improve one’s understanding of something than helping somebody else that is a beginner. That is why I am still responding since it keeps my understanding “alive” and sometimes even improves it, even though I am not doing research on my own anymore (no time left for it).
Doing good(!) MD simulations is much more a craft than a science. There are lots of best practices and conventions and commonly used practices that are far too many to explain them over the internet or write them into text books etc. In many cases those are specific of the area of research and not at all specific to the software in use. Add to that, that there is a very large pool of knowledge in the published literature, too large to digest and process this on your own (MD simulations on computers go back over 50 years!), but also without that kind of knowledge you are bound to repeat all the mistakes and errors that in some cases took people years or even decades to figure out. Thus the most effective way to enter the MD simulation field - like in other crafts - is through some kind of apprenticeship. In the pre-internet times, that was the common procedure: to learn how to do MD simulations one would join a group that had the tools and residual knowledge (either as a visitor or as a student/postdoc) do some collaborative research and then move on (and perhaps learn additional skills from another research groups like a “journeyman” in the crafts).
Bottom line, you will be able to do some level of calculations with the guidance you can (technically) get through forums and emails, and they may be acceptable for publication in less discriminating journals, but without proper guidance and tutoring, you are doomed to become at best some kind of “expert beginner” because you are cut off from the pool of residual knowledge that is vital for doing competitive research in this field. You have to realize that a simulation that runs to the end without warnings or error may still be completely bogus and useless. Comparing to the “craftsman” example, this is like being, say, a plumber that knows how to use all the required tools and still creating a mess by using the wrong materials, making the wrong connections and laying pipe where they shouldn’t be (e.g. where they can easily freeze in winter).
Thanks for taking the time to answer.
All you said makes sense and I guess there are some misunderstandings.
I didn’t mean that I intend to learn LAMMPS personally forever. I am trying to do some hands-on MD simulation projects to assess both my ability and interest, so that I can decide better whether to start a post-doc on simulation or not.
I try forums so that I may, by chance, find people haveing time to give a quick guide, or having the same problems as me, along with other ways available.
With my 25+ years experience in the field of MD simulations with 20+ years of responding to questions on multiple mailing lists and forums, I still think that you may get the wrong impression and learn bad habits without personal tutoring, and that may taint your work for a long time regardless of whether you will join a group later on or not. The lessons learned at the beginning have the most impact, so the importance of learning the basics “the right way™” cannot be underestimated. Try googling the term “expert beginner”. There are some very good essays on that topic online that somebody pointed out to me a long time ago, and that describe the dangers of that approach. Those essays are in the context of software development, but they translate well to research, IMO.
I should add to that, however, that my success rate in convincing people to change their ways is extremely low. There seems to be a level of optimism in people believing that my predictions don’t apply to them, that is beyond my ability to convey my observations. In fact, I would say that - if anything - it is getting even lower as software becomes more accessible (and thus people with less domain specific knowledge use it), has more features (and thus becomes more attractive), and people are migrating from maintaining small, domain specific in-house software to adopting and extending large package software (and thus saving themselves a lot of time maintaining the in-house software and implementing non-domain specific features that other software already has).
Thanks for mentioning it. It was interesting. It seems the world is filled up with the “expert beginners”.
Hahaha.
I cannot judge other people you interacted with, but for me, it was very valuable to be aware of such a hazard and take your advice into account.
