But some results from the simulations in NPT ensemble using my and hybrid pair styles are different.
Ecouple increases linearly over time during the simulation for my pair style and fluctuates around a certain value during the simulation for the hybrid pair style.
Density, potential enegy and enthalpy from the simulations using my and hybrid pair styles differ but each fluctuates around a certain value during the simulations.
Temperature, kinetic energy and pressure from the simulations using my and hybrid pair styles don’t differ much and fluctuate around a certain value during the simulations.
I would appreciate if someone could explain the differences and suggest some ways to fix my issue. I supply more details and materials about force fields and simulations below.
I used the force field code pair_born_coul_long.cpp as a starting point for my code and edited the Van der Waals force and energy expressions and coefficients in the corresponding sections. Then I moved the header and body files in the src/KSPACE directory and compiled LAMMPS (29 Aug 2024 - Update 1) with EXTRA-PAIR and KSPACE packages.
I tested my force field on Al2O3 melt, the initial configuration containing 1500 atoms was obtained by Packmol package. All simulations were carried out in NPT ensemble at 3400 K.
At first, I compared the pairwise van der Waals energy and force as a function of the interatomic distance for my and hybrid pair styles and that part of the equation from the article. I used the pair_write command to output the force and energy values from simulations. I found that all three curves are the same (see two figures below).
Then I compared the values of the thermodynamic output from the simulations with my and hybrid pair styles, and found that they behave as mentioned above (see figures below).
Are you aware how big of an ask it is that you are requesting here??
Basically, this kind of effort is what you could expect from a collaborator that is closely familiar with your research, not from some random folks that volunteer their time to help people get over simple syntax or conceptual errors to get going with their simulations.
It is not like somebody can just look at your materials for a couple of minutes and then immediately point out the mistake and give you instructions for how to correct it. At best we can give some general suggestions.
If you want to compare results, you should always compare them with fix nve and without a thermostat. In fact, you should always use this first to confirm that your simulation settings correctly conserve energy.
You can confirm inconsistencies between the potential energy and the derived forces by using fix numdiff. One common misconception is that the “fpair” variable in LAMMPS pair styles holds the force. It is force/r, so that it only requires multiplication by dx, dy, dz to project the forces on the x-, y-, and z-direction. This step is a big win for the 12-6 Lennard-Jones potential, since it avoids costly divisions.
But as a whole, you have to accept that this is research and thus you are on your own or at best can get help from (local?) collaborators that have a vested interest in your work and your results.
Only if you can provide convincing proof that there is an issue with the distributed version of LAMMPS itself, we would be eager to learn about it and correct it.
Are you aware how big of an ask it is that you are requesting here??
Yes, I understand that but I have no idea what could I do with my problem therefore I decide to write here.
If you want to compare results, you should always compare them with fix nve and without a thermostat. In fact, you should always use this first to confirm that your simulation settings correctly conserve energy.
You can confirm inconsistencies between the potential energy and the derived forces by using fix numdiff. One common misconception is that the “fpair” variable in LAMMPS pair styles holds the force. It is force/r, so that it only requires multiplication by dx, dy, dz to project the forces on the x-, y-, and z-direction. This step is a big win for the 12-6 Lennard-Jones potential, since it avoids costly divisions.
Thank you for your suggestions.
I certainly perform the steps you suggested.
Bottom line is that you (usually with the help of your adviser) will likely need to find a proper collaborator that is sufficiently experienced and interested in what you do. For somebody from the outside, it just takes too much effort to understand the science before you can actually think about tracking down the issues. After all, your issue is first and foremost one of representing the correct science. If you had memory corruption issues or segmentation faults, then some generic debugging techniques can be applied.
While LAMMPS is open source and available at no cost to users, there are some costs involved in solving issues or developing new features. The most sought after commodity is time. And open source software developers are usually volunteering their time or are spending far more time on the project than they are officially entitled to do. At that point, you have to decide what is the most effective way to spend your time. In general, we all want to help everybody and make their experience with LAMMPS a success (otherwise we would not bother in the first place), but we have to make decisions on how to spend the time. This is particularly true, since there is no commonly accepted viable funding model. Neither a “freemium model” or continuously “begging” for donations seems to be working well enough. There are not enough “customers”, and most of them are cash starved in the first place, or else they would probably have looked into using a commercial package.
You did not mention in your initial post whether you had verified your code using nve and nvt integrators first.
I would suggest that you start there. If both of those give identical results then a possible error in your code is incorrect specification of the virial terms. This would mess up the barostat, which would in turn affect ecouple directly, and then the density and potential energy (since the “potential energy” of a condensed system depends on the average interparticle separation and thus the density – scare quotes because the potential energy usually isn’t that meaningful in an MD simulation, but it should be equal between two different implementations of an identical potential function.)
Note that this comment in no way guarantees that there isn’t a problem anywhere else in your code. I’m just thinking through where I might start if I were you.
