diff --git a/src/KSPACE/Install.sh b/src/KSPACE/Install.sh index 41c8849536..7f73650b53 100644 --- a/src/KSPACE/Install.sh +++ b/src/KSPACE/Install.sh @@ -3,17 +3,23 @@ if (test $1 = 1) then cp ewald.cpp .. + cp ewald_n.cpp .. cp pppm.cpp .. cp pppm_old.cpp .. cp pppm_cg.cpp .. + cp pppm_disp.cpp .. + cp pppm_disp_tip4p.cpp .. cp pppm_tip4p.cpp .. cp msm.cpp .. cp pair_born_coul_long.cpp .. cp pair_buck_coul_long.cpp .. + cp pair_buck_coul.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 .. + cp pair_lj_coul.cpp .. + cp pair_lj_coul_tip4p.cpp .. cp pair_born_coul_msm.cpp .. cp pair_buck_coul_msm.cpp .. cp pair_coul_msm.cpp .. @@ -25,18 +31,26 @@ if (test $1 = 1) then cp remap_wrap.cpp .. cp ewald.h .. + cp ewald_n.h .. cp kissfft.h .. + cp math_complex.h .. + cp math_vector.h .. cp msm.h .. cp pppm.h .. cp pppm_old.h .. cp pppm_cg.h .. + cp pppm_disp.h .. + cp pppm_disp_tip4p.h .. cp pppm_tip4p.h .. cp pair_born_coul_long.h .. cp pair_buck_coul_long.h .. + cp pair_buck_coul.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 .. + cp pair_lj_coul.h .. + cp pair_lj_coul_tip4p.h .. cp pair_born_coul_msm.h .. cp pair_buck_coul_msm.h .. cp pair_coul_msm.h .. @@ -50,17 +64,23 @@ if (test $1 = 1) then elif (test $1 = 0) then rm -f ../ewald.cpp + rm -f ../ewald_n.cpp rm -f ../msm.cpp rm -f ../pppm.cpp rm -f ../pppm_old.cpp rm -f ../pppm_cg.cpp + rm -f ../pppm_disp.cpp + rm -f ../pppm_disp_tip4p.cpp rm -f ../pppm_tip4p.cpp rm -f ../pair_born_coul_long.cpp rm -f ../pair_buck_coul_long.cpp + rm -f ../pair_buck_coul.cpp rm -f ../pair_coul_long.cpp rm -f ../pair_lj_cut_coul_long.cpp rm -f ../pair_lj_cut_coul_long_tip4p.cpp rm -f ../pair_lj_charmm_coul_long.cpp + rm -f ../pair_lj_coul.cpp + rm -f ../pair_lj_coul_tip4p.cpp rm -f ../pair_born_coul_msm.cpp rm -f ../pair_buck_coul_msm.cpp rm -f ../pair_coul_msm.cpp @@ -72,18 +92,26 @@ elif (test $1 = 0) then rm -f ../remap_wrap.cpp rm -f ../ewald.h + rm -f ../ewald_n.h rm -f ../kissfft.h + rm -f ../math_vector.h + rm -f ../math_complex.h rm -f ../msm.h rm -f ../pppm.h rm -f ../pppm_old.h rm -f ../pppm_cg.h + rm -f ../pppm_disp.h + rm -f ../pppm_disp_tip4p.h rm -f ../pppm_tip4p.h rm -f ../pair_born_coul_long.h rm -f ../pair_buck_coul_long.h + rm -f ../pair_buck_coul.h rm -f ../pair_coul_long.h rm -f ../pair_lj_cut_coul_long.h rm -f ../pair_lj_cut_coul_long_tip4p.h rm -f ../pair_lj_charmm_coul_long.h + rm -f ../pair_lj_coul.h + rm -f ../pair_lj_coul_tip4p.h rm -f ../pair_born_coul_msm.h rm -f ../pair_buck_coul_msm.h rm -f ../pair_coul_msm.h diff --git a/src/Makefile b/src/Makefile index 5d3116f7ad..036da146b7 100755 --- a/src/Makefile +++ b/src/Makefile @@ -18,8 +18,7 @@ PACKAGE = asphere class2 colloid dipole fld gpu granular kim \ rigid shock srd xtc PACKUSER = user-misc user-atc user-awpmd user-cg-cmm user-colvars \ - user-cuda user-eff user-ewaldn user-omp user-molfile \ - user-reaxc user-sph + user-cuda user-eff user-omp user-molfile user-reaxc user-sph PACKLIB = gpu kim meam poems reax user-atc user-awpmd user-colvars \ user-cuda user-molfile diff --git a/src/USER-EWALDN/Install.sh b/src/USER-EWALDN/Install.sh deleted file mode 100644 index 1cea246393..0000000000 --- a/src/USER-EWALDN/Install.sh +++ /dev/null @@ -1,33 +0,0 @@ -# Install/unInstall package files in LAMMPS - -if (test $1 = 1) then - - cp -p ewald_n.cpp .. - cp -p pair_buck_coul.cpp .. - #cp -p pair_lj_dipole.cpp .. - cp -p pair_lj_coul.cpp .. - - cp -p ewald_n.h .. - cp -p pair_buck_coul.h .. - #cp -p pair_lj_dipole.h .. - cp -p pair_lj_coul.h .. - - cp -p math_vector.h .. - cp -p math_complex.h .. - -elif (test $1 = 0) then - - rm -f ../ewald_n.cpp - rm -f ../pair_buck_coul.cpp - #rm -f ../pair_lj_dipole.cpp - rm -f ../pair_lj_coul.cpp - - rm -f ../ewald_n.h - rm -f ../pair_buck_coul.h - #rm -f ../pair_lj_dipole.h - rm -f ../pair_lj_coul.h - - rm -f ../math_vector.h - rm -f ../math_complex.h - -fi diff --git a/src/USER-EWALDN/README b/src/USER-EWALDN/README deleted file mode 100644 index 71e2da0f07..0000000000 --- a/src/USER-EWALDN/README +++ /dev/null @@ -1,29 +0,0 @@ -This package implements 3 commands which can be used in a LAMMPS input -script: pair_style lj/coul, pair_style buck/coul, and kspace_style -ewald/n. - -The "kspace_style ewald/n" command is similar to standard Ewald for -charges, but also enables the Lennard-Jones interaction, or any 1/r^N -interaction to be of infinite extent, instead of being cutoff. LAMMPS -pair potentials for long-range Coulombic interactions, such as -lj/cut/coul/long can be used with ewald/n. The two new pair_style -commands provide the modifications for the short-range LJ and -Buckingham interactions that can also be used with ewald/n. - -Two other advantages of kspace_style ewald/n are that - -a) it can be used with non-orthogonal (triclinic symmetry) simulation -boxes - -b) it can include long-range summations not just for Coulombic -interactions (1/r), but also for dispersion interactions (1/r^6) and -dipole interactions (1/r^3). - -Neither of these options is currently possible for other kspace styles -such as PPPM and ewald. - -See the doc pages for these commands for details. - -The person who created these files is Pieter in' t Veld while at -Sandia. He is now at BASF (pieter.intveld at basf.com). Contact him -directly if you have questions. diff --git a/src/USER-EWALDN/ewald_n.cpp b/src/USER-EWALDN/ewald_n.cpp deleted file mode 100644 index 8dd262a6e1..0000000000 --- a/src/USER-EWALDN/ewald_n.cpp +++ /dev/null @@ -1,1221 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Pieter in 't Veld (SNL) -------------------------------------------------------------------------- */ - -#include "mpi.h" -#include "string.h" -#include "stdio.h" -#include "stdlib.h" -#include "math.h" -#include "ewald_n.h" -#include "math_vector.h" -#include "math_const.h" -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "pair.h" -#include "domain.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define SMALL 0.00001 - -#define KSPACE_ILLEGAL "Illegal kspace_style ewald/n command" -#define KSPACE_ORDER "Unsupported order in kspace_style ewald/n for" -#define KSPACE_MIX "Unsupported mixing rule in kspace_style ewald/n for" - -enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // same as in pair.h - -//#define DEBUG - -/* ---------------------------------------------------------------------- */ - -EwaldN::EwaldN(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg) -{ - if (narg!=1) error->all(FLERR,KSPACE_ILLEGAL); - accuracy_relative = fabs(atof(arg[0])); - memset(function, 0, EWALD_NORDER*sizeof(int)); - kenergy = kvirial = NULL; - cek_local = cek_global = NULL; - ekr_local = NULL; - hvec = NULL; - kvec = NULL; - B = NULL; - first_output = 0; - energy_self_peratom = NULL; - virial_self_peratom = NULL; - nmax = 0; - q2 = 0; - b2 = 0; -} - -EwaldN::~EwaldN() -{ - deallocate(); - deallocate_peratom(); - delete [] ekr_local; - delete [] B; -} - - -/* --------------------------------------------------------------------- */ - -void EwaldN::init() -{ - nkvec = nkvec_max = nevec = nevec_max = 0; - nfunctions = nsums = sums = 0; - nbox = -1; - bytes = 0.0; - - if (!comm->me) { // output message - if (screen) fprintf(screen,"EwaldN initialization ...\n"); - if (logfile) fprintf(logfile,"EwaldN initialization ...\n"); - } - - if (domain->dimension == 2) // check for errors - error->all(FLERR,"Cannot use EwaldN with 2d simulation"); - if (slabflag == 0 && domain->nonperiodic > 0) - error->all(FLERR,"Cannot use nonperiodic boundaries with EwaldN"); - if (slabflag == 1) { - if (domain->xperiodic != 1 || domain->yperiodic != 1 || - domain->boundary[2][0] != 1 || domain->boundary[2][1] != 1) - error->all(FLERR,"Incorrect boundaries with slab EwaldN"); - } - - scale = 1.0; - //mumurd2e = force->mumurd2e; - //dielectric = force->dielectric; - mumurd2e = dielectric = 1.0; - - int tmp; - Pair *pair = force->pair; - int *ptr = pair ? (int *) pair->extract("ewald_order",tmp) : NULL; - double *cutoff = pair ? (double *) pair->extract("cut_coul",tmp) : NULL; - if (!(ptr||cutoff)) - error->all(FLERR,"KSpace style is incompatible with Pair style"); - int ewald_order = ptr ? *((int *) ptr) : 1<<1; - int ewald_mix = ptr ? *((int *) pair->extract("ewald_mix",tmp)) : GEOMETRIC; - memset(function, 0, EWALD_NFUNCS*sizeof(int)); - for (int i=0; i<=EWALD_NORDER; ++i) // transcribe order - if (ewald_order&(1<pair_style); - error->all(FLERR,str); - default: - sprintf(str, "%s pair_style %s", KSPACE_ORDER, force->pair_style); - error->all(FLERR,str); - } - nfunctions += function[k] = 1; - nsums += n[k]; - } - - - g_ewald = 0; - pair->init(); // so B is defined - init_coeffs(); - init_coeff_sums(); - - double qsum, qsqsum, bsbsum; - qsum = qsqsum = bsbsum = 0.0; - if (function[0]) { - qsum = sum[0].x; - qsqsum = sum[0].x2; - } - if (function[1]) { - bsbsum = sum[1].x2; - } - if (function[2]) { - bsbsum = sum[2].x2; - } - - - if (qsqsum == 0.0 && bsbsum == 0.0) - error->all(FLERR,"Cannot use Ewald/n solver on system with no charge or LJ particles"); - if (fabs(qsum) > SMALL && comm->me == 0) { - char str[128]; - sprintf(str,"System is not charge neutral, net charge = %g",qsum); - error->warning(FLERR,str); - } - - //set accuracy (force units) from accuracy_relative or accuracy_absolute - - if (accuracy_absolute >= 0.0) accuracy = accuracy_absolute; - else accuracy = accuracy_relative * two_charge_force; - - // setup K-space resolution - - q2 = qsqsum * force->qqrd2e / force->dielectric; - b2 = bsbsum; //Are these units right? - bigint natoms = atom->natoms; - - if (function[0]) { //Coulombic - g_ewald = accuracy*sqrt(natoms*(*cutoff)*shape_det(domain->h)) / (2.0*q2); - if (g_ewald >= 1.0) - error->all(FLERR,"KSpace accuracy too large to estimate G vector"); - g_ewald = sqrt(-log(g_ewald)) / *cutoff; - } - else if (function[1] || function[2]) { //Only LJ - - double *cutoffLJ = pair ? (double *) pair->extract("cut_LJ",tmp) : NULL; - - //Try Newton Solver - - //Use old method to get guess - g_ewald = (1.35 - 0.15*log(accuracy))/ *cutoffLJ; - - double g_ewald_new = NewtonSolve(g_ewald,(*cutoffLJ),natoms,shape_det(domain->h),b2); - - if (g_ewald_new > 0.0) g_ewald = g_ewald_new; - else error->warning(FLERR,"Ewald/n Newton solver failed, using old method to estimate g_ewald"); - - if (g_ewald >= 1.0) - error->all(FLERR,"KSpace accuracy too large to estimate G vector"); - } - - if (!comm->me) { // output results - if (screen) fprintf(screen, " G vector = %g\n", g_ewald); - if (logfile) fprintf(logfile, " G vector = %g\n", g_ewald); - } - - - deallocate_peratom(); - peratom_allocate_flag = 0; -} - - -/* ---------------------------------------------------------------------- - adjust EwaldN coeffs, called initially and whenever volume has changed -------------------------------------------------------------------------- */ - -void EwaldN::setup() -{ - volume = shape_det(domain->h)*slab_volfactor; // cell volume - memcpy(unit, domain->h_inv, sizeof(shape)); // wave vector units - shape_scalar_mult(unit, 2.0*MY_PI); - unit[2] /= slab_volfactor; - - //int nbox_old = nbox, nkvec_old = nkvec; - if (accuracy>=1) { - nbox = 0; - error->all(FLERR,"KSpace accuracy too low"); - } - - bigint natoms = atom->natoms; - double err; - int kxmax = 1; - int kymax = 1; - int kzmax = 1; - err = rms(kxmax,domain->h[0],natoms,q2,b2); - while (err > accuracy) { - kxmax++; - err = rms(kxmax,domain->h[0],natoms,q2,b2); - } - err = rms(kymax,domain->h[1],natoms,q2,b2); - while (err > accuracy) { - kymax++; - err = rms(kymax,domain->h[1],natoms,q2,b2); - } - err = rms(kzmax,domain->h[2]*slab_volfactor,natoms,q2,b2); - while (err > accuracy) { - kzmax++; - err = rms(kzmax,domain->h[2]*slab_volfactor,natoms,q2,b2); - } - nbox = MAX(kxmax,kymax); - nbox = MAX(nbox,kzmax); - double gsqxmx = unit[0]*unit[0]*kxmax*kxmax; - double gsqymx = unit[1]*unit[1]*kymax*kymax; - double gsqzmx = unit[2]*unit[2]*kzmax*kzmax; - gsqmx = MAX(gsqxmx,gsqymx); - gsqmx = MAX(gsqmx,gsqzmx); - gsqmx *= 1.00001; - - reallocate(); - coefficients(); // compute coeffs - init_coeffs(); - init_coeff_sums(); - init_self(); - - if (!(first_output||comm->me)) { // output on first - first_output = 1; - if (screen) fprintf(screen, - " vectors: nbox = %d, nkvec = %d\n", nbox, nkvec); - if (logfile) fprintf(logfile, - " vectors: nbox = %d, nkvec = %d\n", nbox, nkvec); - } -} - -/* ---------------------------------------------------------------------- - compute RMS accuracy for a dimension -------------------------------------------------------------------------- */ - -double EwaldN::rms(int km, double prd, bigint natoms, double q2, double b2) -{ - double value = 0.0; - - //Coulombic - - double g2 = g_ewald*g_ewald; - - value += 2.0*q2*g_ewald/prd * - sqrt(1.0/(MY_PI*km*natoms)) * - exp(-MY_PI*MY_PI*km*km/(g2*prd*prd)); - - //Lennard-Jones - - double g7 = g2*g2*g2*g_ewald; - - value += 4.0*b2*g7/3.0 * - sqrt(1.0/(MY_PI*natoms)) * - (exp(-MY_PI*MY_PI*km*km/(g2*prd*prd)) * - (MY_PI*km/(g_ewald*prd) + 1)); - - return value; -} - -void EwaldN::reallocate() // allocate memory -{ - int ix, iy, iz; - int nkvec_max = nkvec; - vector h; - - nkvec = 0; // determine size(kvec) - int *kflag = new int[(nbox+1)*(2*nbox+1)*(2*nbox+1)]; - int *flag = kflag; - - for (ix=0; ix<=nbox; ++ix) - for (iy=-nbox; iy<=nbox; ++iy) - for (iz=-nbox; iz<=nbox; ++iz) - if (!(ix||iy||iz)) *(flag++) = 0; - else if ((!ix)&&(iy<0)) *(flag++) = 0; - else if ((!(ix||iy))&&(iz<0)) *(flag++) = 0; // use symmetry - else { - h[0] = unit[0]*ix; - h[1] = unit[5]*ix+unit[1]*iy; - h[2] = unit[4]*ix+unit[3]*iy+unit[2]*iz; - if ((*(flag++) = h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=gsqmx)) ++nkvec; - } - - if (nkvec>nkvec_max) { - deallocate(); // free memory - hvec = new hvector[nkvec]; // hvec - bytes += (nkvec-nkvec_max)*sizeof(hvector); - kvec = new kvector[nkvec]; // kvec - bytes += (nkvec-nkvec_max)*sizeof(kvector); - kenergy = new double[nkvec*nfunctions]; // kenergy - bytes += (nkvec-nkvec_max)*nfunctions*sizeof(double); - kvirial = new double[6*nkvec*nfunctions]; // kvirial - bytes += 6*(nkvec-nkvec_max)*nfunctions*sizeof(double); - cek_local = new complex[nkvec*nsums]; // cek_local - bytes += (nkvec-nkvec_max)*nsums*sizeof(complex); - cek_global = new complex[nkvec*nsums]; // cek_global - bytes += (nkvec-nkvec_max)*nsums*sizeof(complex); - nkvec_max = nkvec; - } - - flag = kflag; // create index and - kvector *k = kvec; // wave vectors - hvector *hi = hvec; - for (ix=0; ix<=nbox; ++ix) - for (iy=-nbox; iy<=nbox; ++iy) - for (iz=-nbox; iz<=nbox; ++iz) - if (*(flag++)) { - hi->x = unit[0]*ix; - hi->y = unit[5]*ix+unit[1]*iy; - (hi++)->z = unit[4]*ix+unit[3]*iy+unit[2]*iz; - k->x = ix+nbox; k->y = iy+nbox; (k++)->z = iz+nbox; } - - delete [] kflag; -} - - -void EwaldN::reallocate_atoms() -{ - if (eflag_atom || vflag_atom) - if (atom->nlocal > nmax) { - deallocate_peratom(); - allocate_peratom(); - nmax = atom->nmax; - } - - if ((nevec = atom->nmax*(2*nbox+1))<=nevec_max) return; - delete [] ekr_local; - ekr_local = new cvector[nevec]; - bytes += (nevec-nevec_max)*sizeof(cvector); - nevec_max = nevec; -} - - -void EwaldN::allocate_peratom() -{ - memory->create(energy_self_peratom, - atom->nmax,EWALD_NFUNCS,"ewald/n:energy_self_peratom"); - memory->create(virial_self_peratom, - atom->nmax,EWALD_NFUNCS,"ewald/n:virial_self_peratom"); -} - - -void EwaldN::deallocate_peratom() // free memory -{ - memory->destroy(energy_self_peratom); - memory->destroy(virial_self_peratom); -} - - -void EwaldN::deallocate() // free memory -{ - delete [] hvec; hvec = NULL; - delete [] kvec; kvec = NULL; - delete [] kenergy; kenergy = NULL; - delete [] kvirial; kvirial = NULL; - delete [] cek_local; cek_local = NULL; - delete [] cek_global; cek_global = NULL; -} - - -void EwaldN::coefficients() // set up pre-factors -{ - vector h; - hvector *hi = hvec, *nh; - double eta2 = 0.25/(g_ewald*g_ewald); - double b1, b2, expb2, h1, h2, c1, c2; - double *ke = kenergy, *kv = kvirial; - int func0 = function[0], func12 = function[1]||function[2], - func3 = function[3]; - - for (nh = (hi = hvec)+nkvec; hintypes; - - if (function[1]) { // geometric 1/r^6 - double **b = (double **) force->pair->extract("B",tmp); - delete [] B; - B = new double[n+1]; - bytes += (n+1)*sizeof(double); - for (int i=0; i<=n; ++i) B[i] = sqrt(fabs(b[i][i])); - } - if (function[2]) { // arithmetic 1/r^6 - double **epsilon = (double **) force->pair->extract("epsilon",tmp); - double **sigma = (double **) force->pair->extract("sigma",tmp); - double eps_i, sigma_i, sigma_n, *bi = B = new double[7*n+7]; - double c[7] = { - 1.0, sqrt(6.0), sqrt(15.0), sqrt(20.0), sqrt(15.0), sqrt(6.0), 1.0}; - - if (!(epsilon&&sigma)) - error->all( - FLERR,"epsilon or sigma reference not set by pair style in ewald/n"); - for (int i=0; i<=n; ++i) { - eps_i = sqrt(epsilon[i][i]); - sigma_i = sigma[i][i]; - sigma_n = 1.0; - for (int j=0; j<7; ++j) { - *(bi++) = sigma_n*eps_i*c[j]; sigma_n *= sigma_i; - } - } - } -} - - -void EwaldN::init_coeff_sums() // sums based on atoms -{ - if (sums) return; // calculated only once - sums = 1; - - Sum sum_local[EWALD_MAX_NSUMS]; - - memset(sum_local, 0, EWALD_MAX_NSUMS*sizeof(Sum)); - if (function[0]) { // 1/r - double *q = atom->q, *qn = q+atom->nlocal; - for (double *i=q; itype, *ntype = type+atom->nlocal; - for (int *i=type; itype, *ntype = type+atom->nlocal; - for (int *i=type; imu) { // dipole - double *mu = atom->mu[0], *nmu = mu+4*atom->nlocal; - for (double *i = mu; i < nmu; i += 4) - sum_local[9].x2 += i[3]*i[3]; - } - MPI_Allreduce(sum_local, sum, 2*EWALD_MAX_NSUMS, MPI_DOUBLE, MPI_SUM, world); -} - - -void EwaldN::init_self() -{ - double g1 = g_ewald, g2 = g1*g1, g3 = g1*g2; - const double qscale = force->qqrd2e * scale; - - memset(energy_self, 0, EWALD_NFUNCS*sizeof(double)); // self energy - memset(virial_self, 0, EWALD_NFUNCS*sizeof(double)); - - if (function[0]) { // 1/r - virial_self[0] = -0.5*MY_PI*qscale/(g2*volume)*sum[0].x*sum[0].x; - energy_self[0] = sum[0].x2*qscale*g1/sqrt(MY_PI)-virial_self[0]; - } - if (function[1]) { // geometric 1/r^6 - virial_self[1] = MY_PI*sqrt(MY_PI)*g3/(6.0*volume)*sum[1].x*sum[1].x; - energy_self[1] = -sum[1].x2*g3*g3/12.0+virial_self[1]; - } - if (function[2]) { // arithmetic 1/r^6 - virial_self[2] = MY_PI*sqrt(MY_PI)*g3/(48.0*volume)*(sum[2].x*sum[8].x+ - sum[3].x*sum[7].x+sum[4].x*sum[6].x+0.5*sum[5].x*sum[5].x); - energy_self[2] = -sum[2].x2*g3*g3/3.0+virial_self[2]; - } - if (function[3]) { // dipole - virial_self[3] = 0; // in surface - energy_self[3] = sum[9].x2*mumurd2e*2.0*g3/3.0/sqrt(MY_PI)-virial_self[3]; - } -} - - -void EwaldN::init_self_peratom() -{ - if (!(vflag_atom || eflag_atom)) return; - - double g1 = g_ewald, g2 = g1*g1, g3 = g1*g2; - const double qscale = force->qqrd2e * scale; - double *energy = energy_self_peratom[0]; - double *virial = virial_self_peratom[0]; - int nlocal = atom->nlocal; - - memset(energy, 0, EWALD_NFUNCS*nlocal*sizeof(double)); - memset(virial, 0, EWALD_NFUNCS*nlocal*sizeof(double)); - - if (function[0]) { // 1/r - double *ei = energy; - double *vi = virial; - double ce = qscale*g1/sqrt(MY_PI); - double cv = -0.5*MY_PI*qscale/(g2*volume); - double *qi = atom->q, *qn = qi + nlocal; - for (; qi < qn; qi++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) { - double q = *qi; - *vi = cv*q*sum[0].x; - *ei = ce*q*q-vi[0]; - } - } - if (function[1]) { // geometric 1/r^6 - double *ei = energy+1; - double *vi = virial+1; - double ce = -g3*g3/12.0; - double cv = MY_PI*sqrt(MY_PI)*g3/(6.0*volume); - int *typei = atom->type, *typen = typei + atom->nlocal; - for (; typei < typen; typei++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) { - double b = B[*typei]; - *vi = cv*b*sum[1].x; - *ei = ce*b*b+vi[0]; - } - } - if (function[2]) { // arithmetic 1/r^6 - double *bi; - double *ei = energy+2; - double *vi = virial+2; - double ce = -g3*g3/3.0; - double cv = 0.5*MY_PI*sqrt(MY_PI)*g3/(48.0*volume); - int *typei = atom->type, *typen = typei + atom->nlocal; - for (; typei < typen; typei++, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) { - bi = B+7*typei[0]+7; - for (int k=2; k<9; ++k) *vi += cv*sum[k].x*(--bi)[0]; - - /* PJV 20120225: - should this be this instead? above implies an inverse dependence - seems to be the above way in original; i recall having tested - arithmetic mixing in the conception phase, but an extra test would - be prudent (pattern repeats in multiple functions below) - - bi = B+7*typei[0]; - for (int k=2; k<9; ++k) *vi += cv*sum[k].x*(bi++)[0]; - - */ - - *ei = ce*bi[0]*bi[6]+vi[0]; - } - } - if (function[3]&&atom->mu) { // dipole - double *ei = energy+3; - double *vi = virial+3; - double *imu = atom->mu[0], *nmu = imu+4*atom->nlocal; - double ce = mumurd2e*2.0*g3/3.0/sqrt(MY_PI); - for (; imu < nmu; imu += 4, vi += EWALD_NFUNCS, ei += EWALD_NFUNCS) { - *vi = 0; // in surface - *ei = ce*imu[3]*imu[3]-vi[0]; - } - } -} - - -/* ---------------------------------------------------------------------- - compute the EwaldN long-range force, energy, virial -------------------------------------------------------------------------- */ - -void EwaldN::compute(int eflag, int vflag) -{ - if (!nbox) return; - - // set energy/virial flags - // invoke allocate_peratom() if needed for first time - - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = eflag_global = vflag_global = eflag_atom = vflag_atom = 0; - - if (!peratom_allocate_flag && (eflag_atom || vflag_atom)) { - allocate_peratom(); - peratom_allocate_flag = 1; - nmax = atom->nmax; - } - - reallocate_atoms(); - init_self_peratom(); - compute_ek(); - compute_force(); - compute_surface(); - compute_energy(); - compute_energy_peratom(); - compute_virial(); - compute_virial_peratom(); -} - - -void EwaldN::compute_ek() -{ - cvector *ekr = ekr_local; - int lbytes = (2*nbox+1)*sizeof(cvector); - hvector *h = NULL; - kvector *k, *nk = kvec+nkvec; - cvector *z = new cvector[2*nbox+1]; - cvector z1, *zx, *zy, *zz, *zn = z+2*nbox; - complex *cek, zxyz, zxy = COMPLEX_NULL, cx = COMPLEX_NULL; - vector mui; - double *x = atom->x[0], *xn = x+3*atom->nlocal, *q = atom->q, qi = 0.0; - double bi = 0.0, ci[7]; - double *mu = atom->mu ? atom->mu[0] : NULL; - int i, kx, ky, n = nkvec*nsums, *type = atom->type, tri = domain->triclinic; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(cek_local, 0, n*sizeof(complex)); // reset sums - while (xx, 1, 0); C_SET(zz->y, 1, 0); C_SET(zz->z, 1, 0); // z[0] - if (tri) { // triclinic z[1] - C_ANGLE(z1.x, unit[0]*x[0]+unit[5]*x[1]+unit[4]*x[2]); - C_ANGLE(z1.y, unit[1]*x[1]+unit[3]*x[2]); - C_ANGLE(z1.z, x[2]*unit[2]); x += 3; - } - else { // orthogonal z[1] - C_ANGLE(z1.x, *(x++)*unit[0]); - C_ANGLE(z1.y, *(x++)*unit[1]); - C_ANGLE(z1.z, *(x++)*unit[2]); - } - for (; zzx, zz->x, z1.x); // 3D k-vector - C_RMULT(zy->y, zz->y, z1.y); C_CONJ(zx->y, zy->y); - C_RMULT(zy->z, zz->z, z1.z); C_CONJ(zx->z, zy->z); - } - kx = ky = -1; - cek = cek_local; - if (func[0]) qi = *(q++); - if (func[1]) bi = B[*type]; - if (func[2]) memcpy(ci, B+7*type[0], 7*sizeof(double)); - if (func[3]) { - memcpy(mui, mu, sizeof(vector)); - vec_scalar_mult(mui, mu[3]); - mu += 4; - h = hvec; - } - for (k=kvec; ky) { // based on order in - if (kx!=k->x) cx = z[kx = k->x].x; // reallocate - C_RMULT(zxy, z[ky = k->y].y, cx); - } - C_RMULT(zxyz, z[k->z].z, zxy); - if (func[0]) { - cek->re += zxyz.re*qi; (cek++)->im += zxyz.im*qi; - } - if (func[1]) { - cek->re += zxyz.re*bi; (cek++)->im += zxyz.im*bi; - } - if (func[2]) for (i=0; i<7; ++i) { - cek->re += zxyz.re*ci[i]; (cek++)->im += zxyz.im*ci[i]; - } - if (func[3]) { - register double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h; - cek->re += zxyz.re*muk; (cek++)->im += zxyz.im*muk; - } - } - ekr = (cvector *) ((char *) memcpy(ekr, z, lbytes)+lbytes); - ++type; - } - MPI_Allreduce(cek_local, cek_global, 2*n, MPI_DOUBLE, MPI_SUM, world); - - delete [] z; -} - - -void EwaldN::compute_force() -{ - kvector *k; - hvector *h, *nh; - cvector *z = ekr_local; - vector sum[EWALD_MAX_NSUMS], mui = COMPLEX_NULL; - complex *cek, zc, zx = COMPLEX_NULL, zxy = COMPLEX_NULL; - double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL; - double *mu = atom->mu ? atom->mu[0] : NULL; - const double qscale = force->qqrd2e * scale; - double *ke, c[EWALD_NFUNCS] = { - 8.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(12.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 8.0*MY_PI*mumurd2e/volume}; - double kt = 4.0*pow(g_ewald, 3.0)/3.0/sqrt(MY_PI)/c[3]; - int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; - int func[EWALD_NFUNCS]; - - if (atom->torque) t = atom->torque[0]; - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector)); // fj = -dE/dr = - for (; fy) { // based on order in - if (kx!=k->x) zx = z[kx = k->x].x; // reallocate - C_RMULT(zxy, z[ky = k->y].y, zx); - } - C_CRMULT(zc, z[k->z].z, zxy); - if (func[0]) { // 1/r - register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek; - sum[0][0] += h->x*im; sum[0][1] += h->y*im; sum[0][2] += h->z*im; - } - if (func[1]) { // geometric 1/r^6 - register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek; - sum[1][0] += h->x*im; sum[1][1] += h->y*im; sum[1][2] += h->z*im; - } - if (func[2]) { // arithmetic 1/r^6 - register double im, c = *(ke++); - for (i=2; i<9; ++i) { - im = c*(zc.im*cek->re+cek->im*zc.re); ++cek; - sum[i][0] += h->x*im; sum[i][1] += h->y*im; sum[i][2] += h->z*im; - } - } - if (func[3]) { // dipole - register double im = *(ke++)*(zc.im*cek->re+ - cek->im*zc.re)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek; - sum[9][0] += h->x*im; sum[9][1] += h->y*im; sum[9][2] += h->z*im; - } - } - if (func[0]) { // 1/r - register double qi = *(q++)*c[0]; - f[0] -= sum[0][0]*qi; f[1] -= sum[0][1]*qi; f[2] -= sum[0][2]*qi; - } - if (func[1]) { // geometric 1/r^6 - register double bi = B[*type]*c[1]; - f[0] -= sum[1][0]*bi; f[1] -= sum[1][1]*bi; f[2] -= sum[1][2]*bi; - } - if (func[2]) { // arithmetic 1/r^6 - register double *bi = B+7*type[0]+7; - for (i=2; i<9; ++i) { - register double c2 = (--bi)[0]*c[2]; - f[0] -= sum[i][0]*c2; f[1] -= sum[i][1]*c2; f[2] -= sum[i][2]*c2; - } - } - if (func[3]) { // dipole - f[0] -= sum[9][0]; f[1] -= sum[9][1]; f[2] -= sum[9][2]; - *(t++) -= mui[1]*sum[0][2]+mui[2]*sum[0][1]-mui[0]*kt; // torque - *(t++) -= mui[2]*sum[0][0]+mui[0]*sum[0][2]-mui[1]*kt; - *(t++) -= mui[0]*sum[0][1]+mui[1]*sum[0][0]-mui[2]*kt; - } - z = (cvector *) ((char *) z+lbytes); - ++type; - } -} - - -void EwaldN::compute_surface() -{ - if (!function[3]) return; - if (!atom->mu) return; - - vector sum_local = VECTOR_NULL, sum_total; - memset(sum_local, 0, sizeof(vector)); - double *i, *n, *mu = atom->mu[0]; - - for (n = (i = mu) + 4*atom->nlocal; i < n; ++i) { - register double di = i[3]; - sum_local[0] += di*(i++)[0]; - sum_local[1] += di*(i++)[0]; - sum_local[2] += di*(i++)[0]; - } - MPI_Allreduce(sum_local, sum_total, 3, MPI_DOUBLE, MPI_SUM, world); - - energy_self[3] += virial_self[3]; - virial_self[3] = - mumurd2e*(2.0*MY_PI*vec_dot(sum_total,sum_total)/(2.0*dielectric+1)/volume); - energy_self[3] -= virial_self[3]; - - if (!(vflag_atom || eflag_atom)) return; - - double *ei = energy_self_peratom[0]+3; - double *vi = virial_self_peratom[0]+3; - double cv = 2.0*mumurd2e*MY_PI/(2.0*dielectric+1)/volume; - - for (i = mu; i < n; i += 4, ei += EWALD_NFUNCS, vi += EWALD_NFUNCS) { - *ei += *vi; - *vi = cv*i[3]*(i[0]*sum_total[0]+i[1]*sum_total[1]+i[2]*sum_total[2]); - *ei -= *vi; - } -} - - -void EwaldN::compute_energy() -{ - energy = 0.0; - if (!eflag_global) return; - - complex *cek = cek_global; - double *ke = kenergy; - const double qscale = force->qqrd2e * scale; - double c[EWALD_NFUNCS] = { - 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; - double sum[EWALD_NFUNCS]; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(sum, 0, EWALD_NFUNCS*sizeof(double)); // reset sums - for (int k=0; kre*cek->re+cek->im*cek->im); ++cek; } - if (func[1]) { // geometric 1/r^6 - sum[1] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; } - if (func[2]) { // arithmetic 1/r^6 - register double r = - (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+ - (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+ - (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+ - 0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7; - sum[2] += *(ke++)*r; - } - if (func[3]) { // dipole - sum[3] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; } - } - for (int k=0; kq; - double *eatomj = eatom; - double *mu = atom->mu ? atom->mu[0] : NULL; - const double qscale = force->qqrd2e * scale; - double *ke = kenergy; - double c[EWALD_NFUNCS] = { - 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; - int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - for (int j = 0; j < atom->nlocal; j++, ++eatomj) { - k = kvec; - kx = ky = -1; - ke = kenergy; - cek = cek_global; - memset(sum, 0, EWALD_MAX_NSUMS*sizeof(double)); - if (func[3]) { - register double di = mu[3] * c[3]; - mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0]; - mu++; - } - for (nh = (h = hvec)+nkvec; hy) { // based on order in - if (kx!=k->x) zx = z[kx = k->x].x; // reallocate - C_RMULT(zxy, z[ky = k->y].y, zx); - } - C_CRMULT(zc, z[k->z].z, zxy); - if (func[0]) { // 1/r - sum[0] += *(ke++)*(cek->re*zc.re - cek->im*zc.im); ++cek; } - if (func[1]) { // geometric 1/r^6 - sum[1] += *(ke++)*(cek->re*zc.re - cek->im*zc.im); ++cek; } - if (func[2]) { // arithmetic 1/r^6 - register double im, c = *(ke++); - for (i=2; i<9; ++i) { - im = c*(cek->re*zc.re - cek->im*zc.im); ++cek; - sum[i] += im; - } - } - if (func[3]) { // dipole - sum[9] += *(ke++)*(cek->re*zc.re + - cek->im*zc.im)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek; - } - } - - if (func[0]) { // 1/r - register double qj = *(q++)*c[0]; - *eatomj += sum[0]*qj - energy_self_peratom[j][0]; - } - if (func[1]) { // geometric 1/r^6 - register double bj = B[*type]*c[1]; - *eatomj += sum[1]*bj - energy_self_peratom[j][1]; - } - if (func[2]) { // arithmetic 1/r^6 - register double *bj = B+7*type[0]+7; - for (i=2; i<9; ++i) { - register double c2 = (--bj)[0]*c[2]; - *eatomj += 0.5*sum[i]*c2; - } - *eatomj -= energy_self_peratom[j][2]; - } - if (func[3]) { // dipole - *eatomj += sum[9] - energy_self_peratom[j][3]; - } - z = (cvector *) ((char *) z+lbytes); - ++type; - } -} - - -#define swap(a, b) { register double t = a; a= b; b = t; } - -void EwaldN::compute_virial() -{ - memset(virial, 0, sizeof(shape)); - if (!vflag_global) return; - - complex *cek = cek_global; - double *kv = kvirial; - const double qscale = force->qqrd2e * scale; - double c[EWALD_NFUNCS] = { - 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; - shape sum[EWALD_NFUNCS]; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - memset(sum, 0, EWALD_NFUNCS*sizeof(shape)); - for (int k=0; kre*cek->re+cek->im*cek->im; ++cek; - sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r; - sum[0][3] += *(kv++)*r; sum[0][4] += *(kv++)*r; sum[0][5] += *(kv++)*r; - } - if (func[1]) { // geometric 1/r^6 - register double r = cek->re*cek->re+cek->im*cek->im; ++cek; - sum[1][0] += *(kv++)*r; sum[1][1] += *(kv++)*r; sum[1][2] += *(kv++)*r; - sum[1][3] += *(kv++)*r; sum[1][4] += *(kv++)*r; sum[1][5] += *(kv++)*r; - } - if (func[2]) { // arithmetic 1/r^6 - register double r = - (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+ - (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+ - (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+ - 0.5*(cek[3].re*cek[3].re+cek[3].im*cek[3].im); cek += 7; - sum[2][0] += *(kv++)*r; sum[2][1] += *(kv++)*r; sum[2][2] += *(kv++)*r; - sum[2][3] += *(kv++)*r; sum[2][4] += *(kv++)*r; sum[2][5] += *(kv++)*r; - } - if (func[3]) { - register double r = cek->re*cek->re+cek->im*cek->im; ++cek; - sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r; - sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r; - } - } - for (int k=0; kq; - double *vatomj = vatom[0]; - double *mu = atom->mu ? atom->mu[0] : NULL; - const double qscale = force->qqrd2e * scale; - double c[EWALD_NFUNCS] = { - 4.0*MY_PI*qscale/volume, 2.0*MY_PI*sqrt(MY_PI)/(24.0*volume), - 2.0*MY_PI*sqrt(MY_PI)/(192.0*volume), 4.0*MY_PI*mumurd2e/volume}; - shape sum[EWALD_MAX_NSUMS]; - int func[EWALD_NFUNCS]; - - memcpy(func, function, EWALD_NFUNCS*sizeof(int)); - int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; - for (int j = 0; j < atom->nlocal; j++, vatomj += 6) { - k = kvec; - kx = ky = -1; - kv = kvirial; - cek = cek_global; - memset(sum, 0, EWALD_MAX_NSUMS*sizeof(shape)); - if (func[3]) { - register double di = mu[3] * c[3]; - mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[1]; mui[2] = di*(mu++)[2]; - mu++; - } - for (nh = (h = hvec)+nkvec; hy) { // based on order in - if (kx!=k->x) zx = z[kx = k->x].x; // reallocate - C_RMULT(zxy, z[ky = k->y].y, zx); - } - C_CRMULT(zc, z[k->z].z, zxy); - if (func[0]) { // 1/r - register double r = cek->re*zc.re - cek->im*zc.im; ++cek; - sum[0][0] += *(kv++)*r; - sum[0][1] += *(kv++)*r; - sum[0][2] += *(kv++)*r; - sum[0][3] += *(kv++)*r; - sum[0][4] += *(kv++)*r; - sum[0][5] += *(kv++)*r; - } - if (func[1]) { // geometric 1/r^6 - register double r = cek->re*zc.re - cek->im*zc.im; ++cek; - sum[1][0] += *(kv++)*r; - sum[1][1] += *(kv++)*r; - sum[1][2] += *(kv++)*r; - sum[1][3] += *(kv++)*r; - sum[1][4] += *(kv++)*r; - sum[1][5] += *(kv++)*r; - } - if (func[2]) { // arithmetic 1/r^6 - register double r; - for (i=2; i<9; ++i) { - r = cek->re*zc.re - cek->im*zc.im; ++cek; - sum[i][0] += *(kv++)*r; - sum[i][1] += *(kv++)*r; - sum[i][2] += *(kv++)*r; - sum[i][3] += *(kv++)*r; - sum[i][4] += *(kv++)*r; - sum[i][5] += *(kv++)*r; - kv -= 6; - } - kv += 6; - } - if (func[3]) { // dipole - register double - r = (cek->re*zc.re - cek->im*zc.im) - *(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z); ++cek; - sum[9][0] += *(kv++)*r; - sum[9][1] += *(kv++)*r; - sum[9][2] += *(kv++)*r; - sum[9][3] += *(kv++)*r; - sum[9][4] += *(kv++)*r; - sum[9][5] += *(kv++)*r; - } - } - - if (func[0]) { // 1/r - register double qi = *(q++)*c[0]; - for (int n = 0; n < 6; n++) vatomj[n] += sum[0][n]*qi; - } - if (func[1]) { // geometric 1/r^6 - register double bi = B[*type]*c[1]; - for (int n = 0; n < 6; n++) vatomj[n] += sum[1][n]*bi; - } - if (func[2]) { // arithmetic 1/r^6 - register double *bj = B+7*type[0]+7; - for (i=2; i<9; ++i) { - register double c2 = (--bj)[0]*c[2]; - for (int n = 0; n < 6; n++) vatomj[n] += 0.5*sum[i][n]*c2; - } - } - if (func[3]) { // dipole - for (int n = 0; n < 6; n++) vatomj[n] += sum[9][n]; - } - - for (int k=0; kq; - double **x = atom->x; - int nlocal = atom->nlocal; - - double dipole = 0.0; - for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2]; - - // sum local contributions to get global dipole moment - - double dipole_all; - MPI_Allreduce(&dipole,&dipole_all,1,MPI_DOUBLE,MPI_SUM,world); - - // compute corrections - - const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; - const double qscale = force->qqrd2e * scale; - - if (eflag_global) energy += qscale * e_slabcorr; - - // per-atom energy - - if (eflag_atom) { - double efact = 2.0*MY_PI*dipole_all/volume; - for (int i = 0; i < nlocal; i++) eatom[i] += qscale * q[i]*x[i][2]*efact; - } - - // add on force corrections - - double ffact = -4.0*MY_PI*dipole_all/volume; - double **f = atom->f; - - for (int i = 0; i < nlocal; i++) f[i][2] += qscale * q[i]*ffact; -} - -/* ---------------------------------------------------------------------- - Newton solver used to find g_ewald for LJ systems - ------------------------------------------------------------------------- */ - -double EwaldN::NewtonSolve(double x, double Rc, bigint natoms, double vol, double b2) -{ - double dx,tol; - int maxit; - - maxit = 10000; //Maximum number of iterations - tol = 0.00001; //Convergence tolerance - - //Begin algorithm - - for (int i = 0; i < maxit; i++) { - dx = f(x,Rc,natoms,vol,b2) / derivf(x,Rc,natoms,vol,b2); - x = x - dx; //Update x - if (fabs(dx) < tol) return x; - } - return -1; -} - -/* ---------------------------------------------------------------------- - Calculate f(x) - ------------------------------------------------------------------------- */ - -double EwaldN::f(double x, double Rc, bigint natoms, double vol, double b2) -{ - double a = Rc*x; - double f = (4.0*MY_PI*b2*pow(x,4.0)/vol/sqrt((double)natoms)*erfc(a) * - (6.0*pow(a,-5.0) + 6.0*pow(a,-3.0) + 3.0/a + a) - accuracy); - return f; -} - -/* ---------------------------------------------------------------------- - Calculate numerical derivative f'(x) - ------------------------------------------------------------------------- */ - -double EwaldN::derivf(double x, double Rc, bigint natoms, double vol, double b2) -{ - double h = 0.000001; //Derivative step-size - return (f(x + h,Rc,natoms,vol,b2) - f(x,Rc,natoms,vol,b2)) / h; -} diff --git a/src/USER-EWALDN/ewald_n.h b/src/USER-EWALDN/ewald_n.h deleted file mode 100644 index 4bc4bc1d83..0000000000 --- a/src/USER-EWALDN/ewald_n.h +++ /dev/null @@ -1,95 +0,0 @@ -/* ---------------------------------------------------------------------- - 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 KSPACE_CLASS - -KSpaceStyle(ewald/n,EwaldN) - -#else - -#ifndef LMP_EWALD_N_H -#define LMP_EWALD_N_H - -#include "kspace.h" -#include "math_complex.h" - -namespace LAMMPS_NS { - -#define EWALD_NORDER 6 -#define EWALD_NFUNCS 4 -#define EWALD_MAX_NSUMS 10 -#define EWALD_NSUMS {1, 1, 7, 1} - -typedef struct cvector { complex x, y, z; } cvector; -typedef struct hvector { double x, y, z; } hvector; -typedef struct kvector { long x, y, z; } kvector; - -class EwaldN : public KSpace { - public: - EwaldN(class LAMMPS *, int, char **); - ~EwaldN(); - void init(); - void setup(); - void compute(int, int); - double memory_usage() {return bytes;} - - private: - double unit[6]; - int function[EWALD_NFUNCS], first_output; - - int nkvec, nkvec_max, nevec, nevec_max, - nbox, nfunctions, nsums, sums; - int peratom_allocate_flag; - int nmax; - double bytes; - double gsqmx,q2,b2; - double *kenergy, energy_self[EWALD_NFUNCS]; - double *kvirial, virial_self[EWALD_NFUNCS]; - double **energy_self_peratom; - double **virial_self_peratom; - cvector *ekr_local; - hvector *hvec; - kvector *kvec; - - double mumurd2e, dielectric, *B, volume; - struct Sum { double x, x2; } sum[EWALD_MAX_NSUMS]; - complex *cek_local, *cek_global; - - double rms(int, double, bigint, double, double); - void reallocate(); - void allocate_peratom(); - void reallocate_atoms(); - void deallocate(); - void deallocate_peratom(); - void coefficients(); - void init_coeffs(); - void init_coeff_sums(); - void init_self(); - void init_self_peratom(); - void compute_ek(); - void compute_force(); - void compute_surface(); - void compute_energy(); - void compute_energy_peratom(); - void compute_virial(); - void compute_virial_peratom(); - void compute_slabcorr(); - double NewtonSolve(double, double, bigint, double, double); - double f(double, double, bigint, double, double); - double derivf(double, double, bigint, double, double); -}; - -} - -#endif -#endif diff --git a/src/USER-EWALDN/math_complex.h b/src/USER-EWALDN/math_complex.h deleted file mode 100644 index 4f53b45288..0000000000 --- a/src/USER-EWALDN/math_complex.h +++ /dev/null @@ -1,72 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Pieter J. in 't Veld (SNL) -------------------------------------------------------------------------- */ - -#ifndef LMP_MATH_COMPLEX_H -#define LMP_MATH_COMPLEX_H - -#define COMPLEX_NULL {0, 0} - -namespace LAMMPS_NS { - -typedef struct complex { - double re, im; } complex; - -} - -#define C_MULT(d, x, y) { \ - d.re = x.re*y.re-x.im*y.im; \ - d.im = x.re*y.im+x.im*y.re; } - -#define C_RMULT(d, x, y) { \ - register complex t = x; \ - d.re = t.re*y.re-t.im*y.im; \ - d.im = t.re*y.im+t.im*y.re; } - -#define C_CRMULT(d, x, y) { \ - register complex t = x; \ - d.re = t.re*y.re-t.im*y.im; \ - d.im = -t.re*y.im-t.im*y.re; } - -#define C_SMULT(d, x, y) { \ - d.re = x.re*y; \ - d.im = x.im*y; } - -#define C_ADD(d, x, y) { \ - d.re = x.re+y.re; \ - d.im = x.im+y.im; } - -#define C_SUBTR(d, x, y) { \ - d.re = x.re-y.re; \ - d.im = x.im-y.im; } - -#define C_CONJ(d, x) { \ - d.re = x.re; \ - d.im = -x.im; } - -#define C_SET(d, x, y) { \ - d.re = x; \ - d.im = y; } - -#define C_ANGLE(d, angle) { \ - register double a = angle; \ - d.re = cos(a); \ - d.im = sin(a); } - -#define C_COPY(d, x) { \ - memcpy(&d, &x, sizeof(complex)); } - -#endif diff --git a/src/USER-EWALDN/math_vector.h b/src/USER-EWALDN/math_vector.h deleted file mode 100644 index da65e069bd..0000000000 --- a/src/USER-EWALDN/math_vector.h +++ /dev/null @@ -1,448 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Pieter J. in 't Veld (SNL) -------------------------------------------------------------------------- */ - -#ifndef LMP_MATH_VECTOR_H -#define LMP_MATH_VECTOR_H - -#include "math.h" -#include "string.h" - -#define VECTOR_NULL {0, 0, 0} -#define SHAPE_NULL {0, 0, 0, 0, 0, 0} -#define FORM_NULL {0, 0, 0, 0, 0, 0} -#define MATRIX_NULL {VECTOR_NULL, VECTOR_NULL, VECTOR_NULL} -#define VECTOR4_NULL {0, 0, 0, 0} -#define QUATERNION_NULL {0, 0, 0, 0} -#define FORM4_NULL {0, 0, 0, 0, 0, 0, 0, 0, 0, 0} - -#define FZERO 1e-15 -#define fzero(x) (((x)>-FZERO) && ((x)1) xtx[4] = xtx[2]; - if (!form_inv(xtx_inv, xtx)) return 0; - memcpy(eqn, xty, sizeof(vector)); - form_vec_dot(eqn, xtx_inv); - return 1; } - -// cubic regression: y = eqn[0] + eqn[1]*x + eqn[2]*x*x + eqn[3]*x*x*x - -inline int regress3(vector4 &eqn, int order, double *x, double *y, int n) { - form4 xtx = FORM4_NULL, xtx_inv; - vector4 xty = VECTOR4_NULL; - double xn, xi, yi; - int i; - - vec4_null(eqn); - xtx[0] = n; - if ((order = order%3)<0) order = -order; // max: cubic regress - if (order<1) xtx[1] = 1.0; - if (order<2) xtx[2] = 1.0; - if (order<3) xtx[3] = 1.0; - for (i=0; i1) xtx[8] = xtx[1]; - if (order>2) { xtx[6] = xtx[7]; xtx[5] = xtx[2]; } - if (!form4_inv(xtx_inv, xtx)) return 0; - memcpy(eqn, xty, sizeof(vector4)); - form4_vec4_dot(eqn, xtx_inv); - return 1; } - -} - -#endif diff --git a/src/USER-EWALDN/pair_buck_coul.cpp b/src/USER-EWALDN/pair_buck_coul.cpp deleted file mode 100644 index dd7d323852..0000000000 --- a/src/USER-EWALDN/pair_buck_coul.cpp +++ /dev/null @@ -1,1168 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Pieter J. in 't Veld (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "math_vector.h" -#include "pair_buck_coul.h" -#include "atom.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "force.h" -#include "kspace.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 - -/* ---------------------------------------------------------------------- */ - -PairBuckCoul::PairBuckCoul(LAMMPS *lmp) : Pair(lmp) -{ - respa_enable = 1; - ftable = NULL; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -#define PAIR_ILLEGAL "Illegal pair_style buck/coul command" -#define PAIR_CUTOFF "Only one cut-off allowed when requesting all long" -#define PAIR_MISSING "Cut-offs missing in pair_style buck/coul" -#define PAIR_LJ_OFF "LJ6 off not supported in pair_style buck/coul" -#define PAIR_COUL_CUT "Coulombic cut not supported in pair_style buck/coul" -#define PAIR_LARGEST "Using largest cut-off for buck/coul long long" -#define PAIR_MIX "Geometric mixing assumed for 1/r^6 coefficients" - -void PairBuckCoul::options(char **arg, int order) -{ - char *option[] = {"long", "cut", "off", NULL}; - int i; - - if (!*arg) error->all(FLERR,PAIR_ILLEGAL); - for (i=0; option[i]&&strcmp(arg[0], option[i]); ++i); - switch (i) { - default: error->all(FLERR,PAIR_ILLEGAL); - case 0: ewald_order |= 1<all(FLERR,"Illegal pair_style command"); - - ewald_order = 0; - ewald_off = 0; - options(arg, 6); - options(++arg, 1); - if (!comm->me && ewald_order&(1<<6)) error->warning(FLERR,PAIR_MIX); - if (!comm->me && ewald_order==((1<<1)|(1<<6))) error->warning(FLERR,PAIR_LARGEST); - if (!*(++arg)) error->all(FLERR,PAIR_MISSING); - if (ewald_off&(1<<6)) error->all(FLERR,PAIR_LJ_OFF); - if (!((ewald_order^ewald_off)&(1<<1))) error->all(FLERR,PAIR_COUL_CUT); - cut_buck_global = force->numeric(*(arg++)); - if (*arg&&(ewald_order&0x42==0x42)) error->all(FLERR,PAIR_CUTOFF); - if (narg == 4) cut_coul = force->numeric(*arg); - else cut_coul = cut_buck_global; - - if (allocated) { // reset explicit cuts - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i+1; j <= atom->ntypes; j++) - if (setflag[i][j]) cut_buck[i][j] = cut_buck_global; - } -} - -/* ---------------------------------------------------------------------- - free all arrays -------------------------------------------------------------------------- */ - -PairBuckCoul::~PairBuckCoul() -{ - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - memory->destroy(cut_buck_read); - memory->destroy(cut_buck); - memory->destroy(cut_bucksq); - memory->destroy(buck_a_read); - memory->destroy(buck_a); - memory->destroy(buck_c_read); - memory->destroy(buck_c); - memory->destroy(buck_rho_read); - memory->destroy(buck_rho); - memory->destroy(buck1); - memory->destroy(buck2); - memory->destroy(rhoinv); - memory->destroy(offset); - } - if (ftable) free_tables(); -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairBuckCoul::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(cut_buck_read,n+1,n+1,"pair:cut_buck_read"); - memory->create(cut_buck,n+1,n+1,"pair:cut_buck"); - memory->create(cut_bucksq,n+1,n+1,"pair:cut_bucksq"); - memory->create(buck_a_read,n+1,n+1,"pair:buck_a_read"); - memory->create(buck_a,n+1,n+1,"pair:buck_a"); - memory->create(buck_c_read,n+1,n+1,"pair:buck_c_read"); - memory->create(buck_c,n+1,n+1,"pair:buck_c"); - memory->create(buck_rho_read,n+1,n+1,"pair:buck_rho_read"); - memory->create(buck_rho,n+1,n+1,"pair:buck_rho"); - memory->create(buck1,n+1,n+1,"pair:buck1"); - memory->create(buck2,n+1,n+1,"pair:buck2"); - memory->create(rhoinv,n+1,n+1,"pair:rhoinv"); - memory->create(offset,n+1,n+1,"pair:offset"); -} - -/* ---------------------------------------------------------------------- - extract protected data from object -------------------------------------------------------------------------- */ - -void *PairBuckCoul::extract(const char *id, int &dim) -{ - char *ids[] = { - "B", "ewald_order", "ewald_cut", "ewald_mix", "cut_coul", NULL}; - void *ptrs[] = { - buck_c, &ewald_order, &cut_coul, &mix_flag, &cut_coul, NULL}; - int i; - - for (i=0; ids[i]&&strcmp(ids[i], id); ++i); - if (i == 0) dim = 2; - else dim = 0; - return ptrs[i]; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairBuckCoul::coeff(int narg, char **arg) -{ - if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients"); - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(*(arg++),atom->ntypes,ilo,ihi); - force->bounds(*(arg++),atom->ntypes,jlo,jhi); - - double buck_a_one = force->numeric(*(arg++)); - double buck_rho_one = force->numeric(*(arg++)); - double buck_c_one = force->numeric(*(arg++)); - - double cut_buck_one = cut_buck_global; - if (narg == 6) cut_buck_one = force->numeric(*(arg++)); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - buck_a_read[i][j] = buck_a_one; - buck_c_read[i][j] = buck_c_one; - buck_rho_read[i][j] = buck_rho_one; - cut_buck_read[i][j] = cut_buck_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairBuckCoul::init_style() -{ - - // require an atom style with charge defined - - if (!atom->q_flag && (ewald_order&(1<<1))) - error->all(FLERR, - "Invoking coulombic in pair style lj/coul requires atom attribute q"); - - // request regular or rRESPA neighbor lists - - int irequest; - - if (update->whichflag == 0 && strstr(update->integrate_style,"respa")) { - 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 rRESPA cutoffs - - if (strstr(update->integrate_style,"respa") && - ((Respa *) update->integrate)->level_inner >= 0) - cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = NULL; - - // ensure use of KSpace long-range solver, set g_ewald - - if (ewald_order&(1<<1)) { // r^-1 kspace - if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - g_ewald = force->kspace->g_ewald; - } - if (ewald_order&(1<<6)) { // r^-6 kspace - if (!force->kspace && strcmp(force->kspace_style,"ewald/n")) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - g_ewald = force->kspace->g_ewald; - } - - // setup force tables - - if (ncoultablebits) init_tables(); -} - -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairBuckCoul::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 PairBuckCoul::init_one(int i, int j) -{ - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); - - cut_buck[i][j] = cut_buck_read[i][j]; - buck_a[i][j] = buck_a_read[i][j]; - buck_c[i][j] = buck_c_read[i][j]; - buck_rho[i][j] = buck_rho_read[i][j]; - - double cut = MAX(cut_buck[i][j],cut_coul); - cutsq[i][j] = cut*cut; - cut_bucksq[i][j] = cut_buck[i][j] * cut_buck[i][j]; - - buck1[i][j] = buck_a[i][j]/buck_rho[i][j]; - buck2[i][j] = 6.0*buck_c[i][j]; - rhoinv[i][j] = 1.0/buck_rho[i][j]; - - // check interior rRESPA cutoff - - if (cut_respa && MIN(cut_buck[i][j],cut_coul) < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); - - if (offset_flag) { - double rexp = exp(-cut_buck[i][j]/buck_rho[i][j]); - offset[i][j] = buck_a[i][j]*rexp - buck_c[i][j]/pow(cut_buck[i][j],6.0); - } else offset[i][j] = 0.0; - - cutsq[j][i] = cutsq[i][j]; - cut_bucksq[j][i] = cut_bucksq[i][j]; - buck_a[j][i] = buck_a[i][j]; - buck_c[j][i] = buck_c[i][j]; - rhoinv[j][i] = rhoinv[i][j]; - buck1[j][i] = buck1[i][j]; - buck2[j][i] = buck2[i][j]; - offset[j][i] = offset[i][j]; - - return cut; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairBuckCoul::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(&buck_a_read[i][j],sizeof(double),1,fp); - fwrite(&buck_rho_read[i][j],sizeof(double),1,fp); - fwrite(&buck_c_read[i][j],sizeof(double),1,fp); - fwrite(&cut_buck_read[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairBuckCoul::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(&buck_a_read[i][j],sizeof(double),1,fp); - fread(&buck_rho_read[i][j],sizeof(double),1,fp); - fread(&buck_c_read[i][j],sizeof(double),1,fp); - fread(&cut_buck_read[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&buck_a_read[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&buck_rho_read[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&buck_c_read[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_buck_read[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairBuckCoul::write_restart_settings(FILE *fp) -{ - fwrite(&cut_buck_global,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&ewald_order,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairBuckCoul::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&cut_buck_global,sizeof(double),1,fp); - fread(&cut_coul,sizeof(double),1,fp); - fread(&offset_flag,sizeof(int),1,fp); - fread(&mix_flag,sizeof(int),1,fp); - fread(&ewald_order,sizeof(int),1,fp); - } - MPI_Bcast(&cut_buck_global,1,MPI_DOUBLE,0,world); - 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); - MPI_Bcast(&ewald_order,1,MPI_INT,0,world); -} - -/* ---------------------------------------------------------------------- - compute pair interactions -------------------------------------------------------------------------- */ - -void PairBuckCoul::compute(int eflag, int vflag) -{ - double evdwl,ecoul,fpair; - evdwl = ecoul = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - double **x = atom->x, *x0 = x[0]; - double **f = atom->f, *f0 = f[0], *fi = f0; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; - double qi = 0.0, qri = 0.0, *cutsqi, *cut_bucksqi, - *buck1i, *buck2i, *buckai, *buckci, *rhoinvi, *offseti; - double r, rsq, r2inv, force_coul, force_buck; - double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2; - vector xi, d; - - ineighn = (ineigh = list->ilist)+list->inum; - - for (; ineighfirstneigh[i])+list->numneigh[i]; - - for (; jneigh= cutsqi[typej = type[j]]) continue; - r2inv = 1.0/rsq; - r = sqrt(rsq); - - if (order1 && (rsq < cut_coulsq)) { // coulombic - if (!ncoultablebits || rsq <= tabinnersq) { // series real space - register double x = g_ewald*r; - register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x); - if (ni == 0) { - s *= g_ewald*exp(-x*x); - force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s; - if (eflag) ecoul = t; - } - else { // special case - register double f = s*(1.0-special_coul[ni])/r; - s *= g_ewald*exp(-x*x); - force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f; - if (eflag) ecoul = t-f; - } - } // table real space - else { - register union_int_float_t t; - t.f = rsq; - register const int k = (t.i & ncoulmask) >> ncoulshiftbits; - register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni == 0) { - force_coul = qiqj*(ftable[k]+f*dftable[k]); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); - } - else { // special case - t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); - force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f); - } - } - } - else force_coul = ecoul = 0.0; - - if (rsq < cut_bucksqi[typej]) { // buckingham - register double rn = r2inv*r2inv*r2inv, - expr = exp(-r*rhoinvi[typej]); - if (order6) { // long-range - register double x2 = g2*rsq, a2 = 1.0/x2; - x2 = a2*exp(-x2)*buckci[typej]; - if (ni == 0) { - force_buck = - r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; - if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2; - } - else { // special case - register double f = special_lj[ni], t = rn*(1.0-f); - force_buck = f*r*expr*buck1i[typej]- - g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]; - if (eflag) evdwl = f*expr*buckai[typej] - - g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej]; - } - } - else { // cut - if (ni == 0) { - force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]; - if (eflag) evdwl = expr*buckai[typej] - - rn*buckci[typej]-offseti[typej]; - } - else { // special case - register double f = special_lj[ni]; - force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]); - if (eflag) - evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]); - } - } - } - else force_buck = evdwl = 0.0; - - fpair = (force_coul+force_buck)*r2inv; - - if (newton_pair || j < nlocal) { - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*fpair; fj[0] -= f; - fi[1] += f = d[1]*fpair; fj[1] -= f; - fi[2] += f = d[2]*fpair; fj[2] -= f; - } - else { - fi[0] += d[0]*fpair; - fi[1] += d[1]*fpair; - fi[2] += d[2]*fpair; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,d[0],d[1],d[2]); - } - } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- */ - -void PairBuckCoul::compute_inner() -{ - double r, rsq, r2inv, force_coul = 0.0, force_buck, fpair; - - int *type = atom->type; - int nlocal = atom->nlocal; - double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - double cut_out_on = cut_respa[0]; - double cut_out_off = cut_respa[1]; - - double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; - int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); - double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi; - vector xi, d; - - ineighn = (ineigh = listinner->ilist)+listinner->inum; - - for (; ineighfirstneigh[i])+listinner->numneigh[i]; - - for (; jneigh= cut_out_off_sq) continue; - r2inv = 1.0/rsq; - r = sqrt(rsq); - - if (order1 && (rsq < cut_coulsq)) // coulombic - force_coul = ni == 0 ? - qri*q[j]/r : qri*q[j]/r*special_coul[ni]; - - if (rsq < cut_bucksqi[typej = type[j]]) { // buckingham - register double rn = r2inv*r2inv*r2inv, - expr = exp(-r*rhoinvi[typej]); - force_buck = ni == 0 ? - (r*expr*buck1i[typej]-rn*buck2i[typej]) : - (r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni]; - } - else force_buck = 0.0; - - fpair = (force_coul + force_buck) * r2inv; - - if (rsq > cut_out_on_sq) { // switching - register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); - } - - if (newton_pair || j < nlocal) { // force update - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*fpair; fj[0] -= f; - fi[1] += f = d[1]*fpair; fj[1] -= f; - fi[2] += f = d[2]*fpair; fj[2] -= f; - } - else { - fi[0] += d[0]*fpair; - fi[1] += d[1]*fpair; - fi[2] += d[2]*fpair; - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void PairBuckCoul::compute_middle() -{ - double r, rsq, r2inv, force_coul = 0.0, force_buck, fpair; - - int *type = atom->type; - int nlocal = atom->nlocal; - double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - double cut_in_off = cut_respa[0]; - double cut_in_on = cut_respa[1]; - double cut_out_on = cut_respa[2]; - double cut_out_off = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; - int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); - double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi; - vector xi, d; - - ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum; - - for (; ineighfirstneigh[i])+listmiddle->numneigh[i]; - - for (; jneigh= cut_out_off_sq) continue; - if (rsq <= cut_in_off_sq) continue; - r2inv = 1.0/rsq; - r = sqrt(rsq); - - if (order1 && (rsq < cut_coulsq)) // coulombic - force_coul = ni == 0 ? - qri*q[j]/r : qri*q[j]/r*special_coul[ni]; - - if (rsq < cut_bucksqi[typej = type[j]]) { // buckingham - register double rn = r2inv*r2inv*r2inv, - expr = exp(-r*rhoinvi[typej]); - force_buck = ni == 0 ? - (r*expr*buck1i[typej]-rn*buck2i[typej]) : - (r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni]; - } - else force_buck = 0.0; - - fpair = (force_coul + force_buck) * r2inv; - - if (rsq < cut_in_on_sq) { // switching - register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); - } - if (rsq > cut_out_on_sq) { - register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); - } - - if (newton_pair || j < nlocal) { // force update - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*fpair; fj[0] -= f; - fi[1] += f = d[1]*fpair; fj[1] -= f; - fi[2] += f = d[2]*fpair; fj[2] -= f; - } - else { - fi[0] += d[0]*fpair; - fi[1] += d[1]*fpair; - fi[2] += d[2]*fpair; - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void PairBuckCoul::compute_outer(int eflag, int vflag) -{ - double evdwl,ecoul,fpair; - evdwl = ecoul = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = 0; - - double **x = atom->x, *x0 = x[0]; - double **f = atom->f, *f0 = f[0], *fi = f0; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni, respa_flag; - double qi = 0.0, qri = 0.0, *cutsqi, *cut_bucksqi, - *buck1i, *buck2i, *buckai, *buckci, *rhoinvi, *offseti; - double r, rsq, r2inv, force_coul, force_buck; - double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2; - double respa_buck = 0.0, respa_coul = 0.0, frespa = 0.0; - vector xi, d; - - double cut_in_off = cut_respa[2]; - double cut_in_on = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - - ineighn = (ineigh = listouter->ilist)+listouter->inum; - - for (; ineighfirstneigh[i])+listouter->numneigh[i]; - - for (; jneigh= cutsqi[typej = type[j]]) continue; - r2inv = 1.0/rsq; - r = sqrt(rsq); - - if ((respa_flag = (rsq>cut_in_off_sq)&&(rsq> ncoulshiftbits; - register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni == 0) { - force_coul = qiqj*(ftable[k]+f*dftable[k]); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); - } - else { // correct for special - t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); - force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f); - } - } - } - else force_coul = respa_coul = ecoul = 0.0; - - if (rsq < cut_bucksqi[typej]) { // buckingham - register double rn = r2inv*r2inv*r2inv, - expr = exp(-r*rhoinvi[typej]); - if (respa_flag) respa_buck = ni == 0 ? // correct for respa - frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) : - frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni]; - if (order6) { // long-range form - register double x2 = g2*rsq, a2 = 1.0/x2; - x2 = a2*exp(-x2)*buckci[typej]; - if (ni == 0) { - force_buck = - r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; - if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2; - } - else { // correct for special - register double f = special_lj[ni], t = rn*(1.0-f); - force_buck = f*r*expr*buck1i[typej]- - g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]; - if (eflag) evdwl = f*expr*buckai[typej] - - g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej]; - } - } - else { // cut form - if (ni == 0) { - force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]; - if (eflag) - evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej]; - } - else { // correct for special - register double f = special_lj[ni]; - force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]); - if (eflag) - evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]); - } - } - } - else force_buck = respa_buck = evdwl = 0.0; - - fpair = (force_coul+force_buck)*r2inv; - frespa = fpair-(respa_coul+respa_buck)*r2inv; - - if (newton_pair || j < nlocal) { - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*frespa; fj[0] -= f; - fi[1] += f = d[1]*frespa; fj[1] -= f; - fi[2] += f = d[2]*frespa; fj[2] -= f; - } - else { - fi[0] += d[0]*frespa; - fi[1] += d[1]*frespa; - fi[2] += d[2]*frespa; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,d[0],d[1],d[2]); - } - } -} - -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairBuckCoul::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(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - rsq_lookup.i = i << ncoulshiftbits; - rsq_lookup.i |= masklo; - if (rsq_lookup.f < tabinnersq) { - rsq_lookup.i = i << ncoulshiftbits; - rsq_lookup.i |= maskhi; - } - r = sqrt(rsq_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - etable[i] = qqrd2e/r * derfc; - } else { - rtable[i] = rsq_lookup.f; - 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_lookup.f > cut_respa[2]*cut_respa[2]) { - if (rsq_lookup.f < 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_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - 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 = 0.0,v_tmp = 0.0; - itablemin = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - if (rsq_lookup.f < cut_coulsq) { - rsq_lookup.f = cut_coulsq; - r = sqrtf(rsq_lookup.f); - 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_lookup.f > cut_respa[2]*cut_respa[2]) { - if (rsq_lookup.f < 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_lookup.f - 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]; - } - } -} - -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairBuckCoul::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - -/* ---------------------------------------------------------------------- */ - -double PairBuckCoul::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_buck, - double &fforce) -{ - double f, r, r2inv, r6inv, force_coul, force_buck; - double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q; - - r = sqrt(rsq); - r2inv = 1.0/rsq; - double eng = 0.0; - - if ((ewald_order&2) && (rsq < cut_coulsq)) { // coulombic - if (!ncoultablebits || rsq <= tabinnersq) { // series real space - register double x = g_ewald*r; - register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x); - f = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x); - force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f; - eng += t-f; - } - else { // table real space - register union_int_float_t t; - t.f = rsq; - register const int k = (t.i & ncoulmask) >> ncoulshiftbits; - register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j]; - t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]); - force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); - eng += qiqj*(etable[k]+f*detable[k]-t.f); - } - } else force_coul = 0.0; - - if (rsq < cut_bucksq[itype][jtype]) { // buckingham - register double expr = factor_buck*exp(-sqrt(rsq)*rhoinv[itype][jtype]); - r6inv = r2inv*r2inv*r2inv; - if (ewald_order&64) { // long-range - register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_buck); - x2 = a2*exp(-x2)*buck_c[itype][jtype]; - force_buck = buck1[itype][jtype]*r*expr- - g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*buck2[itype][jtype]; - eng += buck_a[itype][jtype]*expr- - g6*((a2+1.0)*a2+0.5)*x2+t*buck_c[itype][jtype]; - } - else { // cut - force_buck = - buck1[itype][jtype]*r*expr-factor_buck*buck_c[itype][jtype]*r6inv; - eng += buck_a[itype][jtype]*expr- - factor_buck*(buck_c[itype][jtype]*r6inv-offset[itype][jtype]); - } - } else force_buck = 0.0; - - fforce = (force_coul+force_buck)*r2inv; - return eng; -} diff --git a/src/USER-EWALDN/pair_buck_coul.h b/src/USER-EWALDN/pair_buck_coul.h deleted file mode 100644 index 6c29b13419..0000000000 --- a/src/USER-EWALDN/pair_buck_coul.h +++ /dev/null @@ -1,76 +0,0 @@ -/* -*- 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(buck/coul,PairBuckCoul) - -#else - -#ifndef LMP_PAIR_BUCK_COUL_H -#define LMP_PAIR_BUCK_COUL_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairBuckCoul : public Pair { - public: - double cut_coul; - - PairBuckCoul(class LAMMPS *); - ~PairBuckCoul(); - virtual void compute(int, int); - - virtual void settings(int, char **); - void coeff(int, char **); - 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 *); - void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); - void *extract(const char *, int &); - - void compute_inner(); - void compute_middle(); - void compute_outer(int, int); - - protected: - double cut_buck_global; - double **cut_buck, **cut_buck_read, **cut_bucksq; - double cut_coulsq; - double **buck_a_read, **buck_a, **buck_c_read, **buck_c; - double **buck1, **buck2, **buck_rho_read, **buck_rho, **rhoinv, **offset; - double *cut_respa; - double g_ewald; - int ewald_order, ewald_off; - - double tabinnersq; - double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable; - double *etable, *detable, *ptable, *dptable, *vtable, *dvtable; - int ncoulshiftbits, ncoulmask; - - void options(char **arg, int order); - void allocate(); - void init_tables(); - void free_tables(); -}; - -} - -#endif -#endif diff --git a/src/USER-EWALDN/pair_lj_coul.cpp b/src/USER-EWALDN/pair_lj_coul.cpp deleted file mode 100644 index c38ad6b6f4..0000000000 --- a/src/USER-EWALDN/pair_lj_coul.cpp +++ /dev/null @@ -1,1160 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Pieter J. in 't Veld (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "math_vector.h" -#include "pair_lj_coul.h" -#include "atom.h" -#include "comm.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "force.h" -#include "kspace.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 - -/* ---------------------------------------------------------------------- */ - -PairLJCoul::PairLJCoul(LAMMPS *lmp) : Pair(lmp) -{ - respa_enable = 1; - ftable = NULL; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -#define PAIR_ILLEGAL "Illegal pair_style lj/coul command" -#define PAIR_CUTOFF "Only one cut-off allowed when requesting all long" -#define PAIR_MISSING "Cut-offs missing in pair_style lj/coul" -#define PAIR_COUL_CUT "Coulombic cut not supported in pair_style lj/coul" -#define PAIR_LARGEST "Using largest cut-off for lj/coul long long" -#define PAIR_MIX "Mixing forced for lj coefficients" - -void PairLJCoul::options(char **arg, int order) -{ - char *option[] = {"long", "cut", "off", NULL}; - int i; - - if (!*arg) error->all(FLERR,PAIR_ILLEGAL); - for (i=0; option[i]&&strcmp(arg[0], option[i]); ++i); - switch (i) { - default: error->all(FLERR,PAIR_ILLEGAL); - case 0: ewald_order |= 1<all(FLERR,"Illegal pair_style command"); - - ewald_off = 0; - ewald_order = 0; - options(arg, 6); - options(++arg, 1); - if (!comm->me && ewald_order&(1<<6)) error->warning(FLERR,PAIR_MIX); - if (!comm->me && ewald_order==((1<<1)|(1<<6))) error->warning(FLERR,PAIR_LARGEST); - if (!*(++arg)) error->all(FLERR,PAIR_MISSING); - if (!((ewald_order^ewald_off)&(1<<1))) error->all(FLERR,PAIR_COUL_CUT); - cut_lj_global = force->numeric(*(arg++)); - if (*arg&&(ewald_order&0x42==0x42)) error->all(FLERR,PAIR_CUTOFF); - if (narg == 4) cut_coul = force->numeric(*arg); - else cut_coul = cut_lj_global; - - if (allocated) { // reset explicit cuts - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i+1; j <= atom->ntypes; j++) - if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; - } -} - -/* ---------------------------------------------------------------------- - free all arrays -------------------------------------------------------------------------- */ - -PairLJCoul::~PairLJCoul() -{ - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - memory->destroy(cut_lj_read); - memory->destroy(cut_lj); - memory->destroy(cut_ljsq); - memory->destroy(epsilon_read); - memory->destroy(epsilon); - memory->destroy(sigma_read); - memory->destroy(sigma); - memory->destroy(lj1); - memory->destroy(lj2); - memory->destroy(lj3); - memory->destroy(lj4); - memory->destroy(offset); - } - if (ftable) free_tables(); -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairLJCoul::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(cut_lj_read,n+1,n+1,"pair:cut_lj_read"); - memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); - memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); - memory->create(epsilon_read,n+1,n+1,"pair:epsilon_read"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma_read,n+1,n+1,"pair:sigma_read"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(offset,n+1,n+1,"pair:offset"); -} - -/* ---------------------------------------------------------------------- - extract protected data from object -------------------------------------------------------------------------- */ - -void *PairLJCoul::extract(const char *id, int &dim) -{ - char *ids[] = { - "B", "sigma", "epsilon", "ewald_order", "ewald_cut", "ewald_mix", - "cut_coul", "cut_LJ", NULL}; - void *ptrs[] = { - lj4, sigma, epsilon, &ewald_order, &cut_coul, &mix_flag, &cut_coul, &cut_lj_global, NULL}; - int i; - - for (i=0; ids[i]&&strcmp(ids[i], id); ++i); - if (i <= 2) dim = 2; - else dim = 0; - return ptrs[i]; -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairLJCoul::coeff(int narg, char **arg) -{ - if (narg < 4 || narg > 5) error->all(FLERR,"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); - - double epsilon_one = force->numeric(arg[2]); - double sigma_one = force->numeric(arg[3]); - - double cut_lj_one = cut_lj_global; - if (narg == 5) cut_lj_one = force->numeric(arg[4]); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - epsilon_read[i][j] = epsilon_one; - sigma_read[i][j] = sigma_one; - cut_lj_read[i][j] = cut_lj_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairLJCoul::init_style() -{ - char *style1[] = {"ewald", "ewald/n", "pppm", NULL}; - char *style6[] = {"ewald/n", NULL}; - int i; - - // require an atom style with charge defined - - if (!atom->q_flag && (ewald_order&(1<<1))) - error->all(FLERR, - "Invoking coulombic in pair style lj/coul requires atom attribute q"); - - // request regular or rRESPA neighbor lists - - int irequest; - - if (update->whichflag == 0 && strstr(update->integrate_style,"respa")) { - 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 rRESPA cutoffs - - if (strstr(update->integrate_style,"respa") && - ((Respa *) update->integrate)->level_inner >= 0) - cut_respa = ((Respa *) update->integrate)->cutoff; - else cut_respa = NULL; - - // ensure use of KSpace long-range solver, set g_ewald - - if (ewald_order&(1<<1)) { // r^-1 kspace - if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - for (i=0; style1[i]&&strcmp(force->kspace_style, style1[i]); ++i); - if (!style1[i]) error->all(FLERR,"Pair style is incompatible with KSpace style"); - } - if (ewald_order&(1<<6)) { // r^-6 kspace - if (force->kspace == NULL) - error->all(FLERR,"Pair style is incompatible with KSpace style"); - for (i=0; style6[i]&&strcmp(force->kspace_style, style6[i]); ++i); - if (!style6[i]) error->all(FLERR,"Pair style is incompatible with KSpace style"); - } - if (force->kspace) g_ewald = force->kspace->g_ewald; - - // setup force tables - - if (ncoultablebits) init_tables(); -} - -/* ---------------------------------------------------------------------- - neighbor callback to inform pair style of neighbor list to use - regular or rRESPA -------------------------------------------------------------------------- */ - -void PairLJCoul::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 PairLJCoul::init_one(int i, int j) -{ - if ((ewald_order&(1<<6))||(setflag[i][j] == 0)) { - epsilon[i][j] = mix_energy(epsilon_read[i][i],epsilon_read[j][j], - sigma_read[i][i],sigma_read[j][j]); - sigma[i][j] = mix_distance(sigma_read[i][i],sigma_read[j][j]); - if (ewald_order&(1<<6)) - cut_lj[i][j] = cut_lj_global; - else - cut_lj[i][j] = mix_distance(cut_lj_read[i][i],cut_lj_read[j][j]); - } - else { - sigma[i][j] = sigma_read[i][j]; - epsilon[i][j] = epsilon_read[i][j]; - cut_lj[i][j] = cut_lj_read[i][j]; - } - - double cut = MAX(cut_lj[i][j], cut_coul); - cutsq[i][j] = cut*cut; - cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; - - 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); - - // check interior rRESPA cutoff - - if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all(FLERR,"Pair cutoff < Respa interior cutoff"); - - if (offset_flag) { - double ratio = sigma[i][j] / cut_lj[i][j]; - offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); - } else offset[i][j] = 0.0; - - cutsq[j][i] = cutsq[i][j]; - cut_ljsq[j][i] = cut_ljsq[i][j]; - lj1[j][i] = lj1[i][j]; - lj2[j][i] = lj2[i][j]; - lj3[j][i] = lj3[i][j]; - lj4[j][i] = lj4[i][j]; - offset[j][i] = offset[i][j]; - - return cut; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairLJCoul::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_read[i][j],sizeof(double),1,fp); - fwrite(&sigma_read[i][j],sizeof(double),1,fp); - fwrite(&cut_lj_read[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairLJCoul::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_read[i][j],sizeof(double),1,fp); - fread(&sigma_read[i][j],sizeof(double),1,fp); - fread(&cut_lj_read[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&epsilon_read[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma_read[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut_lj_read[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairLJCoul::write_restart_settings(FILE *fp) -{ - fwrite(&cut_lj_global,sizeof(double),1,fp); - fwrite(&cut_coul,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); - fwrite(&ewald_order,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairLJCoul::read_restart_settings(FILE *fp) -{ - if (comm->me == 0) { - fread(&cut_lj_global,sizeof(double),1,fp); - fread(&cut_coul,sizeof(double),1,fp); - fread(&offset_flag,sizeof(int),1,fp); - fread(&mix_flag,sizeof(int),1,fp); - fread(&ewald_order,sizeof(int),1,fp); - } - MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); - 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); - MPI_Bcast(&ewald_order,1,MPI_INT,0,world); -} - -/* ---------------------------------------------------------------------- - compute pair interactions -------------------------------------------------------------------------- */ - -void PairLJCoul::compute(int eflag, int vflag) -{ - double evdwl,ecoul,fpair; - evdwl = ecoul = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - double **x = atom->x, *x0 = x[0]; - double **f = atom->f, *f0 = f[0], *fi = f0; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; - double qi = 0.0, qri = 0.0; - double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti; - double rsq, r2inv, force_coul, force_lj; - double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2; - vector xi, d; - - ineighn = (ineigh = list->ilist)+list->inum; - - for (; ineighfirstneigh[i])+list->numneigh[i]; - - for (; jneigh= cutsqi[typej = type[j]]) continue; - r2inv = 1.0/rsq; - - if (order1 && (rsq < cut_coulsq)) { // coulombic - if (!ncoultablebits || rsq <= tabinnersq) { // series real space - register double r = sqrt(rsq), x = g_ewald*r; - register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x); - if (ni == 0) { - s *= g_ewald*exp(-x*x); - force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s; - if (eflag) ecoul = t; - } - else { // special case - r = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x); - force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r; - if (eflag) ecoul = t-r; - } - } // table real space - else { - register union_int_float_t t; - t.f = rsq; - register const int k = (t.i & ncoulmask)>>ncoulshiftbits; - register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni == 0) { - force_coul = qiqj*(ftable[k]+f*dftable[k]); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); - } - else { // special case - t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); - force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f); - } - } - } - else force_coul = ecoul = 0.0; - - if (rsq < cut_ljsqi[typej]) { // lj - if (order6) { // long-range lj - register double rn = r2inv*r2inv*r2inv; - register double x2 = g2*rsq, a2 = 1.0/x2; - x2 = a2*exp(-x2)*lj4i[typej]; - if (ni == 0) { - force_lj = - (rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; - if (eflag) - evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2; - } - else { // special case - register double f = special_lj[ni], t = rn*(1.0-f); - force_lj = f*(rn *= rn)*lj1i[typej]- - g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]; - if (eflag) - evdwl = f*rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[typej]; - } - } - else { // cut lj - register double rn = r2inv*r2inv*r2inv; - if (ni == 0) { - force_lj = rn*(rn*lj1i[typej]-lj2i[typej]); - if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]; - } - else { // special case - register double f = special_lj[ni]; - force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]); - if (eflag) - evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]); - } - } - } - else force_lj = evdwl = 0.0; - - fpair = (force_coul+force_lj)*r2inv; - - if (newton_pair || j < nlocal) { - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*fpair; fj[0] -= f; - fi[1] += f = d[1]*fpair; fj[1] -= f; - fi[2] += f = d[2]*fpair; fj[2] -= f; - } - else { - fi[0] += d[0]*fpair; - fi[1] += d[1]*fpair; - fi[2] += d[2]*fpair; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,d[0],d[1],d[2]); - } - } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- */ - -void PairLJCoul::compute_inner() -{ - double rsq, r2inv, force_coul = 0.0, force_lj, fpair; - - int *type = atom->type; - int nlocal = atom->nlocal; - double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - double cut_out_on = cut_respa[0]; - double cut_out_off = cut_respa[1]; - - double cut_out_diff = cut_out_off - cut_out_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; - int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); - double qri, *cut_ljsqi, *lj1i, *lj2i; - vector xi, d; - - ineighn = (ineigh = list->ilist)+list->inum; - - for (; ineighfirstneigh[i])+list->numneigh[i]; - - for (; jneigh= cut_out_off_sq) continue; - r2inv = 1.0/rsq; - - if (order1 && (rsq < cut_coulsq)) // coulombic - force_coul = ni == 0 ? - qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni]; - - if (rsq < cut_ljsqi[typej = type[j]]) { // lennard-jones - register double rn = r2inv*r2inv*r2inv; - force_lj = ni == 0 ? - rn*(rn*lj1i[typej]-lj2i[typej]) : - rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni]; - } - else force_lj = 0.0; - - fpair = (force_coul + force_lj) * r2inv; - - if (rsq > cut_out_on_sq) { // switching - register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); - } - - if (newton_pair || j < nlocal) { // force update - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*fpair; fj[0] -= f; - fi[1] += f = d[1]*fpair; fj[1] -= f; - fi[2] += f = d[2]*fpair; fj[2] -= f; - } - else { - fi[0] += d[0]*fpair; - fi[1] += d[1]*fpair; - fi[2] += d[2]*fpair; - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void PairLJCoul::compute_middle() -{ - double rsq, r2inv, force_coul = 0.0, force_lj, fpair; - - int *type = atom->type; - int nlocal = atom->nlocal; - double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - double cut_in_off = cut_respa[0]; - double cut_in_on = cut_respa[1]; - double cut_out_on = cut_respa[2]; - double cut_out_off = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_out_diff = cut_out_off - cut_out_on; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - double cut_out_on_sq = cut_out_on*cut_out_on; - double cut_out_off_sq = cut_out_off*cut_out_off; - - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; - int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); - double qri, *cut_ljsqi, *lj1i, *lj2i; - vector xi, d; - - ineighn = (ineigh = list->ilist)+list->inum; - - for (; ineighfirstneigh[i])+list->numneigh[i]; - - for (; jneigh= cut_out_off_sq) continue; - if (rsq <= cut_in_off_sq) continue; - r2inv = 1.0/rsq; - - if (order1 && (rsq < cut_coulsq)) // coulombic - force_coul = ni == 0 ? - qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni]; - - if (rsq < cut_ljsqi[typej = type[j]]) { // lennard-jones - register double rn = r2inv*r2inv*r2inv; - force_lj = ni == 0 ? - rn*(rn*lj1i[typej]-lj2i[typej]) : - rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni]; - } - else force_lj = 0.0; - - fpair = (force_coul + force_lj) * r2inv; - - if (rsq < cut_in_on_sq) { // switching - register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; - fpair *= rsw*rsw*(3.0 - 2.0*rsw); - } - if (rsq > cut_out_on_sq) { - register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; - fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); - } - - if (newton_pair || j < nlocal) { // force update - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*fpair; fj[0] -= f; - fi[1] += f = d[1]*fpair; fj[1] -= f; - fi[2] += f = d[2]*fpair; fj[2] -= f; - } - else { - fi[0] += d[0]*fpair; - fi[1] += d[1]*fpair; - fi[2] += d[2]*fpair; - } - } - } -} - -/* ---------------------------------------------------------------------- */ - -void PairLJCoul::compute_outer(int eflag, int vflag) -{ - double evdwl,ecoul,fpair; - evdwl = ecoul = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = 0; - - double **x = atom->x, *x0 = x[0]; - double **f = atom->f, *f0 = f[0], *fi = f0; - double *q = atom->q; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_coul = force->special_coul; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - double qqrd2e = force->qqrd2e; - - int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); - int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni, respa_flag; - double qi = 0.0, qri = 0.0; - double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti; - double rsq, r2inv, force_coul, force_lj; - double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2; - double respa_lj = 0.0, respa_coul = 0.0, frespa = 0.0; - vector xi, d; - - double cut_in_off = cut_respa[2]; - double cut_in_on = cut_respa[3]; - - double cut_in_diff = cut_in_on - cut_in_off; - double cut_in_off_sq = cut_in_off*cut_in_off; - double cut_in_on_sq = cut_in_on*cut_in_on; - - ineighn = (ineigh = list->ilist)+list->inum; - - for (; ineighfirstneigh[i])+list->numneigh[i]; - - for (; jneigh= cutsqi[typej = type[j]]) continue; - r2inv = 1.0/rsq; - - if ((respa_flag = (rsq>cut_in_off_sq)&&(rsq> ncoulshiftbits; - register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni == 0) { - force_coul = qiqj*(ftable[k]+f*dftable[k]); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); - } - else { // correct for special - t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); - force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); - if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f); - } - } - } - else force_coul = respa_coul = ecoul = 0.0; - - if (rsq < cut_ljsqi[typej]) { // lennard-jones - register double rn = r2inv*r2inv*r2inv; - if (respa_flag) respa_lj = ni == 0 ? // correct for respa - frespa*rn*(rn*lj1i[typej]-lj2i[typej]) : - frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni]; - if (order6) { // long-range form - register double x2 = g2*rsq, a2 = 1.0/x2; - x2 = a2*exp(-x2)*lj4i[typej]; - if (ni == 0) { - force_lj = - (rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; - if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2; - } - else { // correct for special - register double f = special_lj[ni], t = rn*(1.0-f); - force_lj = f*(rn *= rn)*lj1i[typej]- - g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]; - if (eflag) - evdwl = f*rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[typej]; - } - } - else { // cut form - if (ni == 0) { - force_lj = rn*(rn*lj1i[typej]-lj2i[typej]); - if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]; - } - else { // correct for special - register double f = special_lj[ni]; - force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]); - if (eflag) - evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]); - } - } - } - else force_lj = respa_lj = evdwl = 0.0; - - fpair = (force_coul+force_lj)*r2inv; - frespa = fpair-(respa_coul+respa_lj)*r2inv; - - if (newton_pair || j < nlocal) { - register double *fj = f0+(j+(j<<1)), f; - fi[0] += f = d[0]*frespa; fj[0] -= f; - fi[1] += f = d[1]*frespa; fj[1] -= f; - fi[2] += f = d[2]*frespa; fj[2] -= f; - } - else { - fi[0] += d[0]*frespa; - fi[1] += d[1]*frespa; - fi[2] += d[2]*frespa; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,d[0],d[1],d[2]); - } - } -} - -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairLJCoul::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(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - rsq_lookup.i = i << ncoulshiftbits; - rsq_lookup.i |= masklo; - if (rsq_lookup.f < tabinnersq) { - rsq_lookup.i = i << ncoulshiftbits; - rsq_lookup.i |= maskhi; - } - r = sqrtf(rsq_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - etable[i] = qqrd2e/r * derfc; - } else { - rtable[i] = rsq_lookup.f; - 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_lookup.f > cut_respa[2]*cut_respa[2]) { - if (rsq_lookup.f < 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_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - 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 = 0.0,v_tmp = 0.0; - itablemin = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - if (rsq_lookup.f < cut_coulsq) { - rsq_lookup.f = cut_coulsq; - r = sqrtf(rsq_lookup.f); - 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_lookup.f > cut_respa[2]*cut_respa[2]) { - if (rsq_lookup.f < 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_lookup.f - 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]; - } - } -} - -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairLJCoul::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - -/* ---------------------------------------------------------------------- */ - -double PairLJCoul::single(int i, int j, int itype, int jtype, - double rsq, double factor_coul, double factor_lj, - double &fforce) -{ - double r2inv, r6inv, force_coul, force_lj; - double g2 = g_ewald*g_ewald, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q; - - double eng = 0.0; - - r2inv = 1.0/rsq; - if ((ewald_order&2) && (rsq < cut_coulsq)) { // coulombic - if (!ncoultablebits || rsq <= tabinnersq) { // series real space - register double r = sqrt(rsq), x = g_ewald*r; - register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x); - r = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x); - force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r; - eng += t-r; - } - else { // table real space - register union_int_float_t t; - t.f = rsq; - register const int k = (t.i & ncoulmask) >> ncoulshiftbits; - register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j]; - t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]); - force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); - eng += qiqj*(etable[k]+f*detable[k]-t.f); - } - } else force_coul = 0.0; - - if (rsq < cut_ljsq[itype][jtype]) { // lennard-jones - r6inv = r2inv*r2inv*r2inv; - if (ewald_order&64) { // long-range - register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj); - x2 = a2*exp(-x2)*lj4[itype][jtype]; - force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]- - g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype]; - eng += factor_lj*r6inv*lj3[itype][jtype]- - g6*((a2+1.0)*a2+0.5)*x2+t*lj4[itype][jtype]; - } - else { // cut - force_lj = factor_lj*r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]); - eng += factor_lj*(r6inv*(r6inv*lj3[itype][jtype]- - lj4[itype][jtype])-offset[itype][jtype]); - } - } else force_lj = 0.0; - - fforce = (force_coul+force_lj)*r2inv; - return eng; -} diff --git a/src/USER-EWALDN/pair_lj_coul.h b/src/USER-EWALDN/pair_lj_coul.h deleted file mode 100644 index b497cd2546..0000000000 --- a/src/USER-EWALDN/pair_lj_coul.h +++ /dev/null @@ -1,75 +0,0 @@ -/* -*- 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(lj/coul,PairLJCoul) - -#else - -#ifndef LMP_PAIR_LJ_COUL_H -#define LMP_PAIR_LJ_COUL_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairLJCoul : public Pair { - public: - double cut_coul; - - PairLJCoul(class LAMMPS *); - virtual ~PairLJCoul(); - virtual void compute(int, int); - virtual void settings(int, char **); - void coeff(int, char **); - 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 *); - void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); - void *extract(const char *, int &); - - void compute_inner(); - void compute_middle(); - void compute_outer(int, int); - - protected: - double cut_lj_global; - double **cut_lj, **cut_lj_read, **cut_ljsq; - double cut_coulsq; - double **epsilon_read, **epsilon, **sigma_read, **sigma; - double **lj1, **lj2, **lj3, **lj4, **offset; - double *cut_respa; - double g_ewald; - int ewald_order, ewald_off; - - double tabinnersq; - double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable; - double *etable, *detable, *ptable, *dptable, *vtable, *dvtable; - int ncoulshiftbits, ncoulmask; - - void options(char **arg, int order); - void allocate(); - void init_tables(); - void free_tables(); -}; - -} - -#endif -#endif diff --git a/src/kspace.cpp b/src/kspace.cpp index 7e14ad329a..ce2249ede0 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -38,6 +38,11 @@ KSpace::KSpace(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp) order = 5; gridflag = 0; gewaldflag = 0; + + order_6 = 5; + gridflag_6 = 0; + gewaldflag_6 = 0; + slabflag = 0; differentiation_flag = 0; slab_volfactor = 1; @@ -272,12 +277,23 @@ void KSpace::modify_params(int narg, char **arg) if (nx_pppm == 0 && ny_pppm == 0 && nz_pppm == 0) gridflag = 0; else gridflag = 1; iarg += 4; + } else if (strcmp(arg[iarg],"mesh/disp") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal kspace_modify command"); + nx_pppm_6 = atoi(arg[iarg+1]); + ny_pppm_6 = atoi(arg[iarg+2]); + nz_pppm_6 = atoi(arg[iarg+3]); + if (nx_pppm_6 == 0 || ny_pppm_6 == 0 || nz_pppm_6 == 0) gridflag_6 = 0; + else gridflag_6 = 1; + iarg += 4; } else if (strcmp(arg[iarg],"order") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); order = atoi(arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"splitorder") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); + } else if (strcmp(arg[iarg],"order/disp") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); + order_6 = atoi(arg[iarg+1]); + } else if (strcmp(arg[iarg],"order/split") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); split_order = atoi(arg[iarg+1]); if (split_order < 2 || split_order > 6) error->all(FLERR,"Illegal kspace_modify command"); @@ -292,6 +308,12 @@ void KSpace::modify_params(int narg, char **arg) if (g_ewald == 0.0) gewaldflag = 0; else gewaldflag = 1; iarg += 2; + } else if (strcmp(arg[iarg],"gewald/disp") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); + g_ewald_6 = atof(arg[iarg+1]); + if (g_ewald_6 == 0.0) gewaldflag_6 = 0; + else gewaldflag_6 = 1; + iarg += 2; } else if (strcmp(arg[iarg],"slab") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal kspace_modify command"); if (strcmp(arg[iarg+1],"nozforce") == 0) { diff --git a/src/kspace.h b/src/kspace.h index 00c178d0b3..05ed26055d 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -23,14 +23,15 @@ class KSpace : protected Pointers { friend class FixOMP; public: double energy; // accumulated energy + double energy_1,energy_6; double virial[6]; // accumlated virial double *eatom,**vatom; // accumulated per-atom energy/virial double e2group; // accumulated group-group energy double f2group[3]; // accumulated group-group force - double g_ewald; + double g_ewald,g_ewald_6; int nx_pppm,ny_pppm,nz_pppm; - + int nx_pppm_6,ny_pppm_6,nz_pppm_6; int nx_msm_max,ny_msm_max,nz_msm_max; int group_group_enable; // 1 if style supports group/group calculation @@ -61,8 +62,10 @@ class KSpace : protected Pointers { double dgamma(const double &); protected: - int gridflag,gewaldflag,differentiation_flag; - int order,split_order; + int gridflag,gridflag_6; + int gewaldflag,gewaldflag_6; + int order,order_6,split_order; + int differentiation_flag; int slabflag; int suffix_flag; // suffix compatibility flag double scale; diff --git a/src/pair_coul_dsf.cpp b/src/pair_coul_dsf.cpp index ff3f73ad92..85700f81f8 100644 --- a/src/pair_coul_dsf.cpp +++ b/src/pair_coul_dsf.cpp @@ -327,7 +327,7 @@ double PairCoulDSF::single(int i, int j, int itype, int jtype, double rsq, /* ---------------------------------------------------------------------- */ -void *PairCoulDSF::extract(char *str, int &dim) +void *PairCoulDSF::extract(const char *str, int &dim) { if (strcmp(str,"cut_coul") == 0) { dim = 0; diff --git a/src/pair_coul_dsf.h b/src/pair_coul_dsf.h index 6a89e622c3..2403a43a26 100644 --- a/src/pair_coul_dsf.h +++ b/src/pair_coul_dsf.h @@ -38,7 +38,7 @@ class PairCoulDSF : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); - void *extract(char *, int &); + void *extract(const char *, int &); protected: double cut_coul,cut_coulsq; diff --git a/src/pair_lj_cut_coul_dsf.cpp b/src/pair_lj_cut_coul_dsf.cpp index 796eb7465a..10447edcb6 100644 --- a/src/pair_lj_cut_coul_dsf.cpp +++ b/src/pair_lj_cut_coul_dsf.cpp @@ -460,7 +460,7 @@ double PairLJCutCoulDSF::single(int i, int j, int itype, int jtype, double rsq, /* ---------------------------------------------------------------------- */ -void *PairLJCutCoulDSF::extract(char *str, int &dim) +void *PairLJCutCoulDSF::extract(const char *str, int &dim) { if (strcmp(str,"cut_coul") == 0) { dim = 0; diff --git a/src/pair_lj_cut_coul_dsf.h b/src/pair_lj_cut_coul_dsf.h index 73b1bd8b81..665ca88ddc 100644 --- a/src/pair_lj_cut_coul_dsf.h +++ b/src/pair_lj_cut_coul_dsf.h @@ -38,7 +38,7 @@ class PairLJCutCoulDSF : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); - void *extract(char *, int &); + void *extract(const char *, int &); protected: double cut_lj_global;