- cutoff = global cutoff for Coulombic interactions +
- kappa = Debye length (inverse distance units)
Examples:
@@ -26,6 +30,10 @@ pair_coeff * * pair_coeff 2 2 3.5 +pair_style coul/debye 1.4 3.0 +pair_coeff * * +pair_coeff 2 2 3.5 +
pair_style coul/long 10.0 pair_coeff * *@@ -41,6 +49,14 @@ by the dielectric command. The cutoff Rc truncates the interaction distance. +
Style coul/debye adds an additional exp() damping factor to the +Coulombic term, given by +
+
+where kappa is the Debye length. This potential is another way to +mimic the screening effect of a polar solvent. +
Style coul/long computes the same Coulombic interactions as style coul/cut except that an additional damping factor is applied so it can be used in conjunction with the kspace_style @@ -63,9 +79,9 @@
- cutoff (distance units)
For coul/cut, the cutoff coefficient is optional. If it is not used -(as in som of the examples above), the default global value specified -in the pair_style command is used. +
For coul/cut and coul/debye, the cutoff coefficient is optional. +If it is not used (as in some of the examples above), the default +global value specified in the pair_style command is used.
For coul/long no cutoff can be specified for an individual I,J type pair via the pair_coeff command. All type pairs use the same global diff -Naur lammps-5Jan08/doc/pair_coul.txt lammps-6Jan08/doc/pair_coul.txt --- lammps-5Jan08/doc/pair_coul.txt 2007-10-11 17:09:49.000000000 -0600 +++ lammps-6Jan08/doc/pair_coul.txt 2008-01-03 18:38:02.000000000 -0700 @@ -7,14 +7,17 @@ :line pair_style coul/cut command :h3 +pair_style coul/debye command :h3 pair_style coul/long command :h3 [Syntax:] pair_style coul/cut cutoff +pair_style coul/debye kappa cutoff pair_style coul/long cutoff :pre -cutoff = global cutoff for Coulombic interactions :ul +cutoff = global cutoff for Coulombic interactions +kappa = Debye length (inverse distance units) :ul [Examples:] @@ -22,6 +25,10 @@ pair_coeff * * pair_coeff 2 2 3.5 :pre +pair_style coul/debye 1.4 3.0 +pair_coeff * * +pair_coeff 2 2 3.5 :pre + pair_style coul/long 10.0 pair_coeff * * :pre @@ -37,6 +44,14 @@ by the "dielectric"_dielectric.html command. The cutoff Rc truncates the interaction distance. +Style {coul/debye} adds an additional exp() damping factor to the +Coulombic term, given by + +:c,image(Eqs/pair_debye.jpg) + +where kappa is the Debye length. This potential is another way to +mimic the screening effect of a polar solvent. + Style {coul/long} computes the same Coulombic interactions as style {coul/cut} except that an additional damping factor is applied so it can be used in conjunction with the "kspace_style"_kspace_style.html @@ -59,9 +74,9 @@ cutoff (distance units) :ul -For {coul/cut}, the cutoff coefficient is optional. If it is not used -(as in som of the examples above), the default global value specified -in the pair_style command is used. +For {coul/cut} and {coul/debye}, the cutoff coefficient is optional. +If it is not used (as in some of the examples above), the default +global value specified in the pair_style command is used. For {coul/long} no cutoff can be specified for an individual I,J type pair via the pair_coeff command. All type pairs use the same global diff -Naur lammps-5Jan08/src/pair_coul_cut.h lammps-6Jan08/src/pair_coul_cut.h --- lammps-5Jan08/src/pair_coul_cut.h 2008-01-02 12:24:46.000000000 -0700 +++ lammps-6Jan08/src/pair_coul_cut.h 2008-01-03 18:20:28.000000000 -0700 @@ -21,19 +21,19 @@ class PairCoulCut : public Pair { public: PairCoulCut(class LAMMPS *); - ~PairCoulCut(); - void compute(int, int); - void settings(int, char **); + virtual ~PairCoulCut(); + virtual void compute(int, int); + virtual void settings(int, char **); void coeff(int, char **); void init_style(); double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); - void write_restart_settings(FILE *); - void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); + virtual void write_restart_settings(FILE *); + virtual void read_restart_settings(FILE *); + virtual double single(int, int, int, int, double, double, double, double &); - private: + protected: double cut_global; double **cut; diff -Naur lammps-5Jan08/src/pair_coul_debye.cpp lammps-6Jan08/src/pair_coul_debye.cpp --- lammps-5Jan08/src/pair_coul_debye.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jan08/src/pair_coul_debye.cpp 2008-01-03 18:32:18.000000000 -0700 @@ -0,0 +1,188 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_coul_debye.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; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairCoulDebye::PairCoulDebye(LAMMPS *lmp) : PairCoulCut(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void PairCoulDebye::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,screening; + 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; + int nall = nlocal + atom->nghost; + 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]; + + if (j < nall) factor_coul = 1.0; + else { + factor_coul = special_coul[j/nall]; + j %= nall; + } + + 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; + screening = exp(-kappa*r); + forcecoul = qqrd2e * qtmp*q[j] * screening * (kappa + 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 * qtmp*q[j] * rinv * screening; + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + 0.0,ecoul,fpair,delx,dely,delz); + } + } + } + + if (vflag_fdotr) virial_compute(); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairCoulDebye::settings(int narg, char **arg) +{ + if (narg != 2) error->all("Illegal pair_style command"); + + kappa = atof(arg[0]); + cut_global = atof(arg[1]); + + // reset cutoffs that have been explicitly set + + if (allocated) { + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i+1; j <= atom->ntypes; j++) + if (setflag[i][j]) cut[i][j] = cut_global; + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulDebye::write_restart_settings(FILE *fp) +{ + fwrite(&cut_global,sizeof(double),1,fp); + fwrite(&kappa,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 PairCoulDebye::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_global,sizeof(double),1,fp); + fread(&kappa,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(&kappa,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- */ + +double PairCoulDebye::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,screening; + + r2inv = 1.0/rsq; + r = sqrt(rsq); + rinv = 1.0/r; + screening = exp(-kappa*r); + forcecoul = force->qqrd2e * atom->q[i]*atom->q[j] * + screening * (kappa + rinv); + fforce = factor_coul*forcecoul * r2inv; + + phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * rinv * screening; + return factor_coul*phicoul; +} diff -Naur lammps-5Jan08/src/pair_coul_debye.h lammps-6Jan08/src/pair_coul_debye.h --- lammps-5Jan08/src/pair_coul_debye.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-6Jan08/src/pair_coul_debye.h 2008-01-03 18:32:18.000000000 -0700 @@ -0,0 +1,36 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifndef PAIR_COUL_DEBYE_H +#define PAIR_COUL_DEBYE_H + +#include "pair_coul_cut.h" + +namespace LAMMPS_NS { + +class PairCoulDebye : public PairCoulCut { + public: + PairCoulDebye(class LAMMPS *); + 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 &); + + private: + double kappa; +}; + +} + +#endif diff -Naur lammps-5Jan08/src/pair_lj_cut_coul_cut.h lammps-6Jan08/src/pair_lj_cut_coul_cut.h --- lammps-5Jan08/src/pair_lj_cut_coul_cut.h 2008-01-02 12:24:46.000000000 -0700 +++ lammps-6Jan08/src/pair_lj_cut_coul_cut.h 2008-01-03 18:20:28.000000000 -0700 @@ -27,8 +27,8 @@ void coeff(int, char **); void init_style(); double init_one(int, int); - virtual void write_restart(FILE *); - virtual void read_restart(FILE *); + 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 &); diff -Naur lammps-5Jan08/src/pair_lj_cut_coul_debye.cpp lammps-6Jan08/src/pair_lj_cut_coul_debye.cpp --- lammps-5Jan08/src/pair_lj_cut_coul_debye.cpp 2008-01-02 12:24:46.000000000 -0700 +++ lammps-6Jan08/src/pair_lj_cut_coul_debye.cpp 2008-01-03 18:20:28.000000000 -0700 @@ -160,58 +160,6 @@ proc 0 writes to restart file ------------------------------------------------------------------------- */ -void PairLJCutCoulDebye::write_restart(FILE *fp) -{ - write_restart_settings(fp); - - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut_lj[i][j],sizeof(double),1,fp); - fwrite(&cut_coul[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairLJCutCoulDebye::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(&epsilon[i][j],sizeof(double),1,fp); - fread(&sigma[i][j],sizeof(double),1,fp); - fread(&cut_lj[i][j],sizeof(double),1,fp); - fread(&cut_coul[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - void PairLJCutCoulDebye::write_restart_settings(FILE *fp) { fwrite(&cut_lj_global,sizeof(double),1,fp); diff -Naur lammps-5Jan08/src/pair_lj_cut_coul_debye.h lammps-6Jan08/src/pair_lj_cut_coul_debye.h --- lammps-5Jan08/src/pair_lj_cut_coul_debye.h 2008-01-02 12:24:46.000000000 -0700 +++ lammps-6Jan08/src/pair_lj_cut_coul_debye.h 2008-01-03 18:20:28.000000000 -0700 @@ -23,8 +23,6 @@ PairLJCutCoulDebye(class LAMMPS *); void compute(int, int); void settings(int, char **); - void write_restart(FILE *); - void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); diff -Naur lammps-5Jan08/src/style.h lammps-6Jan08/src/style.h --- lammps-5Jan08/src/style.h 2008-01-03 17:12:32.000000000 -0700 +++ lammps-6Jan08/src/style.h 2008-01-03 18:20:28.000000000 -0700 @@ -277,6 +277,7 @@ #include "pair_buck.h" #include "pair_buck_coul_cut.h" #include "pair_coul_cut.h" +#include "pair_coul_debye.h" #include "pair_hybrid.h" #include "pair_hybrid_overlay.h" #include "pair_lj_cut.h" @@ -294,6 +295,7 @@ PairStyle(buck,PairBuck) PairStyle(buck/coul/cut,PairBuckCoulCut) PairStyle(coul/cut,PairCoulCut) +PairStyle(coul/debye,PairCoulDebye) PairStyle(hybrid,PairHybrid) PairStyle(hybrid/overlay,PairHybridOverlay) PairStyle(lj/cut,PairLJCut)