Recently I have been aiming to implement a LJ wall in a python package unrelated to LAMMPS, the wall seemed to work fine but to ensure the validity of my implementation I started to run NVE simulation using a reactive forcefield (MACE), that’s when problems started to kick in.
I realised that when interacting with the wall the total energy is not conserved, even when adding the LJ energy term the energy would fluctuate in the order of ~ 10^-5 eV. Worse, this would still continue to happen after the interaction with the wall, when molecules are far away. To be fair, energy conservation with MACE seems to be particularly poor.
I ran similar simulations (3 water with initial velocities to hit the wall) with LAMMPS (ReaxFF) and it seems that the problem is the same: energy is poorly conserved even after the interactions.
I have a few questions:
Is the energy interaction between wall/particles included in the “TotEng” printed in log.lammps? I am 99% sure the answer is no, just making sure.
Why is the energy conservation perturbed even after the interactions?
My python implementation can be summarised as: adding the additional LJ force terms in verlet (F(t) and F(t + dt)), any comment would be welcomed.
See fix wall documentation:
The fix_modify energy option is supported by this fix to add the energy of interaction between atoms and all the specified walls to the global potential energy of the system as part of thermodynamic output. The default setting for this fix is fix_modify energy no.
Impossible to tell without more insights on your inputs. Especially with ReaxFF which is a very complex forcefield and very sensitive to parameters change.
I don’t understand how your Verlet implementation works. F(t) is used to compute r(t+dt) from which F(t+dt) is computed. I don’t understand why you mention adding the LJ force to both F(t) and F(t+dt)
Additionally reactive forcefields are already kind of a complex stuff to “ensure the validity of your implementation”. The energy conservation is a complicated issue in MD in general because of energy drifts due to numerical integration. Reactive forcefields require notoriously low timestep because of all the complex instantaneous local evaluations of forces that have to be done during reactions to take into account the very fast pace of bond changes (normally at the electron scales). You have to be sure that your implementation and your parameters (like timestep amplitude) already conserve energy with reasonable drift before going into complex comparisons.
Very sorry, I usually make sure to read the documentation thoroughly before posting, this one escaped my attention.
The input for testing lookw like this:
boundary p p m
units real
atom_style charge
read_data mc.lmp_data
pair_style reaxff NULL
pair_coeff * * ./ffield.reax.pt_o_h H O
fix q_equi all qeq/reax 1 0.0 10.0 1e-6 reax/c
fix wallhi all wall/lj126 zhi 24 0.01 1.0 3.0
thermo_style custom step temp density etotal
thermo 1
timestep 0.25
fix 1 all nve
velocity all set 0.0 0.0 0.005
dump 50 all custom 100 dump.lammpstrj id type x y z
run 100000
The system is a box of 3 water molecules, each will hit the wall, reflect and go away.
Yes sorry, due to the specifics of this implementation I worded my sentence wrong, what I meant is that LJ forces are added when computing the forces and nothing else. (I did not modify Verlet in any way)
I’m not sure energy is conserved to high precision (or at all) in charge equilibration MD.
After all, if the charges don’t change much, then charge equilibration is unnecessary. If the charges do change a lot, then the changed charge must appreciably change the energy of the system (since the new charges are calculated to minimise electrostatic energy, albeit “dressed” by electronegativities).