ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

potflag and phi are not used

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-10-10 05:35:51 -04:00
parent 482a6632e9
commit 948ee80169
6 changed files with 12 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
src/DIELECTRIC/msm_dielectric.cpp
View File

@ -31,11 +31,7 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
MSMDielectric::MSMDielectric(LAMMPS *_lmp) : MSM(_lmp)
{
efield = nullptr;
phi = nullptr;
}
MSMDielectric::MSMDielectric(LAMMPS *_lmp) : MSM(_lmp), efield(nullptr) {}
/* ----------------------------------------------------------------------
free all memory
@ -44,7 +40,6 @@ MSMDielectric::MSMDielectric(LAMMPS *_lmp) : MSM(_lmp)
MSMDielectric::~MSMDielectric()
{
memory->destroy(efield);
memory->destroy(phi);
}
/* ----------------------------------------------------------------------
@ -111,11 +106,9 @@ void MSMDielectric::compute(int eflag, int vflag)
if (atom->nmax > nmax) {
memory->destroy(part2grid);
memory->destroy(efield);
memory->destroy(phi);
nmax = atom->nmax;
memory->create(part2grid,nmax,3,"msm:part2grid");
memory->create(efield,nmax,3,"msm:efield");
memory->create(phi,nmax,"msm:phi");
}
// find grid points for all my particles
@ -276,7 +269,7 @@ void MSMDielectric::fieldforce()
int i,l,m,n,nx,ny,nz,mx,my,mz;
double dx,dy,dz;
double phi_x,phi_y,phi_z,u;
double phi_x,phi_y,phi_z;
double dphi_x,dphi_y,dphi_z;
double ekx,eky,ekz,etmp;
@ -303,7 +296,7 @@ void MSMDielectric::fieldforce()
compute_phis_and_dphis(dx,dy,dz);
u = ekx = eky = ekz = 0.0;
ekx = eky = ekz = 0.0;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
phi_z = phi1d[2][n];
@ -317,7 +310,6 @@ void MSMDielectric::fieldforce()
phi_x = phi1d[0][l];
dphi_x = dphi1d[0][l];
etmp = egridn[mz][my][mx];
u += phi_z*phi_y*phi_x*etmp;
ekx += dphi_x*phi_y*phi_z*etmp;
eky += phi_x*dphi_y*phi_z*etmp;
ekz += phi_x*phi_y*dphi_z*etmp;
@ -329,10 +321,6 @@ void MSMDielectric::fieldforce()
eky *= delyinv[0];
ekz *= delzinv[0];
// electrical potential
phi[i] = u;
// effectively divide by length for a triclinic system
if (triclinic) {

1
src/DIELECTRIC/msm_dielectric.h
View File

@ -32,7 +32,6 @@ class MSMDielectric : public MSM {
void compute(int, int) override;
double **efield;
double *phi;
protected:
void fieldforce() override;

26
src/DIELECTRIC/pppm_dielectric.cpp
View File

@ -40,13 +40,10 @@ static constexpr FFT_SCALAR ZEROF = 0.0;
/* ---------------------------------------------------------------------- */
PPPMDielectric::PPPMDielectric(LAMMPS *_lmp) : PPPM(_lmp)
PPPMDielectric::PPPMDielectric(LAMMPS *_lmp) : PPPM(_lmp), efield(nullptr)
{
group_group_enable = 0;
efield = nullptr;
phi = nullptr;
potflag = 0;
use_qscaled = true;
// no warnings about non-neutral systems from qsum_qsq()
@ -61,7 +58,6 @@ PPPMDielectric::PPPMDielectric(LAMMPS *_lmp) : PPPM(_lmp)
PPPMDielectric::~PPPMDielectric()
{
memory->destroy(efield);
memory->destroy(phi);
}
/* ----------------------------------------------------------------------
@ -107,11 +103,9 @@ void PPPMDielectric::compute(int eflag, int vflag)
if (atom->nmax > nmax) {
memory->destroy(part2grid);
memory->destroy(efield);
memory->destroy(phi);
nmax = atom->nmax;
memory->create(part2grid,nmax,3,"pppm/dielectric:part2grid");
memory->create(efield,nmax,3,"pppm/dielectric:efield");
memory->create(phi,nmax,"pppm/dielectric:phi");
}
// find grid points for all my particles
@ -384,7 +378,7 @@ void PPPMDielectric::fieldforce_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;
FFT_SCALAR ekx,eky,ekz;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
@ -408,7 +402,7 @@ void PPPMDielectric::fieldforce_ik()
compute_rho1d(dx,dy,dz);
u = ekx = eky = ekz = ZEROF;
ekx = eky = ekz = ZEROF;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
z0 = rho1d[2][n];
@ -418,7 +412,6 @@ void PPPMDielectric::fieldforce_ik()
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];
@ -426,10 +419,6 @@ void PPPMDielectric::fieldforce_ik()
}
}
// electrostatic potential
if (potflag) phi[i] = u;
// convert E-field to force
const double efactor = scale * eps[i];
efield[i][0] = efactor*ekx;
@ -451,7 +440,7 @@ void PPPMDielectric::fieldforce_ad()
{
int i,l,m,n,nx,ny,nz,mx,my,mz;
FFT_SCALAR dx,dy,dz;
FFT_SCALAR ekx,eky,ekz,u;
FFT_SCALAR ekx,eky,ekz;
double s1,s2,s3;
double sf = 0.0;
@ -489,14 +478,13 @@ void PPPMDielectric::fieldforce_ad()
compute_rho1d(dx,dy,dz);
compute_drho1d(dx,dy,dz);
u = ekx = eky = ekz = ZEROF;
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];
@ -507,10 +495,6 @@ void PPPMDielectric::fieldforce_ad()
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 * eps[i];

2
src/DIELECTRIC/pppm_dielectric.h
View File

@ -31,8 +31,6 @@ class PPPMDielectric : public PPPM {
void compute(int, int) override;
double **efield;
double *phi;
int potflag; // 1/0 if per-atom electrostatic potential phi is needed
protected:
void slabcorr() override;

29
src/DIELECTRIC/pppm_disp_dielectric.cpp
View File

@ -40,7 +40,7 @@ static constexpr FFT_SCALAR ZEROF = 0.0;
/* ---------------------------------------------------------------------- */
PPPMDispDielectric::PPPMDispDielectric(LAMMPS *_lmp) : PPPMDisp(_lmp)
PPPMDispDielectric::PPPMDispDielectric(LAMMPS *_lmp) : PPPMDisp(_lmp), efield(nullptr)
{
dipoleflag = 0; // turned off for now, until dipole works
group_group_enable = 0;
@ -51,10 +51,6 @@ PPPMDispDielectric::PPPMDispDielectric(LAMMPS *_lmp) : PPPMDisp(_lmp)
// no warnings about non-neutral systems from qsum_qsq()
warn_nonneutral = 2;
efield = nullptr;
phi = nullptr;
potflag = 0;
avec = dynamic_cast<AtomVecDielectric *>(atom->style_match("dielectric"));
if (!avec) error->all(FLERR,"pppm/dielectric requires atom style dielectric");
}
@ -64,7 +60,6 @@ PPPMDispDielectric::PPPMDispDielectric(LAMMPS *_lmp) : PPPMDisp(_lmp)
PPPMDispDielectric::~PPPMDispDielectric()
{
memory->destroy(efield);
memory->destroy(phi);
}
/* ----------------------------------------------------------------------
@ -97,14 +92,12 @@ void PPPMDispDielectric::compute(int eflag, int vflag)
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");
@ -657,7 +650,7 @@ 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;
FFT_SCALAR ekx,eky,ekz;
// loop over my charges, interpolate electric field from nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
@ -681,7 +674,7 @@ void PPPMDispDielectric::fieldforce_c_ik()
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
u = ekx = eky = ekz = ZEROF;
ekx = eky = ekz = ZEROF;
for (n = nlower; n <= nupper; n++) {
mz = n+nz;
z0 = rho1d[2][n];
@ -691,7 +684,6 @@ void PPPMDispDielectric::fieldforce_c_ik()
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];
@ -699,10 +691,6 @@ void PPPMDispDielectric::fieldforce_c_ik()
}
}
// electrostatic potential
if (potflag) phi[i] = u;
// convert E-field to force
const double efactor = scale * eps[i];
@ -767,14 +755,13 @@ void PPPMDispDielectric::fieldforce_c_ad()
compute_rho1d(dx,dy,dz, order, rho_coeff, rho1d);
compute_drho1d(dx,dy,dz, order, drho_coeff, drho1d);
u = ekx = eky = ekz = ZEROF;
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];
@ -785,10 +772,6 @@ void PPPMDispDielectric::fieldforce_c_ad()
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;
@ -866,10 +849,6 @@ void PPPMDispDielectric::fieldforce_c_peratom()
}
}
// electrostatic potential
phi[i] = u_pa;
// convert E-field to force
const double qfactor = 0.5*force->qqrd2e * scale * q[i];

2
src/DIELECTRIC/pppm_disp_dielectric.h
View File

@ -33,8 +33,6 @@ class PPPMDispDielectric : public PPPMDisp {
void slabcorr(int) override;
double **efield;
double *phi;
int potflag; // 1/0 if per-atom electrostatic potential phi is needed
protected:
void make_rho_c() override;