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

Bringing USER-DIELECTRIC up-to-date with latest changes in upstream LAMMPS

Browse Source
This commit is contained in:
Trung Nguyen
2021-05-28 12:41:52 -05:00
parent 2a33e3674e
commit 454e11f7a5
15 changed files with 255 additions and 2383 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
src/GPU/pair_lj_cut_gpu.cpp
View File

@ -135,29 +135,24 @@ void PairLJCutGPU::init_style()
{
cut_respa = nullptr;
//if (force->newton_pair)
// error->all(FLERR,"Cannot use newton pair with lj/cut/gpu pair style");
if (force->newton_pair)
error->all(FLERR,"Cannot use newton pair with lj/cut/gpu pair style");
// Repeat cutsq calculation because done after call to init_style
double maxcut = -1.0;
double cut;
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
cut = init_one(i,j);
if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0)) {
cut = init_one(i,j);
//printf("lj/cut/gpu: i = %d; j = %d: setflag = %d cut = %f\n", i, j, setflag[i][j], cut);
cut *= cut;
if (cut > maxcut)
maxcut = cut;
cutsq[i][j] = cutsq[j][i] = cut;
} else {
//printf("lj/cut/gpu: i = %d; j = %d: setflag = %d cut = %f\n", i, j, setflag[i][j], cut);
} else
cutsq[i][j] = cutsq[j][i] = 0.0;
}
}
}
double cell_size = sqrt(maxcut) + neighbor->skin;
int maxspecial=0;

1254
src/USER-DIELECTRIC/atom_vec_dielectric.cpp
View File

File diff suppressed because it is too large Load Diff

52
src/USER-DIELECTRIC/atom_vec_dielectric.h
View File

@ -20,54 +20,24 @@ AtomStyle(dielectric,AtomVecDielectric)
#ifndef LMP_ATOM_VEC_DIELECTRIC_H
#define LMP_ATOM_VEC_DIELECTRIC_H
#include "atom_vec_full.h"
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecDielectric : public AtomVecFull {
class AtomVecDielectric : public AtomVec {
public:
AtomVecDielectric(class LAMMPS *);
~AtomVecDielectric();
void grow(int);
void grow_reset();
void copy(int, int, int);
int pack_comm(int, int *, double *, int, int *);
int pack_comm_vel(int, int *, double *, int, int *);
int pack_comm_hybrid(int, int *, double *);
void unpack_comm(int, int, double *);
void unpack_comm_vel(int, int, double *);
int unpack_comm_hybrid(int, int, double *);
int pack_reverse(int, int, double *);
void unpack_reverse(int, int *, double *);
int pack_border(int, int *, double *, int, int *);
int pack_border_vel(int, int *, double *, int, int *);
int pack_border_hybrid(int, int *, double *);
void unpack_border(int, int, double *);
void unpack_border_vel(int, int, double *);
int unpack_border_hybrid(int, int, double *);
int pack_exchange(int, double *);
int unpack_exchange(double *);
int size_restart();
int pack_restart(int, double *);
int unpack_restart(double *);
void create_atom(int, double *);
void data_atom(double *, imageint, char **);
int data_atom_hybrid(int, char **);
void pack_data(double **);
int pack_data_hybrid(int, double *);
void write_data(FILE *, int, double **);
int write_data_hybrid(FILE *, double *);
bigint memory_usage();
void grow_pointers();
void create_atom_post(int);
void data_atom_post(int);
void pack_data_post(int);
double **mu;
double *area,*ed,*em,*epsilon,*curvature,*q_unscaled;
private:
public:
double **mu; // normal vector at the patch
double *area; // patch area
double *em; // mean dielectric constant at the patch
double *ed; // difference in dielectric constants at the patch
double *epsilon; // dielectric at the patch and real charges
double *curvature; // curvature at the patch
double *q_real; // unscaled charge: q_real = value read in from data file
};
}

61
src/USER-DIELECTRIC/fix_polarize_bem_icc.cpp
View File

@ -21,21 +21,21 @@
Barros, Sinkovits, Luijten, J. Chem. Phys 2014, 140, 064903
------------------------------------------------------------------------- */
#include <cstring>
#include <cstdlib>
#include "fix_polarize_bem_icc.h"
#include "atom_vec_dielectric.h"
#include "update.h"
#include "atom.h"
#include "atom_vec_dielectric.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "kspace.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "math_const.h"
#include "msm_dielectric.h"
#include "neighbor.h"
#include "pair_coul_cut_dielectric.h"
#include "pair_coul_long_dielectric.h"
#include "pair_lj_cut_coul_cut_dielectric.h"
@ -45,7 +45,10 @@
#include "msm_dielectric.h"
#include "random_park.h"
#include "timer.h"
#include "error.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -65,9 +68,9 @@ FixPolarizeBEMICC::FixPolarizeBEMICC(LAMMPS *lmp, int narg, char **arg) :
// parse required arguments
nevery = force->inumeric(FLERR,arg[3]);
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery < 0) error->all(FLERR,"Illegal fix polarize/bem/icc command");
double tol = force->numeric(FLERR,arg[4]);
double tol = utils::numeric(FLERR,arg[4],false,lmp);
tol_abs = tol_rel = tol;
itr_max = 20;
@ -214,12 +217,12 @@ void FixPolarizeBEMICC::pre_force(int)
void FixPolarizeBEMICC::compute_induced_charges()
{
double *q = atom->q;
double *q_real = avec->q_real;
double **norm = avec->mu;
double *area = avec->area;
double *ed = avec->ed;
double *em = avec->em;
double *epsilon = avec->epsilon;
double *q_real = atom->q_unscaled;
double **norm = atom->mu;
double *area = atom->area;
double *ed = atom->ed;
double *em = atom->em;
double *epsilon = atom->epsilon;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double epsilon0 = force->dielectric;
@ -354,11 +357,11 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
while (iarg < narg) {
if (strcmp(arg[iarg],"itr_max") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
itr_max = force->numeric(FLERR,arg[iarg+1]);
itr_max = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"omega") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
omega = force->numeric(FLERR,arg[iarg+1]);
omega = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix_modify command");
@ -371,12 +374,12 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
double epsiloni=-1, areai=-1;
double qreali=0;
int set_charge=0;
double ediff = force->numeric(FLERR,arg[iarg+1]);
double emean = force->numeric(FLERR,arg[iarg+2]);
if (strcmp(arg[iarg+3],"NULL") != 0) epsiloni = force->numeric(FLERR,arg[iarg+3]);
if (strcmp(arg[iarg+4],"NULL") != 0) areai = force->numeric(FLERR,arg[iarg+4]);
double ediff = utils::numeric(FLERR,arg[iarg+1],false,lmp);
double emean = utils::numeric(FLERR,arg[iarg+2],false,lmp);
if (strcmp(arg[iarg+3],"NULL") != 0) epsiloni = utils::numeric(FLERR,arg[iarg+3],false,lmp);
if (strcmp(arg[iarg+4],"NULL") != 0) areai = utils::numeric(FLERR,arg[iarg+4],false,lmp);
if (strcmp(arg[iarg+5],"NULL") != 0) {
qreali = force->numeric(FLERR,arg[iarg+5]);
qreali = utils::numeric(FLERR,arg[iarg+5],false,lmp);
set_charge = 1;
}
set_dielectric_params(ediff, emean, epsiloni, areai, set_charge, qreali);
@ -384,8 +387,8 @@ int FixPolarizeBEMICC::modify_param(int narg, char **arg)
iarg += 6;
} else if (strcmp(arg[iarg],"rand") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal fix_modify command");
ave_charge = force->numeric(FLERR,arg[iarg+1]);
seed_charge = force->numeric(FLERR,arg[iarg+2]);
ave_charge = utils::numeric(FLERR,arg[iarg+1],false,lmp);
seed_charge = utils::numeric(FLERR,arg[iarg+2],false,lmp);
randomized = 1;
iarg += 3;
} else error->all(FLERR,"Illegal fix_modify command");
@ -419,11 +422,11 @@ void FixPolarizeBEMICC::unpack_forward_comm(int n, int first, double *buf)
void FixPolarizeBEMICC::set_dielectric_params(double ediff, double emean,
double epsiloni, double areai, int set_charge, double qreali)
{
double *area = avec->area;
double *ed = avec->ed;
double *em = avec->em;
double *q_real = avec->q_real;
double *epsilon = avec->epsilon;
double *area = atom->area;
double *ed = atom->ed;
double *em = atom->em;
double *q_real = atom->q_unscaled;
double *epsilon = atom->epsilon;
int *mask = atom->mask;
int nlocal = atom->nlocal;

78
src/USER-DIELECTRIC/msm_dielectric.cpp
View File

@ -15,24 +15,22 @@
Contributing authors: Trung Nguyen (Northwestern)
------------------------------------------------------------------------- */
#include <mpi.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "msm_dielectric.h"
#include "atom.h"
#include "atom_vec_dielectric.h"
#include "comm.h"
#include "gridcomm.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "force.h"
#include "gridcomm.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "pair.h"
#include <cstring>
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -86,7 +84,7 @@ void MSMDielectric::compute(int eflag, int vflag)
if (scalar_pressure_flag && vflag_either) {
if (vflag_atom)
error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' to obtain "
"per-atom virial with kspace_style MSMDielectric");
"per-atom virial with kspace_style msm/dielectric");
// must switch on global energy computation if not already on
@ -109,16 +107,7 @@ void MSMDielectric::compute(int eflag, int vflag)
// invoke allocate_peratom() if needed for first time
if (vflag_atom && !peratom_allocate_flag) {
allocate_peratom();
cg_peratom_all->ghost_notify();
cg_peratom_all->setup();
for (int n=0; n<levels; n++) {
if (!active_flag[n]) continue;
cg_peratom[n]->ghost_notify();
cg_peratom[n]->setup();
}
}
if (vflag_atom && !peratom_allocate_flag) allocate_peratom();
// convert atoms from box to lamda coords
@ -146,7 +135,8 @@ void MSMDielectric::compute(int eflag, int vflag)
// to fully sum contribution in their 3d grid
current_level = 0;
cg_all->reverse_comm(this,REVERSE_RHO);
gcall->reverse_comm_kspace(this,1,sizeof(double),REVERSE_RHO,
gcall_buf1,gcall_buf2,MPI_DOUBLE);
// forward communicate charge density values to fill ghost grid points
// compute direct sum interaction and then restrict to coarser grid
@ -154,23 +144,30 @@ void MSMDielectric::compute(int eflag, int vflag)
for (int n=0; n<=levels-2; n++) {
if (!active_flag[n]) continue;
current_level = n;
cg[n]->forward_comm(this,FORWARD_RHO);
gc[n]->forward_comm_kspace(this,1,sizeof(double),FORWARD_RHO,
gc_buf1[n],gc_buf2[n],MPI_DOUBLE);
direct(n);
restriction(n);
}
// compute direct interation for top grid level for nonperiodic
// compute direct interation for top grid level for non-periodic
// and for second from top grid level for periodic
if (active_flag[levels-1]) {
if (domain->nonperiodic) {
current_level = levels-1;
cg[levels-1]->forward_comm(this,FORWARD_RHO);
gc[levels-1]->
forward_comm_kspace(this,1,sizeof(double),FORWARD_RHO,
gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
direct_top(levels-1);
cg[levels-1]->reverse_comm(this,REVERSE_AD);
gc[levels-1]->
reverse_comm_kspace(this,1,sizeof(double),REVERSE_AD,
gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
if (vflag_atom)
cg_peratom[levels-1]->reverse_comm(this,REVERSE_AD_PERATOM);
gc[levels-1]->
reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM,
gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
} else {
// Here using MPI_Allreduce is cheaper than using commgrid
grid_swap_forward(levels-1,qgrid[levels-1]);
@ -178,7 +175,9 @@ void MSMDielectric::compute(int eflag, int vflag)
grid_swap_reverse(levels-1,egrid[levels-1]);
current_level = levels-1;
if (vflag_atom)
cg_peratom[levels-1]->reverse_comm(this,REVERSE_AD_PERATOM);
gc[levels-1]->
reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM,
gc_buf1[levels-1],gc_buf2[levels-1],MPI_DOUBLE);
}
}
@ -190,24 +189,28 @@ void MSMDielectric::compute(int eflag, int vflag)
prolongation(n);
current_level = n;
cg[n]->reverse_comm(this,REVERSE_AD);
gc[n]->reverse_comm_kspace(this,1,sizeof(double),REVERSE_AD,
gc_buf1[n],gc_buf2[n],MPI_DOUBLE);
// extra per-atom virial communication
if (vflag_atom)
cg_peratom[n]->reverse_comm(this,REVERSE_AD_PERATOM);
gc[n]->reverse_comm_kspace(this,6,sizeof(double),REVERSE_AD_PERATOM,
gc_buf1[n],gc_buf2[n],MPI_DOUBLE);
}
// all procs communicate E-field values
// to fill ghost cells surrounding their 3d bricks
current_level = 0;
cg_all->forward_comm(this,FORWARD_AD);
gcall->forward_comm_kspace(this,1,sizeof(double),FORWARD_AD,
gcall_buf1,gcall_buf2,MPI_DOUBLE);
// extra per-atom energy/virial communication
if (vflag_atom)
cg_peratom_all->forward_comm(this,FORWARD_AD_PERATOM);
gcall->forward_comm_kspace(this,6,sizeof(double),FORWARD_AD_PERATOM,
gcall_buf1,gcall_buf2,MPI_DOUBLE);
// calculate the force on my particles (interpolation)
@ -266,8 +269,7 @@ void MSMDielectric::compute(int eflag, int vflag)
// convert atoms back from lamda to box coords
if (triclinic)
domain->lamda2x(atom->nlocal);
if (triclinic) domain->lamda2x(atom->nlocal);
}
/* ----------------------------------------------------------------------
@ -276,8 +278,6 @@ void MSMDielectric::compute(int eflag, int vflag)
void MSMDielectric::fieldforce()
{
//fprintf(screen,"MSM interpolation\n\n");
double ***egridn = egrid[0];
int i,l,m,n,nx,ny,nz,mx,my,mz;
@ -296,7 +296,7 @@ void MSMDielectric::fieldforce()
double *q = atom->q;
double **x = atom->x;
double **f = atom->f;
double *eps = avec->epsilon;
double *eps = atom->epsilon;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {

28
src/USER-DIELECTRIC/pair_coul_cut_dielectric.cpp
View File

@ -74,7 +74,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
double *q_real = avec->q_real;
double *q_real = avec->q_unscaled;
double* eps = avec->epsilon;
double** norm = avec->mu;
double* curvature = avec->curvature;
@ -128,7 +128,7 @@ void PairCoulCutDielectric::compute(int eflag, int vflag)
if (rsq < cutsq[itype][jtype] && rsq > EPSILON) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
efield_i = qqrd2e * scale[itype][jtype] * q[j]*rinv;
efield_i = scale[itype][jtype] * q[j]*rinv;
forcecoul = qtmp*efield_i;
fpair_i = factor_coul*etmp*forcecoul*r2inv;
@ -175,3 +175,27 @@ void PairCoulCutDielectric::init_style()
neighbor->requests[irequest]->full = 1;
}
/* ---------------------------------------------------------------------- */
double PairCoulCutDielectric::single(int i, int j, int itype, int jtype,
double rsq,
double factor_coul, double factor_lj,
double &fforce)
{
double r2inv,forcecoul,phicoul,ei,ej;
double* eps = avec->epsilon;
r2inv = 1.0/rsq;
forcecoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv)*eps[i];
double eng = 0.0;
if (eps[i] == 1) ei = 0;
else ei = eps[i];
if (eps[j] == 1) ej = 0;
else ej = eps[j];
phicoul = force->qqrd2e * atom->q[i]*atom->q[j]*sqrt(r2inv);
phicoul *= 0.5*(ei+ej);
eng += factor_coul*phicoul;
return eng;
}

1
src/USER-DIELECTRIC/pair_coul_cut_dielectric.h
View File

@ -29,6 +29,7 @@ class PairCoulCutDielectric : public PairCoulCut {
PairCoulCutDielectric(class LAMMPS *);
virtual ~PairCoulCutDielectric();
virtual void compute(int, int);
virtual double single(int, int, int, int, double, double, double, double &);
void init_style();
double** efield;

21
src/USER-DIELECTRIC/pair_coul_long_dielectric.cpp
View File

@ -15,25 +15,22 @@
Contributing author: Trung Nguyen (Northwestern)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_coul_long_dielectric.h"
#include "atom.h"
#include "atom_vec_dielectric.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "kspace.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -87,10 +84,10 @@ void PairCoulLongDielectric::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
double** norm = avec->mu;
double* curvature = avec->curvature;
double* area = avec->area;
double *eps = avec->epsilon;
double** norm = atom->mu;
double* curvature = atom->curvature;
double* area = atom->area;
double *eps = atom->epsilon;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;

20
src/USER-DIELECTRIC/pair_lj_cut_coul_cut_dielectric.cpp
View File

@ -15,21 +15,21 @@
Contributing author: Trung Nguyen (Northwestern)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_lj_cut_coul_cut_dielectric.h"
#include "atom.h"
#include "atom_vec_dielectric.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -78,11 +78,11 @@ void PairLJCutCoulCutDielectric::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
double *q_real = avec->q_real;
double* eps = avec->epsilon;
double** norm = avec->mu;
double* curvature = avec->curvature;
double* area = avec->area;
double *q_real = atom->q_unscaled;
double* eps = atom->epsilon;
double** norm = atom->mu;
double* curvature = atom->curvature;
double* area = atom->area;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;

24
src/USER-DIELECTRIC/pair_lj_cut_coul_long_dielectric.cpp
View File

@ -15,25 +15,25 @@
Contributing author: Trung Nguyen (Northwestern)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_lj_cut_coul_long_dielectric.h"
#include "atom.h"
#include "atom_vec_dielectric.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "kspace.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "kspace.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -95,10 +95,10 @@ void PairLJCutCoulLongDielectric::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
double *eps = avec->epsilon;
double** norm = avec->mu;
double* curvature = avec->curvature;
double* area = avec->area;
double *eps = atom->epsilon;
double** norm = atom->mu;
double* curvature = atom->curvature;
double* area = atom->area;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;

22
src/USER-DIELECTRIC/pair_lj_cut_coul_msm_dielectric.cpp
View File

@ -15,25 +15,25 @@
Contributing author: Trung Nguyen (Northwestern)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_lj_cut_coul_msm_dielectric.h"
#include "atom.h"
#include "atom_vec_dielectric.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "kspace.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -108,10 +108,10 @@ void PairLJCutCoulMSMDielectric::compute(int eflag, int vflag)
double **x = atom->x;
double **f = atom->f;
double *q = atom->q;
double *eps = avec->epsilon;
double** norm = avec->mu;
double* curvature = avec->curvature;
double* area = avec->area;
double *eps = atom->epsilon;
double** norm = atom->mu;
double* curvature = atom->curvature;
double* area = atom->area;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;

1017
src/USER-DIELECTRIC/pppm_dielectric.cpp
View File

File diff suppressed because it is too large Load Diff

32
src/USER-DIELECTRIC/pppm_dielectric.h
View File

@ -28,10 +28,7 @@ class PPPMDielectric : public PPPM {
public:
PPPMDielectric(class LAMMPS *);
virtual ~PPPMDielectric();
virtual void init();
void setup_grid();
virtual void compute(int, int);
virtual double memory_usage();
double** efield;
double* phi;
@ -40,40 +37,11 @@ class PPPMDielectric : public PPPM {
void qsum_qsq();
protected:
FFT_SCALAR ***densityx_brick_dipole,***densityy_brick_dipole,***densityz_brick_dipole;
FFT_SCALAR ***vdxx_brick_dipole,***vdyy_brick_dipole,***vdzz_brick_dipole;
FFT_SCALAR ***vdxy_brick_dipole,***vdxz_brick_dipole,***vdyz_brick_dipole;
FFT_SCALAR ***u_brick_dipole;
FFT_SCALAR ***ux_brick_dipole,***uy_brick_dipole,***uz_brick_dipole;
FFT_SCALAR *densityx_fft_dipole,*densityy_fft_dipole,*densityz_fft_dipole;
FFT_SCALAR *work3,*work4;
class GridComm *cg_mu;
virtual void allocate();
virtual void deallocate();
void slabcorr();
void fieldforce_ik();
void fieldforce_ad();
// grid communication
virtual void pack_forward(int, FFT_SCALAR *, int, int *);
virtual void unpack_forward(int, FFT_SCALAR *, int, int *);
virtual void pack_reverse(int, FFT_SCALAR *, int, int *);
virtual void unpack_reverse(int, FFT_SCALAR *, int, int *);
// dipole
int mu_flag;
double musqsum,musum,mu2;
void make_rho_dipole();
void brick2fft_dipole();
void poisson_ik_dipole();
void fieldforce_ik_dipole();
void musum_musq();
class AtomVecDielectric* avec;
};

13
src/atom.cpp
View File

@ -197,6 +197,10 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp)
rho = drho = esph = desph = cv = nullptr;
vest = nullptr;
// USER-DIELECTRIC package
area = ed = em = epsilon = curvature = q_unscaled = nullptr;
// end of customization section
// --------------------------------------------------------------------
@ -509,6 +513,15 @@ void Atom::peratom_create()
add_peratom("eff_plastic_strain_rate",&eff_plastic_strain_rate,DOUBLE,0);
add_peratom("damage",&damage,DOUBLE,0);
// USER-DIELECTRIC package
add_peratom("area",&area,DOUBLE,0);
add_peratom("ed",&ed,DOUBLE,0);
add_peratom("em",&em,DOUBLE,0);
add_peratom("epsilon",&epsilon,DOUBLE,0);
add_peratom("curvature",&curvature,DOUBLE,0);
add_peratom("q_unscaled",&curvature,DOUBLE,0);
// end of customization section
// --------------------------------------------------------------------
}

4
src/atom.h
View File

@ -157,6 +157,10 @@ class Atom : protected Pointers {
double *rho,*drho,*esph,*desph,*cv;
double **vest;
// USER-DIELECTRIC package
double *area,*ed,*em,*epsilon,*curvature,*q_unscaled;
// end of customization section
// --------------------------------------------------------------------