Dear LAMMPS Users,
I am currently trying to make a data file which contains proteins and TIP4P/2005 waters. I know there are TIP4P-specific options such as lj/cut/tip4p/long or pppm/tip4p without explicitly specifying the coordinates of dummy site in the date file. However, for a specific reason, I would like to include the coordinates of dummy sites in the data file and just use regular options of lj/cut/coul/long and pppm. And for my purpose, I just need to generate a configuration for an analysis purpose (not planning on running simulation with this configuration).
In order to check if I have generated a correct file, I check energies (short-range LJ, short-range coulomb, long-range coulomb, etc.). I check by comparing with the energies calculated by using GROMACS and in-house code. The only discrepancy I am seeing is with the electrostatic energies when I use the Ewald or PPPM. If I just use cut-off at a certain distance (for electrostatics), the energies match very well, but only when I use Ewald or PPPM there are discrepancies in electrostatic energies (both short and long ranges).
I specifically need to use lj/cut/coul/long option with either ewald or pppm, and I was wondering what I am setting wrong. The proteins are no problem. Any suggestions on this issue would be greatly appreciated. The following is the input I am using:
Dear LAMMPS Users,
[...]
I specifically need to use lj/cut/coul/long option with either ewald or
pppm, and I was wondering what I am setting wrong. The proteins are no
problem. Any suggestions on this issue would be greatly appreciated. The
following is the input I am using:
make a less complex input, e.g. one that you can compute by hand and
then you should *much* more carefully investigate how to synchronize
the settings for long-range electrostatics to make those codes match
(close enough). there may also be implementation specific differences
(gromacs applies its cutoff to charge groups, not individual atoms)
also PPPM is different from (S)PME and different codes use different
estimators for kspace parameters. 1.-e-4 is very sloppy. you are
likely to only get agreement, if you choose a very narrow kspace
convergence in all codes.
-------------------------------------------------------------------------------------------------------------
special_bonds amber
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style none
pair_style lj/cut/coul/long 10.00000 10.000000
pair_modify mix arithmetic
kspace_style ewald 1.0e-4
dielectric 1.000000
# neighbor list update
neighbor 0.000000 bin
neigh_modify every 1 delay 0 check no once no
# water (20 – oxygen, 21 – hydrogen, 22 – dummy site)
group water type 20 21 22
neigh_modify exclude molecule water
fix 1 water shake 0.0001 20 0 b 26 a 38
shake doesn't work for tip4p with explicit M. and for a rerun it is
pointless altogether.
axel