There are three type of atoms A, B and C which are uniformly distributed in the system. I would like to compute force between A and B as a function of coordinates. I used “compute property/atom” to get the force of atoms A and used “ave/chunk” to get the the force distribution in three directions. However, the force is the total force of atoms A which includes the unwanted component between A and C. I can used the “compte group/group” to get the force btween A and B, but I am not able to use “ave/chunck” to get the distribution as coordinate.
In post-processing this would be relatively easy – just “turn off” the interactions of atoms C (set coefficients and charges to zero or use neigh_modify exclude or something else) and rerun the trajectory, and you will have per-atom forces that are purely A-B interactions.
During the simulation itself? Doesn’t seem possible. But you could use compute coord/atom to calculate the average number of B atoms around each A atom and their average distances, and even try to chunk that, and then work out the average forces as a function of coordination-density.
The rerun option seems a good choice for me.
To get the force distribution and trajectory, I used below command in the running script
compute friction sodium property/atom fx fy fz
compute friction_reservoir1 sodium chunk/atom bin/3d x lower 3 y lower 3 z -42 3
fix profile_Fxyz sodium ave/chunk 10 100000 1000000 friction_reservoir1 c_friction[*] ave running file reservoir_Fxyz.txt
dump rerunDump all custom 1000 rerun.lammpstrj id x y z vx vy vz fx fy fz
In the rerun script, I set the unwanted coefficient to zero and keep all other commands same, including the “compute property/atom”, “compute chunk/atom” and “fix ave/chunk”. The simulation trajectory read by “rerun rerun.lammpstrj dump x y z vx vy vz box no format native”. However, I found there was no output in the reservoir_Fxyz.txt where force distribution btween A and B is supposed to be found.
If nobody will want to read it, then nobody will be helping you, either. So you have to simplify it until it is readable and so anybody wants to help you. In the current situation, it is impossible.
Thus narrow down your input to the essential bits, take out everything that is not relevant to the issue you are experiencing and use a small test system. That will also make it easier for yourself to debug it and perhaps you can understand the problem and figure out the solution by yourself.
How should one provide assistance, if there is not way to confirm your issue? The only option is to wildly guess, which can then also be wildly wrong.
As @akohlmey said, you really should share your script here. And if you are not comfortable with that you should pause and spend a day (or a week) organising your script so you think it would make sense to someone else. Not only is that important for you to get help on this forum, but it is important for two other reasons.
Firstly, I guess at some point you want to write and publish this work as a paper. When you do, some reviewers (unfortunately not enough) will want to look over your scripts and confirm that they actually do what you said they do. If you say at that point “sorry but I can’t show you my scripts”, they have every right to reject your paper, because science should be reproducible. Even if a reviewer does not ask for your scripts, future collaborators may be interested in your work and want to know how to replicate it. How will you help them replicate your work if you can’t share your scripts? Tidying up your scripts is a simple and necessary step to help other researchers recognise and build on your work.
Secondly, if your script is very complex, how do you know it doesn’t have other bugs in it? After all you already know that some lines in your script are not doing what you expect. How do you know that all the other lines in your script are correct? If you take a step back now and double-check everything, in the process of tidying up your script, you may have the chance to spot and correct other errors. On the other hand, if you have written the paper and drawn the graphs and then only clean up your script, you may find that you have done three months of simulations that were wrong to begin with. That is not a good outcome.
I have simplified the input and attached them here. I would be appreciate if someone could find the reason why there is no output for the rerun. 50000atm_TIP3P.lmp (8.0 KB) 50000atm_TIP3P_rerun.lmp (8.0 KB)
If that is simplified, I really don’t want to see the “complicated” version. This is still very convoluted and complex and not at all a minimal input version that just reproduces the specific problem for just a minimal and simple input deck. It is not even possible to just rerun your input, since the data file is missing.
Sorry for the complex input, however those command are necessary to run my model. please find new the attached input and data file, so someone can reproduce the my problem that there is no output in the reservoir_Fxyz.txt after rerun. 50000atm_TIP3P.lmp (6.3 KB) 50000atm_TIP3P_rerun.lmp (6.3 KB) eqm_eqm_4.0ns.data (3.6 MB)
Using comments for version control (keeping track of past versions of a file) is fine for very quick experiments, but makes your file incredibly confusing after even a few rounds. (For example, you will not know which combinations of commented commands were used in previous runs.)
To really do better you must learn something like git. But the most basic level of version control is: every time you move on to a new version of a file, save the old one with (ISO formatted) date in its name (so 50000atm_TIP3P_20102023.lammps). Then your new file will only hold necessary comments.
fix ave/chunk knows how to take fx, fy and fz as variables. You do not need separate compute property/atom commands to track them.
many names can be shortened, such as c for carbon and na for sodium. (Use lowercase letters for those group names, since you use uppercase letters like Na as your element labels later.)
similarly something like
region upper_cnt ...
group upper_cnt ...
will work instead of what you have. LAMMPS knows when (for example) upper_cnt makes sense as an input group or as an input region, you usually shouldn’t have to track that manually.
if pair_coeffs and so on are specified in a data file, they don’t need to be re-specified in the input.
your fix ave/chunk timestep settings may be incompatible with the rerun trajectory length.
finally, I really, really hope you have some method of externally validating your results, because you likely have a big thermostat problem. This seems to be rarely appreciated (probably because it’s rarely relevant), but using compute temp/partial with rigid molecules that are also partially ordered will often measure incorrect temperatures, because the rigid molecules’ degrees of freedom will either over-correlate or under-correlate with individual particles’ Cartesian axes. (As a first-level example, in a system of nematically-aligned rigid dumbbells, the kinetic energy will not be equipartitioned between the x, y and z-kinetic energies of the dumbbell ends.)
You didn’t follow what I was asking. I asked to narrow it down to the essential bits that reproduce the problem (lack of output). That does not automatically mean that you have to simulate the entire system with all bells and whistles. In many cases, you could even substitute a much simpler input, e.g. the LAMMPS developers often try to reproduce issues using the melt example input or the peptide example input or similar. Those are small, run fast, and one can often easily make a few small modifications to recreate the problematic command, e.g. average over a different property.
More generally speaking, I share the concerns of @srtee, that when you do something this complex and convoluted, you will have a hard time explaining and justifying your simulation setup. I have serious doubts that it has to be this complex and given your recent questions t the forum, I have even bigger doubts that you have verified and validated each bit to give meaningful results, and thus putting your entire simulation setup at risk of being a complete waste of your time. As was already suggested, I also recommend to take a step back, look at all you have got right now, look through all the comments and suggestions you have received and then start all over. Build your system in steps and check at each step that you are getting meaningful results. Start with a simple setup and verify the analysis, then make your simulation more complex and confirm that what you get as results is still meaningful. You need very much try to keep things simple, since most things are simple. Going for complexity is often not improving the answer you can get, it just makes it harder to understand and explain to others.
Maybe we should take a page from GROMACS and have LAMMPS remind us at the end of every simulation:
Since all models are wrong the scientist cannot obtain
a “correct” one by excessive elaboration. … Just as the ability to devise simple but evocative models is the signature of the great scientist so overelaboration and overparameterization is often the mark of mediocrity. Since all models are wrong the scientist must be alert to what is importantly wrong. It is inappropriate to be concerned about mice when there are tigers abroad. – George Box, Science and Statistics, JASA 1976.
