Hello:
I do not know how to solve this problem by myself.
I try to simulate a oil–water–surfactant system, but the water–oil beads radial pair distribution functions and interfacial tension are wrong. I cannot find the reason, the only thing I know when the system is without bonds the results are agree with the review.
The parameters and results are from 2013literature- Interfacial tension in oil–water–surfactant systems: On the role of intra-molecular forces on interfacial tension values using DPD simulations.
Can anyone please guide me where is the error, Thank you.
the Script file is
in.tensor (2.0 KB)
the data file is
modeldata (395.5 KB)
Please do not post the same question twice. If you want to change your question, you can always edit your post.
Sorry, I delete one.
When the simulation runs without a syntax error or the simulation crashing due to incorrect input data, then the problem is more likely a problem of not correctly representing the science.
Unfortunately, that is beyond the scope of this forum and a problem that you can best resolve discussion with people that know about your research. Those you won’t likely find online. It would rather be your adviser(s), colleague(s), or collaborator(s).
If you cannot reproduce the results from the reference publication, you need to more carefully compare the description of the simulation details with what you have input.
There is only one detail that springs to my mind based on your description, and that would be the choice of settings for the special_bonds command — LAMMPS documentation
By default non-bonded interactions with atoms up to 4 steps away within the bond topology are excluded (as is required by most molecular all-atom force fields). However, for DPD simulations, that is often not the case. I have seen different settings being used like including all non-bonded (e.g. special_bonds lj/coul 1.0 1.0 1.0) or excluding only 1-2 neighbors, i.e. the immediate bonds (e.g. special_bonds lj/coul 0.0 1.0 1.0) which also applies to FENE bonded model polymers.
Please also note, that these settings also affect the computation of the radial distribution functions g(r). If you need to include the excluded pairs in the rdf calculations, you may need to use 1.0e-100 instead of 0.0 to have those pairs included without contributing significantly to the interactions. Please also see the documentation for the compute rdf command — LAMMPS documentation
As mentioned above, this is not really a topic for a forum discussing the LAMMPS software itself. There is the “Science Talk” category here that is more appropriate, but likely talking to local people that know about your research and can help you compare the reference publication to your input and results is going to be the most productive. Except for the rare case that you come across somebody that has done research on the exact same subject using LAMMPS, you are not likely to get more specific input. To provide more detailed insight would require actually learning about your research and repeating at least part of your work and none of the people responding online have time to do that.
1 Like
Thank you very much for your suggestions. I use special_bonds lj 1.0 1.0 1.0, then the value of interfacial tension is changed. Could you tell what is the cpp file，it seems that there is no special_bonds.cpp
Of course it changes, this setting changes the model. Question is, was this the setting that was used in the paper you are trying to reproduce?
The cpp file for what? and to what purpose?
I can’t quite understand the content from special_bonds command — LAMMPS documentation. So I want to know the cpp file.
There is no single .cpp file. The command itself is executed from the Input class in input.cpp.
But the setting applies “everywhere” in pair styles and in the neighbor list code. There is the NEIGHMASK define in the lmptype.h header that is related to it and the sbmask() inline function in pair.h and the special list in the Atom class, in atom.h/atom.cpp.
If you need better understanding, you have to revert to a text book on classical MD and how molecular force fields are constructed. For consistency, non-bonded interactions (e.g. Lennard-Jones and Coulomb) are not computed for pairs of atoms that have an explicit bond, or angle or dihedral. That way computing the force constant for a harmonic bond from computing the vibrational spectrum with a quantum code can be directly determined without having to subtract the non-bonded forces first.
There also are four publications listed at the bottom of the documentation page that should consult.