Thermostat and barostat with fix npt

Dear all,

I met a problem with thermostat. I developed a fix style wich adds inter-atom force accounting for polarization effect by post-force method. Induced energy is added to the system by “fix_modify energy yes” command and compute_scalar method. Additional virial is also added to the system. I’m pretty sure there is no wrong formula or bug in my code. (I have verified the code by running in nve ensemble, It ran well except the energy is not conserved. I think it is because of the cut-off method to calculate the interaction added by my fix style which doesn’t account for the long range force about 1/r^4 ). When I ran it in NPT ensemble with fix npt and set T=298K p=1atm rho=1g/cm^3, the temperature was always around 315 K which is much higher than 298K and the density is about 1.25 g/cm^3. I ran it for 1ns with 1fs time step. I think it should be enough to get equilibrium, or at least the temperature shouldn’t be like this. I know my problem is not completely related to LAMMPS but I’m curious if I forget something in order to add such inter-atom interaction. I check the code of pair style and think force\energy\virial are enough to achieve such a cut-off inter-atom interaction force field.

Thanks in advance

​if you cannot conserve energy in NVE, how can you expect to get meaningful
results with NPT?

in general, you have to​:

Thanks, Axel. I can conserve the energy in NVE ensemble without my fix. But as I said, since my fix introduces cut-off model to calculate the added force, I don’t think I can conserve the energy with my fix. I have tried to reduce the time step from1 to 0.5 and 0.05 for nve but the energy keeps growing. I think it’s just like the energy is not conserved with cut-off model for coulomb energy without considering long range fore by kspace methods. If my understanding is incorrect, please point out. It’s currently hard for me to implement long range calculation for my fix. But any advice on such problem is appreciated.

The energy is almost identical except small difference in decimal digits when running on different # of procs. I have printed the communicated information to files and they are correct.

Thanks, Axel. I can conserve the energy in NVE ensemble without my fix.
But as I said, since my fix introduces cut-off model to calculate the added
force, I don't think I can conserve the energy with my fix. I have tried to
reduce the time step from1 to 0.5 and 0.05 for nve but the energy keeps
growing. I think it's just like the energy is not conserved with cut-off
model for coulomb energy without considering long range fore by kspace
methods.

​what is the justification for this statement? ​

If my understanding is incorrect, please point out. It's currently hard for

me to implement long range calculation for my fix. But any advice on such
problem is appreciated.

​i cannot comment on your code.​ i know from personal experience, that
implementing polarization is a tricky thing (i did that with my adviser's
code 20 years ago) and there are many subtle pitfalls. however, 20 years is
a long time and a lot has happened since then, so my recollections are
hazy. but it is easy to have a sign error or to do the manipulations in
your fix in the wrong fix method, or it may simply be required to define
fixes in a particular order.

​if you cannot handle kspace, the you should test your system completely
without.​

​in general, it is almost impossible to give specific advice when knowing
next to no details. my only advice that is independent of this is:
testing and debugging is 90% of the time of developing code​. don't try to
skip out of it and remain humble and paranoid. bugs are lurking everywhere,
even if they are only a simple flip of the sign.

​axel.​