From 7b7cb297e6beca11414c494e041b5e820dcede77 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Oct 2016 22:26:47 -0400 Subject: [PATCH] add mixing for pair styles gauss and gauss/cut contributed by andrew jewett. also add support for write_data. --- doc/src/pair_gauss.txt | 27 +++++++++++++++-- src/USER-MISC/pair_gauss_cut.cpp | 39 ++++++++++++++++++++++-- src/USER-MISC/pair_gauss_cut.h | 2 ++ src/pair_gauss.cpp | 52 +++++++++++++++++++++++++++----- src/pair_gauss.h | 2 ++ 5 files changed, 110 insertions(+), 12 deletions(-) diff --git a/doc/src/pair_gauss.txt b/doc/src/pair_gauss.txt index 32fc5f6149..8f64122866 100644 --- a/doc/src/pair_gauss.txt +++ b/doc/src/pair_gauss.txt @@ -106,8 +106,31 @@ more instructions on how to use the accelerated styles effectively. [Mixing, shift, table, tail correction, restart, rRESPA info]: -These pair styles do not support mixing. Thus, coefficients for all -I,J pairs must be specified explicitly. +For atom type pairs I,J and I != J, the A, B, H, sigma_h, r_mh +parameters, and the cutoff distance for these pair styles can be mixed: +A (energy units) +sqrt(1/B) (distance units, see below) +H (energy units) +sigma_h (distance units) +r_mh (distance units) +cutoff (distance units):ul + +The default mix value is {geometric}. +Only {arithmetic} and {geometric} mix values are supported. +See the "pair_modify" command for details. + +The A and H parameters are mixed using the same rules normally +used to mix the "epsilon" parameter in a Lennard Jones interaction. +The sigma_h, r_mh, and the cutoff distance are mixed using the same +rules used to mix the "sigma" parameter in a Lennard Jones interaction. +The B parameter is converted to a distance (sigma), before mixing +(using sigma=B^-0.5), and converted back to a coefficient +afterwards (using B=sigma^2). +Negative A values are converted to positive A values (using abs(A)) +before mixing, and converted back after mixing +(by multiplying by sign(Ai)*sign(Aj)). +This way, if either particle is repulsive (if Ai<0 or Aj<0), +then the default interaction between both particles will be repulsive. The {gauss} style does not support the "pair_modify"_pair_modify.html shift option. There is no effect due to the Gaussian well beyond the diff --git a/src/USER-MISC/pair_gauss_cut.cpp b/src/USER-MISC/pair_gauss_cut.cpp index 0df6734849..f9cdcb96e6 100644 --- a/src/USER-MISC/pair_gauss_cut.cpp +++ b/src/USER-MISC/pair_gauss_cut.cpp @@ -39,6 +39,7 @@ using namespace MathConst; PairGaussCut::PairGaussCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 0; + writedata = 1; } /* ---------------------------------------------------------------------- */ @@ -96,7 +97,6 @@ void PairGaussCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; @@ -221,7 +221,21 @@ void PairGaussCut::coeff(int narg, char **arg) double PairGaussCut::init_one(int i, int j) { - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + if (setflag[i][j] == 0) { + hgauss[i][j] = mix_energy(fabs(hgauss[i][i]), fabs(hgauss[j][j]), + fabs(sigmah[i][i]), fabs(sigmah[j][j])); + + // If either of the particles is repulsive (ie, if hgauss > 0), + // then the interaction between both is repulsive. + double sign_hi = (hgauss[i][i] >= 0.0) ? 1.0 : -1.0; + double sign_hj = (hgauss[j][j] >= 0.0) ? 1.0 : -1.0; + hgauss[i][j] *= MAX(sign_hi, sign_hj); + + sigmah[i][j] = mix_distance(sigmah[i][i], sigmah[j][j]); + rmh[i][j] = mix_distance(rmh[i][i], rmh[j][j]); + cut[i][j] = mix_distance(cut[i][i], cut[j][j]); + } + pgauss[i][j] = hgauss[i][j] / sqrt(MY_2PI) / sigmah[i][j]; if (offset_flag) { @@ -334,6 +348,27 @@ void PairGaussCut::read_restart_settings(FILE *fp) MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairGaussCut::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g %g\n",i,hgauss[i][i],rmh[i][i],sigmah[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairGaussCut::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g %g\n",i,j,hgauss[i][j],rmh[i][j],sigmah[i][j],cut[i][j]); +} + /* ---------------------------------------------------------------------- */ double PairGaussCut::single(int i, int j, int itype, int jtype, double rsq, diff --git a/src/USER-MISC/pair_gauss_cut.h b/src/USER-MISC/pair_gauss_cut.h index b43276fe7d..9173f83bc0 100644 --- a/src/USER-MISC/pair_gauss_cut.h +++ b/src/USER-MISC/pair_gauss_cut.h @@ -42,6 +42,8 @@ class PairGaussCut : public Pair { virtual void read_restart(FILE *); virtual void write_restart_settings(FILE *); virtual void read_restart_settings(FILE *); + virtual void write_data(FILE *fp); + virtual void write_data_all(FILE *fp); virtual double memory_usage(); diff --git a/src/pair_gauss.cpp b/src/pair_gauss.cpp index 58373021e2..1b922a215f 100644 --- a/src/pair_gauss.cpp +++ b/src/pair_gauss.cpp @@ -34,10 +34,11 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairGauss::PairGauss(LAMMPS *lmp) :Pair(lmp) +PairGauss::PairGauss(LAMMPS *lmp) : Pair(lmp) { nextra = 1; pvector = new double[1]; + writedata = 1; } /* ---------------------------------------------------------------------- */ @@ -187,7 +188,7 @@ void PairGauss::coeff(int narg, char **arg) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); - int ilo, ihi, jlo, jhi; + int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); @@ -204,7 +205,7 @@ void PairGauss::coeff(int narg, char **arg) b[i][j] = b_one; cut[i][j] = cut_one; setflag[i][j] = 1; - count++ ; + count++; } } @@ -217,12 +218,27 @@ void PairGauss::coeff(int narg, char **arg) double PairGauss::init_one(int i, int j) { + if (setflag[i][j] == 0) { + double sign_bi = (b[i][i] >= 0.0) ? 1.0 : -1.0; + double sign_bj = (b[j][j] >= 0.0) ? 1.0 : -1.0; + double si = sqrt(0.5/fabs(b[i][i])); + double sj = sqrt(0.5/fabs(b[j][j])); + double sij = mix_distance(si, sj); + b[i][j] = 0.5 / (sij*sij); + b[i][j] *= MAX(sign_bi, sign_bj); - // This error is triggered when ti is performed on lj/cut tail - // in presence of extra atom type for tether sites - // "i = 2 j = 1 ERROR: All pair coeffs are not set (pair_gauss.cpp:223)" - // if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + // Negative "a" values are useful for simulating repulsive particles. + // If either of the particles is repulsive (a<0), then the + // interaction between both is repulsive. + double sign_ai = (a[i][i] >= 0.0) ? 1.0 : -1.0; + double sign_aj = (a[j][j] >= 0.0) ? 1.0 : -1.0; + a[i][j] = mix_energy(fabs(a[i][i]), fabs(a[j][j]), si, sj); + a[i][j] *= MIN(sign_ai, sign_aj); + cut[i][j] = mix_distance(cut[i][i],cut[j][j]); + } + + // cutoff correction to energy if (offset_flag) offset[i][j] = a[i][j]*exp(-b[i][j]*cut[i][j]*cut[i][j]); else offset[i][j] = 0.0; @@ -260,7 +276,6 @@ void PairGauss::write_restart(FILE *fp) void PairGauss::read_restart(FILE *fp) { read_restart_settings(fp); - allocate(); int i,j; @@ -309,6 +324,27 @@ void PairGauss::read_restart_settings(FILE *fp) MPI_Bcast(&mix_flag,1,MPI_INT,0,world); } +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void PairGauss::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + fprintf(fp,"%d %g %g\n",i,a[i][i],b[i][i]); +} + +/* ---------------------------------------------------------------------- + proc 0 writes all pairs to data file +------------------------------------------------------------------------- */ + +void PairGauss::write_data_all(FILE *fp) +{ + for (int i = 1; i <= atom->ntypes; i++) + for (int j = i; j <= atom->ntypes; j++) + fprintf(fp,"%d %d %g %g %g\n",i,j,a[i][j],b[i][j],cut[i][j]); +} + /* ---------------------------------------------------------------------- */ double PairGauss::single(int i, int j, int itype, int jtype, double rsq, diff --git a/src/pair_gauss.h b/src/pair_gauss.h index 54ba29214e..901658a4c2 100644 --- a/src/pair_gauss.h +++ b/src/pair_gauss.h @@ -36,6 +36,8 @@ class PairGauss : public Pair { void read_restart(FILE *); void write_restart_settings(FILE *); void read_restart_settings(FILE *); + void write_data(FILE *fp); + void write_data_all(FILE *fp); double single(int, int, int, int, double, double, double, double &); void *extract(const char *, int &);