diff --git a/src/USER-DIELECTRIC/Install.sh b/src/USER-DIELECTRIC/Install.sh index fd540f3c7e..be6c4f6c9e 100644 --- a/src/USER-DIELECTRIC/Install.sh +++ b/src/USER-DIELECTRIC/Install.sh @@ -30,7 +30,7 @@ action () { # are installed, which in turn requires KSPACE if (test $1 = 1) then - if (test ! -e ../ppp.cpp) then + if (test ! -e ../pppm.cpp) then echo "Must install KSPACE package with USER-DIELECTRIC" exit 1 fi diff --git a/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp b/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp index 564a156a04..d43a552f15 100644 --- a/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp +++ b/src/USER-DIELECTRIC/pair_lj_long_coul_long_dielectric.cpp @@ -53,9 +53,9 @@ using namespace MathExtra; PairLJLongCoulLongDielectric::PairLJLongCoulLongDielectric(LAMMPS *lmp) : PairLJLongCoulLong(lmp) { respa_enable = 0; - cut_respa = NULL; - efield = NULL; - epot = NULL; + cut_respa = nullptr; + efield = nullptr; + epot = nullptr; nmax = 0; } diff --git a/src/USER-DIELECTRIC/pppm_disp_dielectric.cpp b/src/USER-DIELECTRIC/pppm_disp_dielectric.cpp new file mode 100644 index 0000000000..cd20fc3c7a --- /dev/null +++ b/src/USER-DIELECTRIC/pppm_disp_dielectric.cpp @@ -0,0 +1,872 @@ +/* ---------------------------------------------------------------------- + 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 authors: Trung Nguyen (Northwestern) +------------------------------------------------------------------------- */ + +#include "pppm_disp_dielectric.h" + +#include "angle.h" +#include "atom.h" +#include "atom_vec_dielectric.h" +#include "bond.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fft3d_wrap.h" +#include "force.h" +#include "gridcomm.h" +#include "math_const.h" +#include "memory.h" +#include "neighbor.h" +#include "pair.h" +#include "remap_wrap.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define MAXORDER 7 +#define OFFSET 16384 +#define SMALL 0.00001 +#define LARGE 10000.0 +#define EPS_HOC 1.0e-7 + +enum{REVERSE_RHO,REVERSE_RHO_GEOM,REVERSE_RHO_ARITH,REVERSE_RHO_NONE}; +enum{FORWARD_IK,FORWARD_AD,FORWARD_IK_PERATOM,FORWARD_AD_PERATOM, + FORWARD_IK_GEOM,FORWARD_AD_GEOM, + FORWARD_IK_PERATOM_GEOM,FORWARD_AD_PERATOM_GEOM, + FORWARD_IK_ARITH,FORWARD_AD_ARITH, + FORWARD_IK_PERATOM_ARITH,FORWARD_AD_PERATOM_ARITH, + FORWARD_IK_NONE,FORWARD_AD_NONE,FORWARD_IK_PERATOM_NONE, + FORWARD_AD_PERATOM_NONE}; + +#ifdef FFT_SINGLE +#define ZEROF 0.0f +#define ONEF 1.0f +#else +#define ZEROF 0.0 +#define ONEF 1.0 +#endif + +/* ---------------------------------------------------------------------- */ + +PPPMDispDielectric::PPPMDispDielectric(LAMMPS *lmp) : PPPMDisp(lmp) +{ + dipoleflag = 0; // turned off for now, until dipole works + group_group_enable = 0; + + mu_flag = 0; + + efield = nullptr; + phi = nullptr; + potflag = 0; + + avec = (AtomVecDielectric *) atom->style_match("dielectric"); + if (!avec) error->all(FLERR,"pppm/dielectric requires atom style dielectric"); +} + +/* ---------------------------------------------------------------------- */ + +PPPMDispDielectric::~PPPMDispDielectric() +{ + memory->destroy(efield); + memory->destroy(phi); +} + +/* ---------------------------------------------------------------------- + compute the PPPM long-range force, energy, virial +------------------------------------------------------------------------- */ + +void PPPMDispDielectric::compute(int eflag, int vflag) +{ + + int i; + + // set energy/virial flags + // invoke allocate_peratom() if needed for first time + + ev_init(eflag,vflag); + + if (evflag_atom && !peratom_allocate_flag) allocate_peratom(); + + // convert atoms from box to lamda coords + + if (triclinic == 0) boxlo = domain->boxlo; + else { + boxlo = domain->boxlo_lamda; + domain->x2lamda(atom->nlocal); + } + // extend size of per-atom arrays if necessary + + if (atom->nmax > nmax) { + + if (function[0]) { + memory->destroy(part2grid); + memory->destroy(efield); + memory->destroy(phi); + } + if (function[1] + function[2] + function[3]) memory->destroy(part2grid_6); + nmax = atom->nmax; + if (function[0]) { + memory->create(part2grid,nmax,3,"pppm/disp:part2grid"); + memory->create(efield,nmax,3,"pppm/disp:efield"); + memory->create(phi,nmax,"pppm/disp:phi"); + } + if (function[1] + function[2] + function[3]) + memory->create(part2grid_6,nmax,3,"pppm/disp:part2grid_6"); + } + + energy = 0.0; + energy_1 = 0.0; + energy_6 = 0.0; + if (vflag) for (i = 0; i < 6; i++) virial_6[i] = virial_1[i] = 0.0; + + // find grid points for all my particles + // distribute partcles' charges/dispersion coefficients on the grid + // communication between processors and remapping two fft + // Solution of poissons equation in k-space and backtransformation + // communication between processors + // calculation of forces + + if (function[0]) { + + // perform calculations for coulomb interactions only + + particle_map_c(delxinv,delyinv,delzinv,shift,part2grid,nupper,nlower, + nxlo_out,nylo_out,nzlo_out,nxhi_out,nyhi_out,nzhi_out); + + make_rho_c(); + + gc->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); + + brick2fft(nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in, + density_brick,density_fft,work1,remap); + + if (differentiation_flag == 1) { + poisson_ad(work1,work2,density_fft,fft1,fft2, + nx_pppm,ny_pppm,nz_pppm,nfft, + nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft, + nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in, + energy_1,greensfn, + virial_1,vg,vg2, + u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick); + + gc->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); + + fieldforce_c_ad(); + + if (vflag_atom) + gc->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); + + } else { + poisson_ik(work1,work2,density_fft,fft1,fft2, + nx_pppm,ny_pppm,nz_pppm,nfft, + nxlo_fft,nylo_fft,nzlo_fft,nxhi_fft,nyhi_fft,nzhi_fft, + nxlo_in,nylo_in,nzlo_in,nxhi_in,nyhi_in,nzhi_in, + energy_1,greensfn, + fkx,fky,fkz,fkx2,fky2,fkz2, + vdx_brick,vdy_brick,vdz_brick,virial_1,vg,vg2, + u_brick,v0_brick,v1_brick,v2_brick,v3_brick,v4_brick,v5_brick); + + gc->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); + + fieldforce_c_ik(); + + if (evflag_atom) + gc->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM, + gc_buf1,gc_buf2,MPI_FFT_SCALAR); + } + + if (evflag_atom) fieldforce_c_peratom(); + } + + if (function[1]) { + + // perform calculations for geometric mixing + + particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6, + nupper_6,nlower_6, + nxlo_out_6,nylo_out_6,nzlo_out_6, + nxhi_out_6,nyhi_out_6,nzhi_out_6); + + make_rho_g(); + + gc6->reverse_comm_kspace(this,1,sizeof(FFT_SCALAR),REVERSE_RHO_GEOM, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + brick2fft(nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6, + density_brick_g,density_fft_g,work1_6,remap_6); + + if (differentiation_flag == 1) { + poisson_ad(work1_6,work2_6,density_fft_g,fft1_6,fft2_6, + nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6, + nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6, + nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6, + energy_6,greensfn_6, + virial_6,vg_6,vg2_6, + u_brick_g,v0_brick_g,v1_brick_g,v2_brick_g, + v3_brick_g,v4_brick_g,v5_brick_g); + + gc6->forward_comm_kspace(this,1,sizeof(FFT_SCALAR),FORWARD_AD_GEOM, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + fieldforce_g_ad(); + + if (vflag_atom) + gc6->forward_comm_kspace(this,6,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_GEOM, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + } else { + poisson_ik(work1_6,work2_6,density_fft_g,fft1_6,fft2_6, + nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6, + nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6, + nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6, + energy_6,greensfn_6, + fkx_6,fky_6,fkz_6,fkx2_6,fky2_6,fkz2_6, + vdx_brick_g,vdy_brick_g,vdz_brick_g,virial_6,vg_6,vg2_6, + u_brick_g,v0_brick_g,v1_brick_g,v2_brick_g, + v3_brick_g,v4_brick_g,v5_brick_g); + + gc6->forward_comm_kspace(this,3,sizeof(FFT_SCALAR),FORWARD_IK_GEOM, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + fieldforce_g_ik(); + + if (evflag_atom) + gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_GEOM, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + } + + if (evflag_atom) fieldforce_g_peratom(); + } + + if (function[2]) { + + // perform calculations for arithmetic mixing + + particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6, + nupper_6,nlower_6, + nxlo_out_6,nylo_out_6,nzlo_out_6, + nxhi_out_6,nyhi_out_6,nzhi_out_6); + + make_rho_a(); + + gc6->reverse_comm_kspace(this,7,sizeof(FFT_SCALAR),REVERSE_RHO_ARITH, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + brick2fft_a(); + + if (differentiation_flag == 1) { + poisson_ad(work1_6,work2_6,density_fft_a3,fft1_6,fft2_6, + nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6, + nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6, + nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6, + energy_6,greensfn_6, + virial_6,vg_6,vg2_6, + u_brick_a3,v0_brick_a3,v1_brick_a3,v2_brick_a3, + v3_brick_a3,v4_brick_a3,v5_brick_a3); + poisson_2s_ad(density_fft_a0,density_fft_a6, + u_brick_a0,v0_brick_a0,v1_brick_a0,v2_brick_a0, + v3_brick_a0,v4_brick_a0,v5_brick_a0, + u_brick_a6,v0_brick_a6,v1_brick_a6,v2_brick_a6, + v3_brick_a6,v4_brick_a6,v5_brick_a6); + poisson_2s_ad(density_fft_a1,density_fft_a5, + u_brick_a1,v0_brick_a1,v1_brick_a1,v2_brick_a1, + v3_brick_a1,v4_brick_a1,v5_brick_a1, + u_brick_a5,v0_brick_a5,v1_brick_a5,v2_brick_a5, + v3_brick_a5,v4_brick_a5,v5_brick_a5); + poisson_2s_ad(density_fft_a2,density_fft_a4, + u_brick_a2,v0_brick_a2,v1_brick_a2,v2_brick_a2, + v3_brick_a2,v4_brick_a2,v5_brick_a2, + u_brick_a4,v0_brick_a4,v1_brick_a4,v2_brick_a4, + v3_brick_a4,v4_brick_a4,v5_brick_a4); + + gc6->forward_comm_kspace(this,7,sizeof(FFT_SCALAR),FORWARD_AD_ARITH, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + fieldforce_a_ad(); + + if (evflag_atom) + gc6->forward_comm_kspace(this,42,sizeof(FFT_SCALAR),FORWARD_AD_PERATOM_ARITH, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + } else { + poisson_ik(work1_6,work2_6,density_fft_a3,fft1_6,fft2_6, + nx_pppm_6,ny_pppm_6,nz_pppm_6,nfft_6, + nxlo_fft_6,nylo_fft_6,nzlo_fft_6,nxhi_fft_6,nyhi_fft_6,nzhi_fft_6, + nxlo_in_6,nylo_in_6,nzlo_in_6,nxhi_in_6,nyhi_in_6,nzhi_in_6, + energy_6,greensfn_6, + fkx_6,fky_6,fkz_6,fkx2_6,fky2_6,fkz2_6, + vdx_brick_a3,vdy_brick_a3,vdz_brick_a3,virial_6,vg_6,vg2_6, + u_brick_a3,v0_brick_a3,v1_brick_a3,v2_brick_a3, + v3_brick_a3,v4_brick_a3,v5_brick_a3); + poisson_2s_ik(density_fft_a0,density_fft_a6, + vdx_brick_a0,vdy_brick_a0,vdz_brick_a0, + vdx_brick_a6,vdy_brick_a6,vdz_brick_a6, + u_brick_a0,v0_brick_a0,v1_brick_a0,v2_brick_a0, + v3_brick_a0,v4_brick_a0,v5_brick_a0, + u_brick_a6,v0_brick_a6,v1_brick_a6,v2_brick_a6, + v3_brick_a6,v4_brick_a6,v5_brick_a6); + poisson_2s_ik(density_fft_a1,density_fft_a5, + vdx_brick_a1,vdy_brick_a1,vdz_brick_a1, + vdx_brick_a5,vdy_brick_a5,vdz_brick_a5, + u_brick_a1,v0_brick_a1,v1_brick_a1,v2_brick_a1, + v3_brick_a1,v4_brick_a1,v5_brick_a1, + u_brick_a5,v0_brick_a5,v1_brick_a5,v2_brick_a5, + v3_brick_a5,v4_brick_a5,v5_brick_a5); + poisson_2s_ik(density_fft_a2,density_fft_a4, + vdx_brick_a2,vdy_brick_a2,vdz_brick_a2, + vdx_brick_a4,vdy_brick_a4,vdz_brick_a4, + u_brick_a2,v0_brick_a2,v1_brick_a2,v2_brick_a2, + v3_brick_a2,v4_brick_a2,v5_brick_a2, + u_brick_a4,v0_brick_a4,v1_brick_a4,v2_brick_a4, + v3_brick_a4,v4_brick_a4,v5_brick_a4); + + gc6->forward_comm_kspace(this,21,sizeof(FFT_SCALAR),FORWARD_IK_ARITH, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + fieldforce_a_ik(); + + if (evflag_atom) + gc6->forward_comm_kspace(this,49,sizeof(FFT_SCALAR),FORWARD_IK_PERATOM_ARITH, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + } + + if (evflag_atom) fieldforce_a_peratom(); + } + + if (function[3]) { + + // perform calculations if no mixing rule applies + + particle_map(delxinv_6,delyinv_6,delzinv_6,shift_6,part2grid_6, + nupper_6,nlower_6, + nxlo_out_6,nylo_out_6,nzlo_out_6, + nxhi_out_6,nyhi_out_6,nzhi_out_6); + + make_rho_none(); + + gc6->reverse_comm_kspace(this,nsplit_alloc,sizeof(FFT_SCALAR),REVERSE_RHO_NONE, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + brick2fft_none(); + + if (differentiation_flag == 1) { + int n = 0; + for (int k = 0; k < nsplit_alloc/2; k++) { + poisson_none_ad(n,n+1,density_fft_none[n],density_fft_none[n+1], + u_brick_none[n],u_brick_none[n+1], + v0_brick_none,v1_brick_none,v2_brick_none, + v3_brick_none,v4_brick_none,v5_brick_none); + n += 2; + } + + gc6->forward_comm_kspace(this,1*nsplit_alloc,sizeof(FFT_SCALAR), + FORWARD_AD_NONE, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + fieldforce_none_ad(); + + if (vflag_atom) + gc6->forward_comm_kspace(this,6*nsplit_alloc,sizeof(FFT_SCALAR), + FORWARD_AD_PERATOM_NONE, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + } else { + int n = 0; + for (int k = 0; k < nsplit_alloc/2; k++) { + poisson_none_ik(n,n+1,density_fft_none[n],density_fft_none[n+1], + vdx_brick_none[n],vdy_brick_none[n],vdz_brick_none[n], + vdx_brick_none[n+1],vdy_brick_none[n+1],vdz_brick_none[n+1], + u_brick_none,v0_brick_none,v1_brick_none,v2_brick_none, + v3_brick_none,v4_brick_none,v5_brick_none); + n += 2; + } + + gc6->forward_comm_kspace(this,3*nsplit_alloc,sizeof(FFT_SCALAR), + FORWARD_IK_NONE, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + + fieldforce_none_ik(); + + if (evflag_atom) + gc6->forward_comm_kspace(this,7*nsplit_alloc,sizeof(FFT_SCALAR), + FORWARD_IK_PERATOM_NONE, + gc6_buf1,gc6_buf2,MPI_FFT_SCALAR); + } + + if (evflag_atom) fieldforce_none_peratom(); + } + + // update qsum and qsqsum, if atom count has changed and energy needed + + if ((eflag_global || eflag_atom) && atom->natoms != natoms_original) { + qsum_qsq(); + natoms_original = atom->natoms; + } + + // sum energy across procs and add in volume-dependent term + + const double qscale = force->qqrd2e * scale; + + if (eflag_global) { + double energy_all; + MPI_Allreduce(&energy_1,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); + energy_1 = energy_all; + MPI_Allreduce(&energy_6,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); + energy_6 = energy_all; + + energy_1 *= 0.5*volume; + energy_6 *= 0.5*volume; + + energy_1 -= g_ewald*qsqsum/MY_PIS + + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); + energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij + + 1.0/12.0*pow(g_ewald_6,6)*csum; + energy_1 *= qscale; + } + + // sum virial across procs + + if (vflag_global) { + double virial_all[6]; + MPI_Allreduce(virial_1,virial_all,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i]; + MPI_Allreduce(virial_6,virial_all,6,MPI_DOUBLE,MPI_SUM,world); + for (i = 0; i < 6; i++) virial[i] += 0.5*volume*virial_all[i]; + if (function[1]+function[2]+function[3]) { + double a = MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij; + virial[0] -= a; + virial[1] -= a; + virial[2] -= a; + } + } + + if (eflag_atom) { + if (function[0]) { + double *q = atom->q; + // coulomb self energy correction + for (i = 0; i < atom->nlocal; i++) { + eatom[i] -= qscale*g_ewald*q[i]*q[i]/MY_PIS + + qscale*MY_PI2*q[i]*qsum / (g_ewald*g_ewald*volume); + } + } + if (function[1] + function[2] + function[3]) { + int tmp; + for (i = 0; i < atom->nlocal; i++) { + tmp = atom->type[i]; + eatom[i] += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp] + + 1.0/12.0*pow(g_ewald_6,6)*cii[tmp]; + } + } + } + + if (vflag_atom) { + if (function[1] + function[2] + function[3]) { + int tmp; + // dispersion self virial correction + for (i = 0; i < atom->nlocal; i++) { + tmp = atom->type[i]; + for (int n = 0; n < 3; n++) + vatom[i][n] -= MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumi[tmp]; + } + } + } + + // 2d slab correction + + if (slabflag) slabcorr(eflag); + if (function[0]) energy += energy_1; + if (function[1] + function[2] + function[3]) energy += energy_6; + + // convert atoms back from lamda to box coords + + if (triclinic) domain->lamda2x(atom->nlocal); +} + +/* ---------------------------------------------------------------------- + interpolate from grid to get electric field & force on my particles + for ik scheme +------------------------------------------------------------------------- */ + +void PPPMDispDielectric::fieldforce_c_ik() +{ + int i,l,m,n,nx,ny,nz,mx,my,mz; + FFT_SCALAR dx,dy,dz,x0,y0,z0; + FFT_SCALAR ekx,eky,ekz,u; + + // loop over my charges, interpolate electric field from nearby grid points + // (nx,ny,nz) = global coords of grid pt to "lower left" of charge + // (dx,dy,dz) = distance to "lower left" grid pt + // (mx,my,mz) = global coords of moving stencil pt + // ek = 3 components of E-field on particle + + double *q = atom->q; + double **x = atom->x; + double **f = atom->f; + double *eps = atom->epsilon; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) { + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; + + compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d); + + u = ekx = eky = ekz = ZEROF; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + z0 = rho1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + y0 = z0*rho1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + x0 = y0*rho1d[0][l]; + if (potflag) u += x0*u_brick[mz][my][mx]; + ekx -= x0*vdx_brick[mz][my][mx]; + eky -= x0*vdy_brick[mz][my][mx]; + ekz -= x0*vdz_brick[mz][my][mx]; + } + } + } + + // electrostatic potential + + if (potflag) phi[i] = u; + + // convert E-field to force + + const double efactor = scale * eps[i]; + efield[i][0] = efactor*ekx; + efield[i][1] = efactor*eky; + efield[i][2] = efactor*ekz; + + // convert E-field to force + + const double qfactor = force->qqrd2e * scale * q[i]; + f[i][0] += qfactor*ekx; + f[i][1] += qfactor*eky; + if (slabflag != 2) f[i][2] += qfactor*ekz; + } +} +/* ---------------------------------------------------------------------- + interpolate from grid to get electric field & force on my particles + for ad scheme +------------------------------------------------------------------------- */ + +void PPPMDispDielectric::fieldforce_c_ad() +{ + int i,l,m,n,nx,ny,nz,mx,my,mz; + FFT_SCALAR dx,dy,dz; + FFT_SCALAR ekx,eky,ekz,u; + double s1,s2,s3; + double sf = 0.0; + + double *prd; + + if (triclinic == 0) prd = domain->prd; + else prd = domain->prd_lamda; + + double xprd = prd[0]; + double yprd = prd[1]; + double zprd = prd[2]; + double zprd_slab = zprd*slab_volfactor; + + double hx_inv = nx_pppm/xprd; + double hy_inv = ny_pppm/yprd; + double hz_inv = nz_pppm/zprd_slab; + + // loop over my charges, interpolate electric field from nearby grid points + // (nx,ny,nz) = global coords of grid pt to "lower left" of charge + // (dx,dy,dz) = distance to "lower left" grid pt + // (mx,my,mz) = global coords of moving stencil pt + // ek = 3 components of E-field on particle + + double *q = atom->q; + double **x = atom->x; + double **f = atom->f; + double *eps = atom->epsilon; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) { + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; + + compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d); + compute_drho1d(dx,dy,dz, order, drho_coeff, drho1d); + + u = ekx = eky = ekz = ZEROF; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + u += rho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx]; + ekx += drho1d[0][l]*rho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx]; + eky += rho1d[0][l]*drho1d[1][m]*rho1d[2][n]*u_brick[mz][my][mx]; + ekz += rho1d[0][l]*rho1d[1][m]*drho1d[2][n]*u_brick[mz][my][mx]; + } + } + } + ekx *= hx_inv; + eky *= hy_inv; + ekz *= hz_inv; + + // electrical potential + + if (potflag) phi[i] = u; + + // convert E-field to force and substract self forces + const double qfactor = qqrd2e * scale; + double qtmp = eps[i]*q[i]; + + s1 = x[i][0]*hx_inv; + s2 = x[i][1]*hy_inv; + s3 = x[i][2]*hz_inv; + sf = sf_coeff[0]*sin(2*MY_PI*s1); + sf += sf_coeff[1]*sin(4*MY_PI*s1); + sf *= 2*q[i]*q[i]; + f[i][0] += qfactor*(ekx*q[i] - sf); + + sf = sf_coeff[2]*sin(2*MY_PI*s2); + sf += sf_coeff[3]*sin(4*MY_PI*s2); + sf *= 2*q[i]*q[i]; + f[i][1] += qfactor*(eky*q[i] - sf); + + sf = sf_coeff[4]*sin(2*MY_PI*s3); + sf += sf_coeff[5]*sin(4*MY_PI*s3); + sf *= 2*q[i]*q[i]; + if (slabflag != 2) f[i][2] += qfactor*(ekz*q[i] - sf); + } +} + +/* ---------------------------------------------------------------------- + interpolate from grid to get electric field & force on my particles +------------------------------------------------------------------------- */ + +void PPPMDispDielectric::fieldforce_c_peratom() +{ + int i,l,m,n,nx,ny,nz,mx,my,mz; + FFT_SCALAR dx,dy,dz,x0,y0,z0; + FFT_SCALAR u_pa,v0,v1,v2,v3,v4,v5; + + // loop over my charges, interpolate electric field from nearby grid points + // (nx,ny,nz) = global coords of grid pt to "lower left" of charge + // (dx,dy,dz) = distance to "lower left" grid pt + // (mx,my,mz) = global coords of moving stencil pt + // ek = 3 components of E-field on particle + + double *q = atom->q; + double **x = atom->x; + + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++) { + nx = part2grid[i][0]; + ny = part2grid[i][1]; + nz = part2grid[i][2]; + dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv; + dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv; + dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv; + + compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d); + + u_pa = v0 = v1 = v2 = v3 = v4 = v5 = ZEROF; + for (n = nlower; n <= nupper; n++) { + mz = n+nz; + z0 = rho1d[2][n]; + for (m = nlower; m <= nupper; m++) { + my = m+ny; + y0 = z0*rho1d[1][m]; + for (l = nlower; l <= nupper; l++) { + mx = l+nx; + x0 = y0*rho1d[0][l]; + if (eflag_atom) u_pa += x0*u_brick[mz][my][mx]; + if (vflag_atom) { + v0 += x0*v0_brick[mz][my][mx]; + v1 += x0*v1_brick[mz][my][mx]; + v2 += x0*v2_brick[mz][my][mx]; + v3 += x0*v3_brick[mz][my][mx]; + v4 += x0*v4_brick[mz][my][mx]; + v5 += x0*v5_brick[mz][my][mx]; + } + } + } + } + + // electrostatic potential + + phi[i] = u_pa; + + // convert E-field to force + + const double qfactor = 0.5*force->qqrd2e * scale * q[i]; + + if (eflag_atom) eatom[i] += u_pa*qfactor; + if (vflag_atom) { + vatom[i][0] += v0*qfactor; + vatom[i][1] += v1*qfactor; + vatom[i][2] += v2*qfactor; + vatom[i][3] += v3*qfactor; + vatom[i][4] += v4*qfactor; + vatom[i][5] += v5*qfactor; + } + } +} + +/* ---------------------------------------------------------------------- + Slab-geometry correction term to dampen inter-slab interactions between + periodically repeating slabs. Yields good approximation to 2D Ewald if + adequate empty space is left between repeating slabs (J. Chem. Phys. + 111, 3155). Slabs defined here to be parallel to the xy plane. Also + extended to non-neutral systems (J. Chem. Phys. 131, 094107). +------------------------------------------------------------------------- */ + +void PPPMDispDielectric::slabcorr(int eflag) +{ + // compute local contribution to global dipole moment + + double *q = atom->q; + double **x = atom->x; + double *eps = atom->epsilon; + double zprd = domain->zprd; + int nlocal = atom->nlocal; + + double dipole = 0.0; + for (int i = 0; i < nlocal; i++) dipole += q[i]*x[i][2]; + + if (mu_flag) { + double **mu = atom->mu; + for (int i = 0; i < nlocal; i++) dipole += mu[i][2]; + } + + // sum local contributions to get global dipole moment + + double dipole_all; + MPI_Allreduce(&dipole,&dipole_all,1,MPI_DOUBLE,MPI_SUM,world); + + // need to make non-neutral systems and/or + // per-atom energy translationally invariant + + double dipole_r2 = 0.0; + if (eflag_atom || fabs(qsum) > SMALL) { + + if (mu_flag) + error->all(FLERR,"Cannot (yet) use kspace slab correction with " + "long-range dipoles and non-neutral systems or per-atom energy"); + + for (int i = 0; i < nlocal; i++) + dipole_r2 += q[i]*x[i][2]*x[i][2]; + + // sum local contributions + + double tmp; + MPI_Allreduce(&dipole_r2,&tmp,1,MPI_DOUBLE,MPI_SUM,world); + dipole_r2 = tmp; + } + + // compute corrections + + const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - + qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + const double qscale = qqrd2e * scale; + + if (eflag_global) energy += qscale * e_slabcorr; + + // per-atom energy + + if (eflag_atom) { + double efact = qscale * MY_2PI/volume; + for (int i = 0; i < nlocal; i++) + eatom[i] += efact * eps[i]*q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + + qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + } + + // add on force corrections + + double ffact = qscale * (-4.0*MY_PI/volume); + double **f = atom->f; + + for (int i = 0; i < nlocal; i++) { + f[i][2] += ffact * eps[i]*q[i]*(dipole_all - qsum*x[i][2]); + efield[i][2] += ffact * eps[i]*(dipole_all - qsum*x[i][2]); + } + + // add on torque corrections + + if (mu_flag && atom->torque) { + double **mu = atom->mu; + double **torque = atom->torque; + for (int i = 0; i < nlocal; i++) { + torque[i][0] += ffact * dipole_all * mu[i][1]; + torque[i][1] += -ffact * dipole_all * mu[i][0]; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local arrays +------------------------------------------------------------------------- */ + +double PPPMDispDielectric::memory_usage() +{ + double bytes = PPPMDisp::memory_usage(); + bytes += nmax*3 * sizeof(double); + bytes += nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + compute qsum,qsqsum,q2 and give error/warning if not charge neutral + called initially, when particle count changes, when charges are changed +------------------------------------------------------------------------- */ + +void PPPMDispDielectric::qsum_qsq() +{ + const double * const q = atom->q; + const double * const eps = atom->epsilon; + const int nlocal = atom->nlocal; + double qsum_local(0.0), qsqsum_local(0.0); + +#if defined(_OPENMP) +#pragma omp parallel for default(none) reduction(+:qsum_local,qsqsum_local) +#endif + for (int i = 0; i < nlocal; i++) { + double qtmp = eps[i]*q[i]; + qsum_local += qtmp; + qsqsum_local += qtmp*qtmp; + } + + MPI_Allreduce(&qsum_local,&qsum,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world); + + q2 = qsqsum * force->qqrd2e; +} diff --git a/src/USER-DIELECTRIC/pppm_disp_dielectric.h b/src/USER-DIELECTRIC/pppm_disp_dielectric.h new file mode 100644 index 0000000000..86623b91c4 --- /dev/null +++ b/src/USER-DIELECTRIC/pppm_disp_dielectric.h @@ -0,0 +1,62 @@ +/* -*- 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 KSPACE_CLASS + +KSpaceStyle(pppm/disp/dielectric,PPPMDispDielectric) + +#else + +#ifndef LMP_PPPM_DISP_DIELECTRIC_H +#define LMP_PPPM_DISP_DIELECTRIC_H + +#include "pppm_disp.h" + +namespace LAMMPS_NS { + +class PPPMDispDielectric : public PPPMDisp { + public: + PPPMDispDielectric(class LAMMPS *); + virtual ~PPPMDispDielectric(); + virtual double memory_usage(); + virtual void compute(int, int); + void qsum_qsq(); + void slabcorr(int); + + double** efield; + double* phi; + int potflag; // 1/0 if per-atom electrostatic potential phi is needed + + protected: + virtual void fieldforce_c_ik(); + virtual void fieldforce_c_ad(); + virtual void fieldforce_c_peratom(); + + class AtomVecDielectric* avec; + int mu_flag; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +*/