diff --git a/src/KSPACE/Install.csh b/src/KSPACE/Install.csh index 5646081304..4e5569c155 100644 --- a/src/KSPACE/Install.csh +++ b/src/KSPACE/Install.csh @@ -8,6 +8,7 @@ if ($1 == 1) then cp pppm.cpp .. cp pppm_tip4p.cpp .. cp pair_buck_coul_long.cpp .. + cp pair_coul_long.cpp .. cp pair_lj_cut_coul_long.cpp .. cp pair_lj_cut_coul_long_tip4p.cpp .. cp pair_lj_charmm_coul_long.cpp .. @@ -20,6 +21,7 @@ if ($1 == 1) then cp pppm.h .. cp pppm_tip4p.h .. cp pair_buck_coul_long.h .. + cp pair_coul_long.h .. cp pair_lj_cut_coul_long.h .. cp pair_lj_cut_coul_long_tip4p.h .. cp pair_lj_charmm_coul_long.h .. @@ -37,6 +39,7 @@ else if ($1 == 0) then rm ../pppm.cpp rm ../pppm_tip4p.cpp rm ../pair_buck_coul_long.cpp + rm ../pair_coul_long.cpp rm ../pair_lj_cut_coul_long.cpp rm ../pair_lj_cut_coul_long_tip4p.cpp rm ../pair_lj_charmm_coul_long.cpp @@ -49,6 +52,7 @@ else if ($1 == 0) then rm ../pppm.h rm ../pppm_tip4p.h rm ../pair_buck_coul_long.h + rm ../pair_coul_long.h rm ../pair_lj_cut_coul_long.h rm ../pair_lj_cut_coul_long_tip4p.h rm ../pair_lj_charmm_coul_long.h diff --git a/src/KSPACE/pair_buck_coul_long.cpp b/src/KSPACE/pair_buck_coul_long.cpp index fabc8c1f1b..a2e998a187 100644 --- a/src/KSPACE/pair_buck_coul_long.cpp +++ b/src/KSPACE/pair_buck_coul_long.cpp @@ -20,9 +20,9 @@ #include "comm.h" #include "force.h" #include "kspace.h" -#include "update.h" -#include "memory.h" #include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -66,20 +66,18 @@ PairBuckCoulLong::~PairBuckCoulLong() void PairBuckCoulLong::compute(int eflag, int vflag) { - int i,j,k,numneigh,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcecoul,forcebuck,fforce,factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double factor,phicoul,phibuck,r,rexp; - int *neighs; - double **f; + int *ilist,*jlist,*numneigh,**firstneigh; eng_vdwl = eng_coul = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag == 2) f = update->f_pair; - else f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -89,19 +87,25 @@ void PairBuckCoulLong::compute(int eflag, int vflag) 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 (i = 0; i < nlocal; i++) { + 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]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -314,6 +318,8 @@ void PairBuckCoulLong::init_style() else if (strcmp(force->kspace_style,"pppm") == 0) g_ewald = force->kspace->g_ewald; else error->all("Pair style is incompatible with KSpace style"); + + int irequest = neighbor->request(this); } /* ---------------------------------------------------------------------- diff --git a/src/KSPACE/pair_buck_coul_long.h b/src/KSPACE/pair_buck_coul_long.h index 5b3886efd8..410ea0640d 100644 --- a/src/KSPACE/pair_buck_coul_long.h +++ b/src/KSPACE/pair_buck_coul_long.h @@ -25,8 +25,8 @@ class PairBuckCoulLong : public Pair { void compute(int, int); void settings(int, char **); void coeff(int, char **); - double init_one(int, int); void init_style(); + double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); diff --git a/src/KSPACE/pair_coul_long.cpp b/src/KSPACE/pair_coul_long.cpp new file mode 100644 index 0000000000..84c9b14c56 --- /dev/null +++ b/src/KSPACE/pair_coul_long.cpp @@ -0,0 +1,578 @@ +/* ---------------------------------------------------------------------- + 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: Paul Crozier (SNL) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "pair_coul_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 MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +#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 + +/* ---------------------------------------------------------------------- */ + +PairCoulLong::PairCoulLong(LAMMPS *lmp) : Pair(lmp) +{ + ftable = NULL; +} + +/* ---------------------------------------------------------------------- */ + +PairCoulLong::~PairCoulLong() +{ + if (allocated) { + memory->destroy_2d_int_array(setflag); + memory->destroy_2d_double_array(cutsq); + } + if (ftable) free_tables(); +} + +/* ---------------------------------------------------------------------- */ + +void PairCoulLong::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itable; + double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fraction,table; + double r,r2inv,forcecoul,fforce,factor_coul; + double grij,expm2,prefactor,t,erfc; + double factor,phicoul; + int *ilist,*jlist,*numneigh,**firstneigh; + float rsq; + int *int_rsq = (int *) &rsq; + + eng_coul = 0.0; + if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; + + double **x = atom->x; + double **f = atom->f; + double *q = atom->q; + int nlocal = atom->nlocal; + int nall = atom->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]; + 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; + + if (rsq < cut_coulsq) { + r2inv = 1.0/rsq; + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrtf(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; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + itable = *int_rsq & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } + + fforce = forcecoul * r2inv; + + f[i][0] += delx*fforce; + f[i][1] += dely*fforce; + f[i][2] += delz*fforce; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fforce; + f[j][1] -= dely*fforce; + f[j][2] -= delz*fforce; + } + + if (eflag) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = qtmp*q[j] * table; + } + if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; + if (newton_pair || j < nlocal) eng_coul += phicoul; + else eng_coul += 0.5*phicoul; + } + + if (vflag == 1) { + if (newton_pair == 0 && j >= nlocal) fforce *= 0.5; + virial[0] += delx*delx*fforce; + virial[1] += dely*dely*fforce; + virial[2] += delz*delz*fforce; + virial[3] += delx*dely*fforce; + virial[4] += delx*delz*fforce; + virial[5] += dely*delz*fforce; + } + } + } + } + if (vflag == 2) virial_compute(); +} +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairCoulLong::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairCoulLong::settings(int narg, char **arg) +{ + if (narg != 1) error->all("Illegal pair_style command"); + + cut_coul = atof(arg[0]); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairCoulLong::coeff(int narg, char **arg) +{ + if (narg != 2) error->all("Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(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++) { + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairCoulLong::init_style() +{ + int i,j; + + if (!atom->q_flag) + error->all("Pair style lj/cut/coul/long requires atom attribute q"); + + int irequest = neighbor->request(this); + + cut_coulsq = cut_coul * cut_coul; + + // set & error check interior rRESPA cutoffs + + if (strcmp(update->integrate_style,"respa") == 0) { + if (((Respa *) update->integrate)->level_inner >= 0) { + cut_respa = ((Respa *) update->integrate)->cutoff; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) + if (cut_coul < cut_respa[3]) + error->all("Pair cutoff < Respa interior cutoff"); + } + } else cut_respa = NULL; + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all("Pair style is incompatible with KSpace style"); + else if (strcmp(force->kspace_style,"ewald") == 0) + g_ewald = force->kspace->g_ewald; + else if (strcmp(force->kspace_style,"pppm") == 0) + g_ewald = force->kspace->g_ewald; + else error->all("Pair style is incompatible with KSpace style"); + + // setup force tables + + if (ncoultablebits) init_tables(); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairCoulLong::init_one(int i, int j) +{ + return cut_coul; +} + +/* ---------------------------------------------------------------------- + setup force tables used in compute routines +------------------------------------------------------------------------- */ + +void PairCoulLong::init_tables() +{ + int masklo,maskhi; + double r,grij,expm2,derfc,rsw; + double qqrd2e = force->qqrd2e; + + tabinnersq = tabinner*tabinner; + init_bitmap(tabinner,cut_coul,ncoultablebits, + masklo,maskhi,ncoulmask,ncoulshiftbits); + + int ntable = 1; + for (int i = 0; i < ncoultablebits; i++) ntable *= 2; + + // linear lookup tables of length N = 2^ncoultablebits + // stored value = value at lower edge of bin + // d values = delta from lower edge to upper edge of bin + + if (ftable) free_tables(); + + rtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:rtable"); + ftable = (double *) memory->smalloc(ntable*sizeof(double),"pair:ftable"); + ctable = (double *) memory->smalloc(ntable*sizeof(double),"pair:ctable"); + etable = (double *) memory->smalloc(ntable*sizeof(double),"pair:etable"); + drtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:drtable"); + dftable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dftable"); + dctable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dctable"); + detable = (double *) memory->smalloc(ntable*sizeof(double),"pair:detable"); + + if (cut_respa == NULL) { + vtable = ptable = dvtable = dptable = NULL; + } else { + vtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:vtable"); + ptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:ptable"); + dvtable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dvtable"); + dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable"); + } + + float rsq; + int *int_rsq = (int *) &rsq; + float minrsq; + int *int_minrsq = (int *) &minrsq; + int itablemin; + *int_minrsq = 0 << ncoulshiftbits; + *int_minrsq = *int_minrsq | maskhi; + + for (int i = 0; i < ntable; i++) { + *int_rsq = i << ncoulshiftbits; + *int_rsq = *int_rsq | masklo; + if (rsq < tabinnersq) { + *int_rsq = i << ncoulshiftbits; + *int_rsq = *int_rsq | maskhi; + } + r = sqrtf(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + derfc = erfc(grij); + if (cut_respa == NULL) { + rtable[i] = rsq; + ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + ctable[i] = qqrd2e/r; + etable[i] = qqrd2e/r * derfc; + } else { + rtable[i] = rsq; + ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); + ctable[i] = 0.0; + etable[i] = qqrd2e/r * derfc; + ptable[i] = qqrd2e/r; + vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + if (rsq > cut_respa[2]*cut_respa[2]) { + if (rsq < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + ctable[i] = qqrd2e/r; + } + } + } + minrsq = MIN(minrsq,rsq); + } + + tabinnersq = minrsq; + + int ntablem1 = ntable - 1; + + for (int i = 0; i < ntablem1; i++) { + drtable[i] = 1.0/(rtable[i+1] - rtable[i]); + dftable[i] = ftable[i+1] - ftable[i]; + dctable[i] = ctable[i+1] - ctable[i]; + detable[i] = etable[i+1] - etable[i]; + } + + if (cut_respa) { + for (int i = 0; i < ntablem1; i++) { + dvtable[i] = vtable[i+1] - vtable[i]; + dptable[i] = ptable[i+1] - ptable[i]; + } + } + + // get the delta values for the last table entries + // tables are connected periodically between 0 and ntablem1 + + drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); + dftable[ntablem1] = ftable[0] - ftable[ntablem1]; + dctable[ntablem1] = ctable[0] - ctable[ntablem1]; + detable[ntablem1] = etable[0] - etable[ntablem1]; + if (cut_respa) { + dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; + dptable[ntablem1] = ptable[0] - ptable[ntablem1]; + } + + // get the correct delta values at itablemax + // smallest r is in bin itablemin + // largest r is in bin itablemax, which is itablemin-1, + // or ntablem1 if itablemin=0 + // deltas at itablemax only needed if corresponding rsq < cut*cut + // if so, compute deltas between rsq and cut*cut + + double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; + itablemin = *int_minrsq & ncoulmask; + itablemin >>= ncoulshiftbits; + int itablemax = itablemin - 1; + if (itablemin == 0) itablemax = ntablem1; + *int_rsq = itablemax << ncoulshiftbits; + *int_rsq = *int_rsq | maskhi; + + if (rsq < cut_coulsq) { + rsq = cut_coulsq; + r = sqrtf(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + derfc = erfc(grij); + + if (cut_respa == NULL) { + f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + c_tmp = qqrd2e/r; + e_tmp = qqrd2e/r * derfc; + } else { + f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); + c_tmp = 0.0; + e_tmp = qqrd2e/r * derfc; + p_tmp = qqrd2e/r; + v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + if (rsq > cut_respa[2]*cut_respa[2]) { + if (rsq < cut_respa[3]*cut_respa[3]) { + rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); + f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); + } else { + f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); + c_tmp = qqrd2e/r; + } + } + } + + drtable[itablemax] = 1.0/(rsq - rtable[itablemax]); + dftable[itablemax] = f_tmp - ftable[itablemax]; + dctable[itablemax] = c_tmp - ctable[itablemax]; + detable[itablemax] = e_tmp - etable[itablemax]; + if (cut_respa) { + dvtable[itablemax] = v_tmp - vtable[itablemax]; + dptable[itablemax] = p_tmp - ptable[itablemax]; + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulLong::write_restart(FILE *fp) +{ + write_restart_settings(fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairCoulLong::read_restart(FILE *fp) +{ + read_restart_settings(fp); + + allocate(); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairCoulLong::write_restart_settings(FILE *fp) +{ + fwrite(&cut_coul,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 PairCoulLong::read_restart_settings(FILE *fp) +{ + if (comm->me == 0) { + fread(&cut_coul,sizeof(double),1,fp); + fread(&offset_flag,sizeof(int),1,fp); + fread(&mix_flag,sizeof(int),1,fp); + } + MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); + MPI_Bcast(&offset_flag,1,MPI_INT,0,world); + MPI_Bcast(&mix_flag,1,MPI_INT,0,world); +} + +/* ---------------------------------------------------------------------- + free memory for tables used in pair computations +------------------------------------------------------------------------- */ + +void PairCoulLong::free_tables() +{ + memory->sfree(rtable); + memory->sfree(drtable); + memory->sfree(ftable); + memory->sfree(dftable); + memory->sfree(ctable); + memory->sfree(dctable); + memory->sfree(etable); + memory->sfree(detable); + memory->sfree(vtable); + memory->sfree(dvtable); + memory->sfree(ptable); + memory->sfree(dptable); +} + +/* ---------------------------------------------------------------------- */ + +void PairCoulLong::single(int i, int j, int itype, int jtype, + double rsq, + double factor_coul, double factor_lj, + int eflag, One &one) +{ + double r2inv,r,grij,expm2,t,erfc,prefactor; + double fraction,table,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; + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + float rsq_single = rsq; + int *int_rsq = (int *) &rsq_single; + itable = *int_rsq & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_single - 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; + } + } + one.fforce = forcecoul * r2inv; + + if (eflag) { + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc; + 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; + one.eng_coul = phicoul; + one.eng_vdwl = 0.0; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairCoulLong::extract_long(double *p_cut_coul) +{ + *p_cut_coul = cut_coul; +} diff --git a/src/KSPACE/pair_coul_long.h b/src/KSPACE/pair_coul_long.h new file mode 100644 index 0000000000..f11cb11271 --- /dev/null +++ b/src/KSPACE/pair_coul_long.h @@ -0,0 +1,55 @@ +/* ---------------------------------------------------------------------- + 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_LONG_H +#define PAIR_COUL_LONG_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairCoulLong : public Pair { + public: + PairCoulLong(class LAMMPS *); + ~PairCoulLong(); + void compute(int, int); + 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 *); + void single(int, int, int, int, double, double, double, int, One &); + + void extract_long(double *); + + private: + double cut_coul,cut_coulsq; + double *cut_respa; + double g_ewald; + + double tabinnersq; + double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable; + double *etable,*detable,*ptable,*dptable,*vtable,*dvtable; + int ncoulshiftbits,ncoulmask; + + void allocate(); + void init_tables(); + void free_tables(); +}; + +} + +#endif diff --git a/src/KSPACE/pair_lj_charmm_coul_long.cpp b/src/KSPACE/pair_lj_charmm_coul_long.cpp index bc158fef0b..ffc0c86fae 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_long.cpp @@ -27,8 +27,10 @@ #include "update.h" #include "integrate.h" #include "respa.h" -#include "memory.h" #include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -80,22 +82,20 @@ PairLJCharmmCoulLong::~PairLJCharmmCoulLong() void PairLJCharmmCoulLong::compute(int eflag, int vflag) { - int i,j,k,numneigh,itype,jtype,itable; + int i,j,ii,jj,inum,jnum,itype,jtype,itable; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fraction,table; double r,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double factor,phicoul,philj,switch1,switch2; - int *neighs; - double **f; + int *ilist,*jlist,*numneigh,**firstneigh; float rsq; int *int_rsq = (int *) &rsq; eng_vdwl = eng_coul = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag == 2) f = update->f_pair; - else f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -105,19 +105,25 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag) 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 (i = 0; i < nlocal; i++) { + 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]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -226,14 +232,14 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag) void PairLJCharmmCoulLong::compute_inner() { - int i,j,k,numneigh,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double rsw; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; - double **f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -243,6 +249,11 @@ void PairLJCharmmCoulLong::compute_inner() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; + inum = listinner->inum; + ilist = listinner->ilist; + numneigh = listinner->numneigh; + firstneigh = listinner->firstneigh; + double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -252,17 +263,18 @@ void PairLJCharmmCoulLong::compute_inner() // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { + 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]; - neighs = neighbor->firstneigh_inner[i]; - numneigh = neighbor->numneigh_inner[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -309,15 +321,15 @@ void PairLJCharmmCoulLong::compute_inner() void PairLJCharmmCoulLong::compute_middle() { - int i,j,k,numneigh,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double philj,switch1,switch2; double rsw; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; - double **f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -327,6 +339,11 @@ void PairLJCharmmCoulLong::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; + inum = listmiddle->inum; + ilist = listmiddle->ilist; + numneigh = listmiddle->numneigh; + firstneigh = listmiddle->firstneigh; + double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; double cut_out_on = cut_respa[2]; @@ -341,18 +358,19 @@ void PairLJCharmmCoulLong::compute_middle() // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { - + 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]; - neighs = neighbor->firstneigh_middle[i]; - numneigh = neighbor->numneigh_middle[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -410,21 +428,21 @@ void PairLJCharmmCoulLong::compute_middle() void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) { - int i,j,k,numneigh,itype,jtype,itable; + int i,j,ii,jj,inum,jnum,itype,jtype,itable; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fraction,table; double r,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double factor,phicoul,philj,switch1,switch2; double rsw; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; float rsq; int *int_rsq = (int *) &rsq; eng_vdwl = eng_coul = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - double **f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -434,6 +452,11 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; + inum = listouter->inum; + ilist = listouter->ilist; + numneigh = listouter->numneigh; + firstneigh = listouter->firstneigh; + double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -443,18 +466,18 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { - + 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]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; - - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -703,44 +726,6 @@ void PairLJCharmmCoulLong::coeff(int narg, char **arg) if (count == 0) error->all("Incorrect args for pair coefficients"); } -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairLJCharmmCoulLong::init_one(int i, int j) -{ - // always mix arithmetically - - if (setflag[i][j] == 0) { - epsilon[i][j] = sqrt(epsilon[i][i]*epsilon[j][j]); - sigma[i][j] = 0.5 * (sigma[i][i] + sigma[j][j]); - eps14[i][j] = sqrt(eps14[i][i]*eps14[j][j]); - sigma14[i][j] = 0.5 * (sigma14[i][i] + sigma14[j][j]); - } - - double cut = MAX(cut_lj,cut_coul); - - lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj14_1[i][j] = 48.0 * eps14[i][j] * pow(sigma14[i][j],12.0); - lj14_2[i][j] = 24.0 * eps14[i][j] * pow(sigma14[i][j],6.0); - lj14_3[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],12.0); - lj14_4[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],6.0); - - lj1[j][i] = lj1[i][j]; - lj2[j][i] = lj2[i][j]; - lj3[j][i] = lj3[i][j]; - lj4[j][i] = lj4[i][j]; - lj14_1[j][i] = lj14_1[i][j]; - lj14_2[j][i] = lj14_2[i][j]; - lj14_3[j][i] = lj14_3[i][j]; - lj14_4[j][i] = lj14_4[i][j]; - - return cut; -} - /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ @@ -750,6 +735,42 @@ void PairLJCharmmCoulLong::init_style() if (!atom->q_flag) error->all("Pair style lj/charmm/coul/long requires atom attribute q"); + // request regular or rRESPA neighbor lists + + int irequest; + + if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) { + int respa = 0; + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + + if (respa == 0) irequest = neighbor->request(this); + else if (respa == 1) { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } else { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 2; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respamiddle = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } + + } else irequest = neighbor->request(this); + // require cut_lj_inner < cut_lj if (cut_lj_inner >= cut_lj) @@ -793,6 +814,57 @@ void PairLJCharmmCoulLong::init_style() if (ncoultablebits) init_tables(); } +/* ---------------------------------------------------------------------- + neighbor callback to inform pair style of neighbor list to use + regular or rRESPA +------------------------------------------------------------------------- */ + +void PairLJCharmmCoulLong::init_list(int id, NeighList *ptr) +{ + if (id == 0) list = ptr; + else if (id == 1) listinner = ptr; + else if (id == 2) listmiddle = ptr; + else if (id == 3) listouter = ptr; +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairLJCharmmCoulLong::init_one(int i, int j) +{ + // always mix arithmetically + + if (setflag[i][j] == 0) { + epsilon[i][j] = sqrt(epsilon[i][i]*epsilon[j][j]); + sigma[i][j] = 0.5 * (sigma[i][i] + sigma[j][j]); + eps14[i][j] = sqrt(eps14[i][i]*eps14[j][j]); + sigma14[i][j] = 0.5 * (sigma14[i][i] + sigma14[j][j]); + } + + double cut = MAX(cut_lj,cut_coul); + + lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); + lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); + lj14_1[i][j] = 48.0 * eps14[i][j] * pow(sigma14[i][j],12.0); + lj14_2[i][j] = 24.0 * eps14[i][j] * pow(sigma14[i][j],6.0); + lj14_3[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],12.0); + lj14_4[i][j] = 4.0 * eps14[i][j] * pow(sigma14[i][j],6.0); + + lj1[j][i] = lj1[i][j]; + lj2[j][i] = lj2[i][j]; + lj3[j][i] = lj3[i][j]; + lj4[j][i] = lj4[i][j]; + lj14_1[j][i] = lj14_1[i][j]; + lj14_2[j][i] = lj14_2[i][j]; + lj14_3[j][i] = lj14_3[i][j]; + lj14_4[j][i] = lj14_4[i][j]; + + return cut; +} + /* ---------------------------------------------------------------------- setup force tables used in compute routines ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_charmm_coul_long.h b/src/KSPACE/pair_lj_charmm_coul_long.h index 13dc9526a9..515c8ee765 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.h +++ b/src/KSPACE/pair_lj_charmm_coul_long.h @@ -25,8 +25,9 @@ class PairLJCharmmCoulLong : public Pair { void compute(int, int); void settings(int, char **); void coeff(int, char **); - double init_one(int, int); void init_style(); + void init_list(int, class NeighList *); + double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index 3c0e008fc6..c80f4ad877 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -27,8 +27,10 @@ #include "update.h" #include "integrate.h" #include "respa.h" -#include "memory.h" #include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -77,22 +79,20 @@ PairLJCutCoulLong::~PairLJCutCoulLong() void PairLJCutCoulLong::compute(int eflag, int vflag) { - int i,j,k,numneigh,itype,jtype,itable; + int i,j,ii,jj,inum,jnum,itype,jtype,itable; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fraction,table; double r,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double factor,phicoul,philj; - int *neighs; - double **f; + int *ilist,*jlist,*numneigh,**firstneigh; float rsq; int *int_rsq = (int *) &rsq; eng_vdwl = eng_coul = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag == 2) f = update->f_pair; - else f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -102,19 +102,25 @@ void PairLJCutCoulLong::compute(int eflag, int vflag) 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 (i = 0; i < nlocal; i++) { + 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]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -211,14 +217,14 @@ void PairLJCutCoulLong::compute(int eflag, int vflag) void PairLJCutCoulLong::compute_inner() { - int i,j,k,numneigh,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double rsw; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; - double **f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -228,6 +234,11 @@ void PairLJCutCoulLong::compute_inner() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; + inum = listinner->inum; + ilist = listinner->ilist; + numneigh = listinner->numneigh; + firstneigh = listinner->firstneigh; + double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; @@ -237,17 +248,18 @@ void PairLJCutCoulLong::compute_inner() // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { + 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]; - neighs = neighbor->firstneigh_inner[i]; - numneigh = neighbor->numneigh_inner[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -295,14 +307,14 @@ void PairLJCutCoulLong::compute_inner() void PairLJCutCoulLong::compute_middle() { - int i,j,k,numneigh,itype,jtype; + int i,j,ii,jj,inum,jnum,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz; double rsq,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double rsw; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; - double **f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -312,6 +324,11 @@ void PairLJCutCoulLong::compute_middle() int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; + inum = listmiddle->inum; + ilist = listmiddle->ilist; + numneigh = listmiddle->numneigh; + firstneigh = listmiddle->firstneigh; + double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; double cut_out_on = cut_respa[2]; @@ -326,17 +343,18 @@ void PairLJCutCoulLong::compute_middle() // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { + 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]; - neighs = neighbor->firstneigh_middle[i]; - numneigh = neighbor->numneigh_middle[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -388,21 +406,21 @@ void PairLJCutCoulLong::compute_middle() void PairLJCutCoulLong::compute_outer(int eflag, int vflag) { - int i,j,k,numneigh,itype,jtype,itable; + int i,j,ii,jj,inum,jnum,itype,jtype,itable; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fraction,table; double r,r2inv,r6inv,forcecoul,forcelj,fforce,factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double factor,phicoul,philj; double rsw; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; float rsq; int *int_rsq = (int *) &rsq; eng_vdwl = eng_coul = 0.0; if (vflag) for (i = 0; i < 6; i++) virial[i] = 0.0; - double **f = atom->f; double **x = atom->x; + double **f = atom->f; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; @@ -412,6 +430,11 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; + inum = listouter->inum; + ilist = listouter->ilist; + numneigh = listouter->numneigh; + firstneigh = listouter->firstneigh; + double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; @@ -421,18 +444,18 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag) // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { - + 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]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -650,6 +673,95 @@ void PairLJCutCoulLong::coeff(int narg, char **arg) if (count == 0) error->all("Incorrect args for pair coefficients"); } +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::init_style() +{ + int i,j; + + if (!atom->q_flag) + error->all("Pair style lj/cut/coul/long requires atom attribute q"); + + // request regular or rRESPA neighbor lists + + int irequest; + + if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) { + int respa = 0; + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + + if (respa == 0) irequest = neighbor->request(this); + else if (respa == 1) { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } else { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 2; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respamiddle = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } + + } else irequest = neighbor->request(this); + + cut_coulsq = cut_coul * cut_coul; + + // set & error check interior rRESPA cutoffs + + if (strcmp(update->integrate_style,"respa") == 0) { + if (((Respa *) update->integrate)->level_inner >= 0) { + cut_respa = ((Respa *) update->integrate)->cutoff; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) + if (MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) + error->all("Pair cutoff < Respa interior cutoff"); + } + } else cut_respa = NULL; + + // insure use of KSpace long-range solver, set g_ewald + + if (force->kspace == NULL) + error->all("Pair style is incompatible with KSpace style"); + else if (strcmp(force->kspace_style,"ewald") == 0) + g_ewald = force->kspace->g_ewald; + else if (strcmp(force->kspace_style,"pppm") == 0) + g_ewald = force->kspace->g_ewald; + else error->all("Pair style is incompatible with KSpace style"); + + // setup force tables + + if (ncoultablebits) init_tables(); +} + +/* ---------------------------------------------------------------------- + neighbor callback to inform pair style of neighbor list to use + regular or rRESPA +------------------------------------------------------------------------- */ + +void PairLJCutCoulLong::init_list(int id, NeighList *ptr) +{ + if (id == 0) list = ptr; + else if (id == 1) listinner = ptr; + else if (id == 2) listmiddle = ptr; + else if (id == 3) listouter = ptr; +} + /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ @@ -713,46 +825,6 @@ double PairLJCutCoulLong::init_one(int i, int j) return cut; } -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::init_style() -{ - int i,j; - - if (!atom->q_flag) - error->all("Pair style lj/cut/coul/long requires atom attribute q"); - - cut_coulsq = cut_coul * cut_coul; - - // set & error check interior rRESPA cutoffs - - if (strcmp(update->integrate_style,"respa") == 0) { - if (((Respa *) update->integrate)->level_inner >= 0) { - cut_respa = ((Respa *) update->integrate)->cutoff; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - if (MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all("Pair cutoff < Respa interior cutoff"); - } - } else cut_respa = NULL; - - // insure use of KSpace long-range solver, set g_ewald - - if (force->kspace == NULL) - error->all("Pair style is incompatible with KSpace style"); - else if (strcmp(force->kspace_style,"ewald") == 0) - g_ewald = force->kspace->g_ewald; - else if (strcmp(force->kspace_style,"pppm") == 0) - g_ewald = force->kspace->g_ewald; - else error->all("Pair style is incompatible with KSpace style"); - - // setup force tables - - if (ncoultablebits) init_tables(); -} - /* ---------------------------------------------------------------------- setup force tables used in compute routines ------------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_cut_coul_long.h b/src/KSPACE/pair_lj_cut_coul_long.h index aaae3b8afb..79ab24b88f 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.h +++ b/src/KSPACE/pair_lj_cut_coul_long.h @@ -25,8 +25,9 @@ class PairLJCutCoulLong : public Pair { virtual void compute(int, int); virtual void settings(int, char **); void coeff(int, char **); - double init_one(int, int); virtual void init_style(); + void init_list(int, class NeighList *); + double init_one(int, int); void write_restart(FILE *); void read_restart(FILE *); virtual void write_restart_settings(FILE *); diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index 1a6a0e7f02..0a18c6410f 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -29,8 +29,10 @@ #include "kspace.h" #include "update.h" #include "respa.h" -#include "memory.h" #include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" #include "error.h" using namespace LAMMPS_NS; @@ -58,7 +60,7 @@ PairLJCutCoulLongTIP4P::PairLJCutCoulLongTIP4P(LAMMPS *lmp) : void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) { - int i,j,k,numneigh,itype,jtype,itable; + int i,j,ii,jj,inum,jnum,itype,jtype,itable; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,fraction,table; double delx1,dely1,delz1,delx2,dely2,delz2,delx3,dely3,delz3; double r,r2inv,r6inv,forcecoul,forcelj,cforce,negforce; @@ -69,7 +71,7 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) double xiM[3],xjM[3]; double *x1,*x2; double fO[3],fH[3]; - int *neighs; + int *ilist,*jlist,*numneigh,**firstneigh; double **f; float rsq; int *int_rsq = (int *) &rsq; @@ -78,7 +80,8 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) if (vflag) for (i = 0; i < 6; i++) virial[i] = tvirial[i] = 0.0; if (vflag == 2) { - f = update->f_pair; + //f = update->f_pair; + f = atom->f; tf = atom->f; } else f = atom->f; @@ -91,9 +94,15 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) double *special_lj = force->special_lj; double qqrd2e = force->qqrd2e; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + // loop over neighbors of my atoms - for (i = 0; i < nlocal; i++) { + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; qtmp = q[i]; xtmp = x[i][0]; ytmp = x[i][1]; @@ -103,11 +112,11 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) find_M(i,iH1,iH2,xiM); x1 = xiM; } else x1 = x[i]; - neighs = neighbor->firstneigh[i]; - numneigh = neighbor->numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (k = 0; k < numneigh; k++) { - j = neighs[k]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; if (j < nall) factor_coul = factor_lj = 1.0; else { @@ -281,12 +290,12 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) tf[j][1] -= dely * cforce; tf[j][2] -= delz * cforce; - tvirial[0] += 0.5 * (delx * delx * cforce); - tvirial[1] += 0.5 * (dely * dely * cforce); - tvirial[2] += 0.5 * (delz * delz * cforce); - tvirial[3] += 0.5 * (dely * delx * cforce); - tvirial[4] += 0.5 * (delz * delx * cforce); - tvirial[5] += 0.5 * (delz * dely * cforce); + tvirial[0] += 0.5 * delx * delx * cforce; + tvirial[1] += 0.5 * dely * dely * cforce; + tvirial[2] += 0.5 * delz * delz * cforce; + tvirial[3] += 0.5 * dely * delx * cforce; + tvirial[4] += 0.5 * delz * delx * cforce; + tvirial[5] += 0.5 * delz * dely * cforce; } } else { @@ -413,6 +422,42 @@ void PairLJCutCoulLongTIP4P::init_style() if (!atom->q_flag) error->all("Pair style lj/cut/coul/long/tip4p requires atom attribute q"); + // request regular or rRESPA neighbor lists + + int irequest; + + if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) { + int respa = 0; + if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; + if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; + + if (respa == 0) irequest = neighbor->request(this); + else if (respa == 1) { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } else { + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 1; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respainner = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 2; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respamiddle = 1; + irequest = neighbor->request(this); + neighbor->requests[irequest]->id = 3; + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->respaouter = 1; + } + + } else irequest = neighbor->request(this); + cut_coulsq = cut_coul * cut_coul; // set & error check interior rRESPA cutoffs diff --git a/src/KSPACE/style_kspace.h b/src/KSPACE/style_kspace.h index 52cae87bdb..8c7244afb9 100644 --- a/src/KSPACE/style_kspace.h +++ b/src/KSPACE/style_kspace.h @@ -25,6 +25,7 @@ KSpaceStyle(pppm/tip4p,PPPMTIP4P) #ifdef PairInclude #include "pair_buck_coul_long.h" +#include "pair_coul_long.h" #include "pair_lj_cut_coul_long.h" #include "pair_lj_cut_coul_long_tip4p.h" #include "pair_lj_charmm_coul_long.h" @@ -32,6 +33,7 @@ KSpaceStyle(pppm/tip4p,PPPMTIP4P) #ifdef PairClass PairStyle(buck/coul/long,PairBuckCoulLong) +PairStyle(coul/long,PairCoulLong) PairStyle(lj/cut/coul/long,PairLJCutCoulLong) PairStyle(lj/cut/coul/long/tip4p,PairLJCutCoulLongTIP4P) PairStyle(lj/charmm/coul/long,PairLJCharmmCoulLong)