forked from lijiext/lammps
Merge pull request #1883 from evoyiatzis/master
Coulomb pair style with smeared out charges (coul/slater)
This commit is contained in:
Axel Kohlmeyer
GitHub
commit
164bf1b60e
|
|
@ -77,6 +77,8 @@ OPT.
|
|||
* :doc:`coul/long/cs (g) <pair_cs>`
|
||||
* :doc:`coul/long/soft (o) <pair_fep_soft>`
|
||||
* :doc:`coul/msm (o) <pair_coul>`
|
||||
* :doc:`coul/slater/cut <pair_coul_slater>`
|
||||
* :doc:`coul/slater/long <pair_coul_slater>`
|
||||
* :doc:`coul/shield <pair_coul_shield>`
|
||||
* :doc:`coul/streitz <pair_coul>`
|
||||
* :doc:`coul/wolf (ko) <pair_coul>`
|
||||
|
|
|
|||
|
|
@ -0,0 +1,121 @@
|
|||
.. index:: pair_style coul/slater
|
||||
|
||||
pair_style coul/slater/cut command
|
||||
==================================
|
||||
|
||||
pair_style coul/slater/long command
|
||||
===================================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style coul/slater/cut lamda cutoff
|
||||
pair_style coul/slater/long lamda cutoff
|
||||
|
||||
lamda = decay length of the charge (distance units)
|
||||
cutoff = cutoff (distance units)
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style coul/slater/cut 1.0 3.5
|
||||
pair_coeff * *
|
||||
pair_coeff 2 2 2.5
|
||||
|
||||
pair_style coul/slater/long 1.0 12.0
|
||||
pair_coeff * *
|
||||
pair_coeff 1 1 5.0
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Styles *coul/slater* compute electrostatic interactions in mesoscopic models
|
||||
which employ potentials without explicit excluded-volume interactions.
|
||||
The goal is to prevent artificial ionic pair formation by including a charge
|
||||
distribution in the Coulomb potential, following the formulation of
|
||||
:ref:`(Melchor) <Melchor>`:
|
||||
|
||||
.. math::
|
||||
|
||||
E = \frac{Cq_iq_j}{\epsilon r} \left( 1- \left( 1 + \frac{r_{ij}}{\lambda} exp\left( -2r_{ij}/\lambda \right) \right) \right) \qquad r < r_c
|
||||
|
||||
|
||||
where :math:`r_c` is the cutoff distance and :math:`\lambda` is the decay length of the charge.
|
||||
C is the same Coulomb conversion factor as in the pair\_styles coul/cut and coul/long. In this way the Coulomb
|
||||
interaction between ions is corrected at small distances r.
|
||||
For the *coul/slater/cut* style, the potential energy for distances larger than the cutoff is zero,
|
||||
while for the *coul/slater/long*, the long-range interactions are computed either by the Ewald or the PPPM technique.
|
||||
|
||||
Phenomena that can be captured at a mesoscopic level using this type of electrostatic
|
||||
interactions include the formation of polyelectrolyte-surfactant aggregates,
|
||||
charge stabilization of colloidal suspensions, and the formation of
|
||||
complexes driven by charged species in biological systems. :ref:`(Vaiwala) <Vaiwala>`.
|
||||
|
||||
The cutoff distance is optional. If it is not used,
|
||||
the default global value specified in the pair_style command is used.
|
||||
For each pair of atom types, a specific cutoff distance can be defined via the :doc:`pair_coeff <pair_coeff>` command as in the example
|
||||
above, or in the data file or restart files read by the
|
||||
:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
|
||||
commands:
|
||||
|
||||
* :math:`r_c` (distance units)
|
||||
|
||||
The global decay length of the charge (:math:`\lambda`) specified in the pair\_style command is used for all pairs.
|
||||
|
||||
|
||||
----------
|
||||
|
||||
|
||||
**Mixing, shift, table, tail correction, restart, rRESPA info**\ :
|
||||
|
||||
For atom type pairs I,J and I != J, the cutoff distance for the
|
||||
*coul/slater* styles can be mixed. The default mix value is *geometric*\ .
|
||||
See the "pair\_modify" command for details.
|
||||
|
||||
The :doc:`pair_modify <pair_modify>` shift and table options are not relevant
|
||||
for these pair styles.
|
||||
|
||||
These pair styles do not support the :doc:`pair_modify <pair_modify>`
|
||||
tail option for adding long-range tail corrections to energy and
|
||||
pressure.
|
||||
|
||||
These pair styles write their information to :doc:`binary restart files <restart>`, so pair\_style and pair\_coeff commands do not need
|
||||
to be specified in an input script that reads a restart file.
|
||||
|
||||
This pair style can only be used via the *pair* keyword of the
|
||||
:doc:`run_style respa <run_style>` command. It does not support the
|
||||
*inner*\ , *middle*\ , *outer* keywords.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
The *coul/slater/long* style requires the long-range solvers included in the KSPACE package.
|
||||
|
||||
These styles are part of the "USER-MISC" package. They are only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`pair_coeff <pair_coeff>`, :doc:`pair_style, hybrid/overlay <pair_hybrid>`, :doc:`kspace_style <kspace_style>`
|
||||
|
||||
**Default:** none
|
||||
|
||||
----------
|
||||
|
||||
|
||||
.. _Melchor:
|
||||
|
||||
**(Melchor)** Gonzalez-Melchor, Mayoral, Velázquez, and Alejandre, J Chem Phys, 125, 224107 (2006).
|
||||
|
||||
.. _Vaiwala:
|
||||
|
||||
**(Vaiwala)** Vaiwala, Jadhav, and Thaokar, J Chem Phys, 146, 124904 (2017).
|
||||
|
||||
|
||||
|
|
@ -1027,6 +1027,7 @@ Gmask
|
|||
gneb
|
||||
GNEB
|
||||
Goldfarb
|
||||
Gonzalez-Melchor
|
||||
googlemail
|
||||
Gordan
|
||||
GPa
|
||||
|
|
@ -1290,6 +1291,7 @@ Izvekov
|
|||
izz
|
||||
Izz
|
||||
Jacobsen
|
||||
Jadhav
|
||||
jagreat
|
||||
Jalalvand
|
||||
james
|
||||
|
|
@ -1663,6 +1665,7 @@ maxwell
|
|||
Maxwellian
|
||||
maxX
|
||||
Mayergoyz
|
||||
Mayoral
|
||||
mbt
|
||||
Mbytes
|
||||
MBytes
|
||||
|
|
@ -1691,6 +1694,7 @@ mediumvioletred
|
|||
Mees
|
||||
Mehl
|
||||
Mei
|
||||
Melchor
|
||||
Meloni
|
||||
Melrose
|
||||
Mem
|
||||
|
|
@ -2260,6 +2264,7 @@ polyA
|
|||
polybond
|
||||
polydisperse
|
||||
polydispersity
|
||||
polyelectrolyte
|
||||
polyhedra
|
||||
popen
|
||||
Popov
|
||||
|
|
@ -2769,6 +2774,7 @@ superset
|
|||
supersphere
|
||||
Supinski
|
||||
Surblys
|
||||
surfactant
|
||||
surfactants
|
||||
Suter
|
||||
Sutmann
|
||||
|
|
@ -2831,6 +2837,7 @@ tfmc
|
|||
tfMC
|
||||
th
|
||||
Thakkar
|
||||
Thaokar
|
||||
thb
|
||||
thei
|
||||
Theodorou
|
||||
|
|
@ -3025,6 +3032,7 @@ uwo
|
|||
Uzdin
|
||||
vacf
|
||||
Vaid
|
||||
Vaiwala
|
||||
valent
|
||||
Valeriu
|
||||
valgrind
|
||||
|
|
@ -3053,6 +3061,7 @@ Vectorization
|
|||
vectorized
|
||||
Vegt
|
||||
vel
|
||||
Velázquez
|
||||
Verlag
|
||||
verlet
|
||||
Verlet
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,60 @@
|
|||
# Bulk polyelectrolyte as described in section IV of J. Chem. Phys. 125, 224107 (2006)
|
||||
|
||||
boundary p p p
|
||||
|
||||
units lj
|
||||
atom_style charge
|
||||
|
||||
region my_sim_box block 0.0 10.0 0.0 10.0 0.0 10.0
|
||||
create_box 3 my_sim_box
|
||||
|
||||
create_atoms 1 random 2804 100 my_sim_box
|
||||
create_atoms 2 random 98 200 my_sim_box
|
||||
create_atoms 3 random 98 300 my_sim_box
|
||||
|
||||
set type 2 charge -1.0
|
||||
set type 3 charge 1.0
|
||||
|
||||
comm_modify mode single vel yes
|
||||
|
||||
mass 1 1.0
|
||||
mass 2 1.0
|
||||
mass 3 1.0
|
||||
|
||||
pair_style hybrid/overlay dpd 1.0 1.0 245455 coul/slater/long 0.929 3.0
|
||||
pair_coeff * * dpd 25.0 4.5
|
||||
pair_coeff * * coul/slater/long
|
||||
|
||||
kspace_style ewald 0.00001
|
||||
dielectric 1.0
|
||||
|
||||
neighbor 2.0 bin
|
||||
neigh_modify every 1 delay 0 check no once no
|
||||
|
||||
timestep 0.02
|
||||
|
||||
fix 2 all nve
|
||||
|
||||
thermo 10
|
||||
thermo_style custom step spcpu temp press etotal pe ke ecoul elong evdwl
|
||||
thermo_modify line one
|
||||
|
||||
run 100000
|
||||
write_data data.after_equilibration
|
||||
|
||||
compute RDF_1_1 all rdf 50 1 1 cutoff 3.0
|
||||
compute RDF_1_2 all rdf 50 1 2 cutoff 3.0
|
||||
compute RDF_1_3 all rdf 50 1 3 cutoff 3.0
|
||||
compute RDF_2_2 all rdf 50 2 2 cutoff 3.0
|
||||
compute RDF_2_3 all rdf 50 2 3 cutoff 3.0
|
||||
compute RDF_3_3 all rdf 50 3 3 cutoff 3.0
|
||||
|
||||
fix 11 all ave/time 50 1 50 c_RDF_1_1[*] file tmp_1_1.rdf mode vector
|
||||
fix 12 all ave/time 50 1 50 c_RDF_1_2[*] file tmp_1_2.rdf mode vector
|
||||
fix 13 all ave/time 50 1 50 c_RDF_1_3[*] file tmp_1_3.rdf mode vector
|
||||
fix 14 all ave/time 50 1 50 c_RDF_2_2[*] file tmp_2_2.rdf mode vector
|
||||
fix 15 all ave/time 50 1 50 c_RDF_2_3[*] file tmp_2_3.rdf mode vector
|
||||
fix 16 all ave/time 50 1 50 c_RDF_3_3[*] file tmp_3_3.rdf mode vector
|
||||
|
||||
run 10000
|
||||
write_data data.after_production_run
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
|
@ -76,6 +76,8 @@ pair_style buck/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15
|
|||
pair_style cosine/squared, Eugen Rozic, eugen.rozic.17 at ucl.ac.uk, 9 Aug 19
|
||||
pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11
|
||||
pair_style coul/shield, Wengen Ouyang (Tel Aviv University), w.g.ouyang at gmail dot com, 30 Mar 18
|
||||
pair_style coul/slater/cut, Evangelos Voyiatzis, evoyiatzis at gmail.com, 26 February 2020
|
||||
pair_style coul/slater/long, Evangelos Voyiatzis, evoyiatzis at gmail.com, 26 February 2020
|
||||
pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11
|
||||
pair_style e3b, Steven Strong (U Chicago), stevene.strong at gmail dot com, 16 Apr 19
|
||||
pair_style drip, Mingjian Wen, University of Minnesota, wenxx151 at umn.edu, 17 Apr 19
|
||||
|
|
|
|||
|
|
@ -0,0 +1,187 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
* Contributing author: Evangelos Voyiatzis (Royal DSM)
|
||||
* ------------------------------------------------------------------------- */
|
||||
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <cstring>
|
||||
#include "pair_coul_slater_cut.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulSlaterCut::PairCoulSlaterCut(LAMMPS *lmp) : PairCoulCut(lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterCut::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
double rsq,r2inv,r,rinv,forcecoul,factor_coul,bracket_term;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
rinv = 1.0/r;
|
||||
bracket_term = 1 - exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
|
||||
forcecoul = qqrd2e * scale[itype][jtype] *
|
||||
qtmp*q[j] * bracket_term * rinv;
|
||||
fpair = factor_coul*forcecoul * r2inv;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) ecoul = factor_coul * qqrd2e *
|
||||
scale[itype][jtype] * qtmp*q[j] * rinv *
|
||||
(1 - (1 + r/lamda)*exp(-2*r/lamda));
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
0.0,ecoul,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterCut::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 2) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
lamda = force->numeric(FLERR,arg[0]);
|
||||
cut_global = force->numeric(FLERR,arg[1]);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i,j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j]) cut[i][j] = cut_global;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterCut::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_global,sizeof(double),1,fp);
|
||||
fwrite(&lamda,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterCut::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_global,sizeof(double),1,fp);
|
||||
fread(&lamda,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulSlaterCut::single(int i, int j, int /*itype*/, int /*jtype*/,
|
||||
double rsq, double factor_coul, double /*factor_lj*/,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r,rinv,forcecoul,phicoul,bracket_term;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
r = sqrt(rsq);
|
||||
rinv = 1.0/r;
|
||||
bracket_term = 1 - exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
|
||||
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] *
|
||||
bracket_term * rinv;
|
||||
fforce = factor_coul*forcecoul * r2inv;
|
||||
|
||||
phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * rinv * (1 - (1 + r/lamda)*exp(-2*r/lamda));
|
||||
return factor_coul*phicoul;
|
||||
}
|
||||
|
|
@ -0,0 +1,53 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(coul/slater/cut,PairCoulSlaterCut)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_COUL_SLATER_CUT_H
|
||||
#define LMP_PAIR_COUL_SLATER_CUT_H
|
||||
|
||||
#include "pair_coul_cut.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairCoulSlaterCut : public PairCoulCut {
|
||||
public:
|
||||
PairCoulSlaterCut(class LAMMPS *);
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
double single(int, int, int, int, double, double, double, double &);
|
||||
|
||||
protected:
|
||||
double lamda;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,422 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
* Contributing author: Evangelos Voyiatzis (Royal DSM)
|
||||
* ------------------------------------------------------------------------- */
|
||||
|
||||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <cstring>
|
||||
#include "pair_coul_slater_long.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "kspace.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "update.h"
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define EWALD_F 1.12837917
|
||||
#define EWALD_P 0.3275911
|
||||
#define A1 0.254829592
|
||||
#define A2 -0.284496736
|
||||
#define A3 1.421413741
|
||||
#define A4 -1.453152027
|
||||
#define A5 1.061405429
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulSlaterLong::PairCoulSlaterLong(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
ewaldflag = pppmflag = 1;
|
||||
//ftable = NULL;
|
||||
qdist = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairCoulSlaterLong::~PairCoulSlaterLong()
|
||||
{
|
||||
if (!copymode) {
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
|
||||
memory->destroy(scale);
|
||||
}
|
||||
//if (ftable) free_tables();
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
|
||||
// double fraction,table;
|
||||
double r,r2inv,forcecoul,factor_coul;
|
||||
double grij,expm2,prefactor,t,erfc;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
double rsq;
|
||||
double slater_term;
|
||||
// int itable;
|
||||
|
||||
ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double *q = atom->q;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_coul = force->special_coul;
|
||||
int newton_pair = force->newton_pair;
|
||||
double qqrd2e = force->qqrd2e;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
qtmp = q[i];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_coul = special_coul[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cut_coulsq) {
|
||||
r2inv = 1.0/rsq;
|
||||
// if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
|
||||
prefactor = qqrd2e * scale[itype][jtype] * qtmp*q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
/*
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
}
|
||||
*/
|
||||
|
||||
fpair = forcecoul * r2inv;
|
||||
|
||||
f[i][0] += delx*fpair;
|
||||
f[i][1] += dely*fpair;
|
||||
f[i][2] += delz*fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx*fpair;
|
||||
f[j][1] -= dely*fpair;
|
||||
f[j][2] -= delz*fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
// if (!ncoultablebits || rsq <= tabinnersq)
|
||||
ecoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
|
||||
/*
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
ecoul = scale[itype][jtype] * qtmp*q[j] * table;
|
||||
}
|
||||
*/
|
||||
if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i,j,nlocal,newton_pair,
|
||||
0.0,ecoul,fpair,delx,dely,delz);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes;
|
||||
|
||||
memory->create(setflag,n+1,n+1,"pair:setflag");
|
||||
for (int i = 1; i <= n; i++)
|
||||
for (int j = i; j <= n; j++)
|
||||
setflag[i][j] = 0;
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
memory->create(scale,n+1,n+1,"pair:scale");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 2) error->all(FLERR,"Illegal pair_style command");
|
||||
|
||||
lamda = force->numeric(FLERR,arg[0]);
|
||||
cut_coul = force->numeric(FLERR,arg[1]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 2) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
scale[i][j] = 1.0;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::init_style()
|
||||
{
|
||||
if (!atom->q_flag)
|
||||
error->all(FLERR,"Pair style coul/slater/long requires atom attribute q");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
cut_coulsq = cut_coul * cut_coul;
|
||||
|
||||
// insure use of KSpace long-range solver, set g_ewald
|
||||
|
||||
if (force->kspace == NULL)
|
||||
error->all(FLERR,"Pair style requires a KSpace style");
|
||||
g_ewald = force->kspace->g_ewald;
|
||||
|
||||
// setup force tables
|
||||
|
||||
// if (ncoultablebits) init_tables(cut_coul,NULL);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulSlaterLong::init_one(int i, int j)
|
||||
{
|
||||
scale[j][i] = scale[i][j];
|
||||
return cut_coul+2.0*qdist;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
for (int j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j],sizeof(int),1,fp);
|
||||
if (setflag[i][j])
|
||||
fwrite(&scale[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
|
||||
allocate();
|
||||
|
||||
int i,j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
|
||||
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
|
||||
if (setflag[i][j]) {
|
||||
if (me == 0) fread(&scale[i][j],sizeof(double),1,fp);
|
||||
MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_coul,sizeof(double),1,fp);
|
||||
fwrite(&lamda,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
//fwrite(&ncoultablebits,sizeof(int),1,fp);
|
||||
//fwrite(&tabinner,sizeof(double),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairCoulSlaterLong::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_coul,sizeof(double),1,fp);
|
||||
fread(&lamda,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
//fread(&ncoultablebits,sizeof(int),1,fp);
|
||||
//fread(&tabinner,sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&lamda,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
//MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world);
|
||||
//MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairCoulSlaterLong::single(int i, int j, int /*itype*/, int /*jtype*/,
|
||||
double rsq,
|
||||
double factor_coul, double /*factor_lj*/,
|
||||
double &fforce)
|
||||
{
|
||||
double r2inv,r,grij,expm2,t,erfc,prefactor;
|
||||
double slater_term;
|
||||
// double fraction,table;
|
||||
double forcecoul,phicoul;
|
||||
// int itable;
|
||||
|
||||
r2inv = 1.0/rsq;
|
||||
// if (!ncoultablebits || rsq <= tabinnersq) {
|
||||
r = sqrt(rsq);
|
||||
grij = g_ewald * r;
|
||||
expm2 = exp(-grij*grij);
|
||||
t = 1.0 / (1.0 + EWALD_P*grij);
|
||||
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
|
||||
slater_term = exp(-2*r/lamda)*(1 + (2*r/lamda*(1+r/lamda)));
|
||||
prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
|
||||
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2 - slater_term);
|
||||
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
/*
|
||||
} else {
|
||||
union_int_float_t rsq_lookup;
|
||||
rsq_lookup.f = rsq;
|
||||
itable = rsq_lookup.i & ncoulmask;
|
||||
itable >>= ncoulshiftbits;
|
||||
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
|
||||
table = ftable[itable] + fraction*dftable[itable];
|
||||
forcecoul = atom->q[i]*atom->q[j] * table;
|
||||
if (factor_coul < 1.0) {
|
||||
table = ctable[itable] + fraction*dctable[itable];
|
||||
prefactor = atom->q[i]*atom->q[j] * table;
|
||||
forcecoul -= (1.0-factor_coul)*prefactor;
|
||||
}
|
||||
}
|
||||
*/
|
||||
fforce = forcecoul * r2inv;
|
||||
|
||||
// if (!ncoultablebits || rsq <= tabinnersq)
|
||||
phicoul = prefactor*(erfc - (1 + r/lamda)*exp(-2*r/lamda));
|
||||
/*
|
||||
else {
|
||||
table = etable[itable] + fraction*detable[itable];
|
||||
phicoul = atom->q[i]*atom->q[j] * table;
|
||||
}
|
||||
*/
|
||||
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
|
||||
|
||||
return phicoul;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairCoulSlaterLong::extract(const char *str, int &dim)
|
||||
{
|
||||
if (strcmp(str,"cut_coul") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &cut_coul;
|
||||
}
|
||||
if (strcmp(str,"lamda") == 0) {
|
||||
dim = 0;
|
||||
return (void *) &lamda;
|
||||
}
|
||||
if (strcmp(str,"scale") == 0) {
|
||||
dim = 2;
|
||||
return (void *) scale;
|
||||
}
|
||||
return NULL;
|
||||
}
|
||||
|
|
@ -0,0 +1,78 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
|
||||
PairStyle(coul/slater/long,PairCoulSlaterLong)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_COUL_SLATER_LONG_H
|
||||
#define LMP_PAIR_COUL_SLATER_LONG_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairCoulSlaterLong : public Pair {
|
||||
public:
|
||||
PairCoulSlaterLong(class LAMMPS *);
|
||||
~PairCoulSlaterLong();
|
||||
virtual void compute(int, int);
|
||||
virtual void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
virtual void init_style();
|
||||
virtual double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
virtual void write_restart_settings(FILE *);
|
||||
virtual void read_restart_settings(FILE *);
|
||||
virtual double single(int, int, int, int, double, double, double, double &);
|
||||
virtual void *extract(const char *, int &);
|
||||
|
||||
protected:
|
||||
double cut_coul,cut_coulsq,qdist;
|
||||
double lamda;
|
||||
//double *cut_respa;
|
||||
double g_ewald;
|
||||
double **scale;
|
||||
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal ... command
|
||||
|
||||
Self-explanatory. Check the input script syntax and compare to the
|
||||
documentation for the command. You can use -echo screen as a
|
||||
command-line option when running LAMMPS to see the offending line.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair style coul/slater/long requires atom attribute q
|
||||
|
||||
The atom style defined does not have this attribute.
|
||||
|
||||
E: Pair style requires a KSpace style
|
||||
|
||||
No kspace style is defined.
|
||||
|
||||
*/
|
||||
Loading…
Reference in New Issue