I am trying to build CG force filed. The parameters of bond potential and angle potential energy have been fitted and reproduced in the simulation, but there are some questions about whether the pair potential has been affected by the potential related to the bond when solving the non-bond interaction?
1.I have tried some 12-6LJ but can’t reproduce the number and hight of peek in RDFs, I wonder if Pair potential is also affected by bond, and RDF cannot be calculated without setting pair potential, I am now using the windows version of lammps and it seems that I can’t view the source code, can you tell me how to get the source code?
2. If I use tabular potentials, I don’t know how to deal with the column of forces in the table of potential functions, and the forces obtained by the central difference method are not continuous or the curve is not smooth enough, will this produce errors? What is the range of potential energy I can set for the part of RDF=0 that appears in an all-atom simulation?
3. Finally, there is the question about the RDF algorithm. RDF should finally converge to 1, but in my many attempts, I have used 1fs as a timestep in time and simulated 100000000 steps, and sometimes there is a problem that RDF does not converge to 1, is this the reason why the table potential is too large near the cutoff?
looking forward to your reply
Many of the questions you are asking are about the general science of specifying the pair potential, fitting the pair potential to the target RDF, and in particular how to account for bond forces vs pair forces. These are general questions not to do with LAMMPS specifically, and unless a domain expert with specific experience in these areas comes upon this post and is willing to help (out of their own good will!) you may have to look elsewhere for guidance.
Having said that, with this part:
(I assume you mean that the RDF does not become 1 at a large enough radius) you may need to make sure the RDF compute has a large enough neighbor distance to work with. Any potential will induce ordering past the pair cutoff, and the compute needs to “scan” a large enough radius to get the convergence you are looking for.
Thank you very much for replying to me, what I want to know is whether the pair potential in lammps is affected by bond potential, whether the potential energy part related to bond is ignored in order to get the right pair potential, and whether the calculation of RDF in lammps ignores the atomic pairs that form bonds
See special_bonds command — LAMMPS documentation (which is also mentioned in the documentation for compute rdf).
The LAMMPS source code is available on the LAMMPS homepage as tarball or on GitHub with multiple download options.
This is really not a LAMMPS question but a question about the method. You need to consult with an expert willing to tutor you on building coarse grain methods. Mind you, g(r) inversion is the oldest and least sophisticated approach to coarse graining with known deficiencies. LAMMPS does support other approaches like force matching by the Voth group in the MSCG package and the SPICA force field in the CG-SPICA package.
There is not enough information here. Whether the g(r) converges to 1.0 or not largely depends on your system geometry and how you compute and average the g(r). You don’t provide any information about this, so it is not possible to give any specific comment or recommendation.
As others said, there is very little we can do to help you here. However as I can relate in the “been there, done that” kind of people here is my two cents:
- It is not said that a LJ potential can reproduce CG structure, also you do not say which system you want to make a CG model of. Large molecule? Polymers? Atomistic liquid? Those are different things that people tried to make CG models for. But they behave differently, and worse, the structure you get can change depending on the coarse-graining level. Lots of tricky stuff hide in there. Building a CG model is a lot of work, and a pain.
- The tabulated forces in LAMMPS are interpolated using linear or spline functions. I suggest you play around with test cases of small number of points and use the
pair_writecommand to see what forces is actually computed. See what happens when you set values near r=0 and how they are interpolated (LAMMPS should crash if particle overlap when dr is exactly 0 or below the minimal value anyway). - As stated the convergence is true for liquid-like structures over large enough range. However due to small fluctuations, the RDF value can oscillate around 1 over long range. What you can do is set this values to 1 for computation purpose until they go far enough from unity as r reduces, or at least beyond a cut-off value. Discuss this with your colleagues, make tests, see what it gives.
- I think the misunderstanding you have is on how the
compute rdffunction works. The computes can use neighbour lists of the atoms for efficiency, that means that, in case of bonded interaction, thespecial_bondsvalue can have an effect: if their value is 0, atoms can be excluded. This is discussed in the manual page of the RDF compute command. You can set the values to tiny values like 10^{-15} to still build the neighbour lists without the interaction between bonded particles being meaningful. The important thing is to know if you want pair interaction plus bonded interaction between your CG beads or only bonded interaction and to be consistent with your choice. But this is model dependent, and again, we don’t know much about what you are working on.
For what it’s worth, if you really want to do IBI to compute CG potential, there is this paper that we published not so long ago from my previous position. The results are clearly not scientifically mind-blowing but you’ll find literature discussion and a detailed procedure that might fuel discussions with your advisor/colleagues. If useful, please cite yadi yadi yada, but as @akohlmey said, this method has… well… known deficiencies. 
A couple more comments on this:
- pair style table will refuse to use tables that start at r=0
- there is some hopefully useful python code for creating tables (including numerical differentiation) in the tools/tabulate folder (added recently). However for doing numerical differentiation of coarse data (as is the case here), you cannot just take the data points as is. But rather need to use some interpolation. I have found that a Savitzky–Golay filter - Wikipedia is very attractive in many cases and particularly in this case, since you can apply Horner’s rule to the interpolation polynomial and thus obtain the numerical derivative from the smoothed interpolation directly. I have used this in the past with great success for computing self-diffusion from MSD (with numerical differentiation, you don’t need to do a linear regression, but you can directly pick an interval of the differentiated data and average over it. This can be very helpful in determining which part of the data is sufficiently well converged).
- rather than just setting the g(r) to 1.0 from a given point, I would recommend multiplying with a switching function (e.g. tanh()). This is a common technique also used in other places like NMR spectroscopy to create smoother data or suppress artifacts from noise. One of the examples in the tools/tabulate folder uses a tanh() switching function to smoothly combine a LJ attraction with a Morse repulsion so that not only a smooth energy function but also a smooth force function is obtained.
I forgot to mention how SG filters also saved the day more than once in my case when doing this kind of stuff.
Perhaps we should add an example to the tools/tabulate folder that generates the forces from tabulated energy data? Unfortunately, I am a bit busy with other things™ right now…
I’m sorry I didn’t describe my work clearly, I’m working on the development of a coarse-grained force field for a polymer chain system（poly BAMO）, unfortunately I’m not a chemistry or polymer physics major before, and there may be some misunderstandings about this kind of research. According to part of the literature I read, Utotal=Ubond+Unonbond, through the IBI method Unonbond=-kbtlng(r), and requires further iteration. On the one hand, I am not very clear that for coarse-grained polymer chains, g(r) here should use the RDF of beads between different polymer chains to completely eliminate the effect of Ubond, or only eliminate the RDF of beads that are connected to the adjacent bonds but within the same chain.On the other hand, I don’t have a good understanding of the RDF calculated by the lammps command, and here RDF excludes the role of Ubond. Anyway, thank you for some suggestions, and I will read your article well.
Thank you very much for your reply
I’m looking forward to the implementation of this feature, and the extra computing force is a bit annoying.
And my own differential method always fails to make potential energy and force continuous, I will try the SC filtering you mentioned. Thank you very much for the reminder.
Another question is how to view the source code of a particular command? It doesn’t seem to be categorized by command in the github you provided earlier, and I’ve been looking for a long time to find how the code ‘compute myRDF1 all rdf 250 2 2’ that computes RDF in lammps works.
You have to improve your looking skills. There is a rather simple and straightforward system in place. If the name of the command is “compute rdf”, then there is a file compute_rdf.cpp and a file compute_rdf.h that define the specific compute class instance and its implementation. Since LAMMPS uses C++ and inheritance, there is common functionality in the corresponding base class (Compute) and the files compute.cpp and compute.h and sometimes the individual commands are a modification of another command and thus they are not derived from the lowest base class, but the class of some other command.
If you want to understand the LAMMPS source code in general you need to spend time reading: 3. Modifying & extending LAMMPS — LAMMPS documentation and 4. Information for Developers — LAMMPS documentation