Hello all,
I have a 35 Angstrom width (z direction) semi-periodic channel (x and y) containing SPC/E water and ions. Channel walls are modeled as discrete atoms on a rigid lattice with partial charges to balance the excess ion charge of the system, and so
Nwall * qwall + Nion * qion = 0,
is satisfied. For simplicity, consider only one ionic species. I wish to compute the partial force on a water molecule due to its interaction with the wall. The force would contain LJ interaction + electrostatics due to the discrete nature of the SPC/E water, and its preferred orientation near a charged wall.
A full trajectory is generated where all interactions are present, which contains the following forces on the water molecule
Wall fluid (WF) + Fluid Fluid (FF) + Ion Fluid (IF)
The form of the interaction is lj/cut/coul/long with pppm style for long-range electrostatics.
A rerun is performed (using lj/cut/coul/long) where WF and IF Lennard-Jones parameters are set 0, along with setting the wall and ion charges 0 simultaneously. This results in only FF component that can be removed to obtain WF + IF forces.
If the ion charges are set to 0 to compute only the WF component, the global charge neutrality breaks, that may result in serious artifacts in the computed forces. Can anyone suggest some appropriate long range correction that can mitigate this issue? The imbalance is +/- 32e at maximum.
Thanks,
Ravi Bhadauria
LAMMPS applies a correction automatically in PPPM for non-neutral systems to the energy, but not to the forces. Either the correction to the forces is not necessary or not easily derived, I’m not sure which. You could also take a look at compute group/group here, but it probably won’t work with SHAKE http://lammps.sandia.gov/doc/compute_group_group.html:
Normally the long-range Coulombic energy converges only when the net charge of the unit cell is zero. However, one can assume the net charge of the system is neutralized by a uniform background plasma, and a correction to the system energy can be applied to reduce artifacts. For more information see (Bogusz). If the boundary keyword is set to yes, which is the default, and kspace contributions are included, then this energy correction term will be added to the total group-group energy. This correction term does not affect the force calculation and will be zero if one or both of the groups are charge neutral. This energy correction term is the same as that included in the regular Ewald and PPPM routines.
Hello all,
I have a 35 Angstrom width (z direction) semi-periodic channel (x and y)
containing SPC/E water and ions. Channel walls are modeled as discrete atoms
on a rigid lattice with partial charges to balance the excess ion charge of
the system, and so
Nwall * qwall + Nion * qion = 0,
is satisfied. For simplicity, consider only one ionic species. I wish to
compute the partial force on a water molecule due to its interaction with
the wall. The force would contain LJ interaction + electrostatics due to the
discrete nature of the SPC/E water, and its preferred orientation near a
charged wall.
A full trajectory is generated where all interactions are present, which
contains the following forces on the water molecule
Wall fluid (WF) + Fluid Fluid (FF) + Ion Fluid (IF)
The form of the interaction is lj/cut/coul/long with pppm style for
long-range electrostatics.
A rerun is performed (using lj/cut/coul/long) where WF and IF Lennard-Jones
parameters are set 0, along with setting the wall and ion charges 0
simultaneously. This results in only FF component that can be removed to
obtain WF + IF forces.
If the ion charges are set to 0 to compute only the WF component, the global
charge neutrality breaks, that may result in serious artifacts in the
computed forces. Can anyone suggest some appropriate long range correction
that can mitigate this issue? The imbalance is +/- 32e at maximum.
since you are doing a re-run to compute this, i would recommend to
rather drop PPPM and use lj/cut/coul/cut with a sufficiently large
cutoff. in that case your "f" boundary is going to be a real "f"
boundary and no inter-slab dipole correction is required.
the effects on forces, that usually keep people from using this
setting for a good reason during regular runs, don't apply when doing
a rerun, and you avoid a whole lot of complications with sources of
errors, that are much easier to control and determine. doing this kind
of partitioning will incur errors, anyway. and you can empirically
check on the convergence with respect to the cutoff easily, too.
axel.
This would work well for z forces, where in the original run I have kept p p p boundaries with large inter-slab separation (250 Angstroms). However, the physical periodicity in x and y (38 Ang X 36 Ang) dimension would restrict cutoff to 18 Ang which may not be sufficient.
One way to go about this would be increasing the periodic box dimensions in x and y directions while doing the full run with PPPM, so that a large cutoff can be utilized and empirical convergence can be checked during reruns. Can you suggest a ballpark on what should be the minimum value of cut-off to begin with?
Thanks,
This would work well for z forces, where in the original run I have kept p p
p boundaries with large inter-slab separation (250 Angstroms). However, the
that means, you wasted a lot of CPU time.
physical periodicity in x and y (38 Ang X 36 Ang) dimension would restrict
cutoff to 18 Ang which may not be sufficient.
no, this restriction doesn't apply with LAMMPS. you can have cutoffs
larger than half the box. but as i already noted, since you are not
using the forces for MD, you may not even need to converge them as
tightly.
and even 20 \AA or more may be overkill. have you made a test for how
well you need to converge your force decomposition?
One way to go about this would be increasing the periodic box dimensions in
x and y directions while doing the full run with PPPM, so that a large
cutoff can be utilized and empirical convergence can be checked during
reruns. Can you suggest a ballpark on what should be the minimum value of
cut-off to begin with?
it looks like you haven't really understood what i was suggesting and
why during rerun a lot of things are acceptable, that won't do during
production run.
axel.