Hi
While this my seem like a strange question it is the final problem we are having before we can submit our new reactive model to lammps so I would be very grateful for some help here as I am completely stuck.
I need to create some sort of bond list within my code to screen out charges. Specifically, I need to remove 1-2 and 1-3 bonded neighbors from the computation of charge. At first I thought this would be a straightforward task make the list, create a correction, and match the forces (and energy). However, in each of the ways I tried there seems to be a problem with the way my code interacts with the larger lammps code.
First I tried a direct subtraction removing the charge and force components from the bond energy. This produced a very strange and uncorrectable behavior in the computation of pressure.
Second, I tried implementing coul/cut directly in my pair style. This worked best of everything but was far too slow requiring cutoffs greater than 15Å for reasonable results.
Then I tried a direct computation using wolf charges first as a contribution to my forces and energy, which produced the same pressure problem as in step 1.
Fourth I tried using it is a separate unit operating independently (its own evtally) but this altered the energy and the forces calculated. The pressures were actually calculating reasonably.
In all these cases the dynamics is problematic, the error in the pressure does seem to reflect a real problem with the computation of forces, even though when I output them in a dump command they appear correct.
So given these are all in retrospect bad ideas, what is a good way to go about this?
In particular does it make sense to access sbmask? If so where in the code that I could find its particulars?
Or does it make sense to go through the bond creation and destruction routines that can be referenced from the input? Is there a way to access them internally (i.e. from within my pair style)?
Clearly I am missing some basic understanding of the operation of the code.
thanks
matthew