diff --git a/source/lmp/fix_dplr.cpp b/source/lmp/fix_dplr.cpp index c0f0bb930c..6120b5069c 100644 --- a/source/lmp/fix_dplr.cpp +++ b/source/lmp/fix_dplr.cpp @@ -34,7 +34,11 @@ is_key (const string& input) FixDPLR::FixDPLR(LAMMPS *lmp, int narg, char **arg) - :Fix(lmp, narg, arg) + :Fix(lmp, narg, arg), + efield(3, 0.0), + efield_fsum(4, 0.0), + efield_fsum_all(4, 0.0), + efield_force_flag(0) { virial_flag = 1; @@ -54,6 +58,13 @@ FixDPLR::FixDPLR(LAMMPS *lmp, int narg, char **arg) model = string(arg[iarg+1]); iarg += 2; } + if (string(arg[iarg]) == string("efield")) { + if (iarg+3 > narg) error->all(FLERR,"Illegal fix adapt command, efield should be provided 3 float numbers"); + efield[0] = atof(arg[iarg+1]); + efield[1] = atof(arg[iarg+2]); + efield[2] = atof(arg[iarg+3]); + iarg += 4; + } if (string(arg[iarg]) == string("type_associate")) { int iend = iarg+1; while (iend < narg && (! is_key(arg[iend]) )) { @@ -105,6 +116,7 @@ FixDPLR::FixDPLR(LAMMPS *lmp, int narg, char **arg) int FixDPLR::setmask() { int mask = 0; + mask |= THERMO_ENERGY; mask |= POST_INTEGRATE; mask |= PRE_FORCE; mask |= POST_FORCE; @@ -347,6 +359,7 @@ void FixDPLR::post_force(int vflag) int nall = nlocal + nghost; vector dcoord(nall*3, 0.0), dbox(9, 0.0), dfele(nlocal*3, 0.0); vector dtype(nall, 0); + // set values for dcoord, dbox, dfele { int *type = atom->type; for (int ii = 0; ii < nall; ++ii){ @@ -366,9 +379,36 @@ void FixDPLR::post_force(int vflag) } } assert(dfele_.size() == nlocal * 3); + // revise force according to efield for (int ii = 0; ii < nlocal*3; ++ii){ dfele[ii] = dfele_[ii]; } + // revise force and virial according to efield + double * q = atom->q; + imageint *image = atom->image; + double unwrap[3]; + double v[6]; + efield_fsum[0] = efield_fsum[1] = efield_fsum[2] = efield_fsum[3] = 0.0; + efield_force_flag = 0; + for (int ii = 0; ii < nlocal; ++ii){ + for (int dd = 0; dd < 3; ++dd){ + dfele[ii*3+dd] += q[ii] * efield[dd]; + } + domain->unmap(x[ii],image[ii],unwrap); + efield_fsum[0] -= efield[0]*unwrap[0]+efield[1]*unwrap[1]+efield[2]*unwrap[2]; + efield_fsum[1] += efield[0]; + efield_fsum[2] += efield[1]; + efield_fsum[3] += efield[2]; + if (evflag) { + v[0] = q[ii] * efield[0] *unwrap[0]; + v[1] = q[ii] * efield[1] *unwrap[1]; + v[2] = q[ii] * efield[2] *unwrap[2]; + v[3] = q[ii] * efield[0] *unwrap[1]; + v[4] = q[ii] * efield[0] *unwrap[2]; + v[5] = q[ii] * efield[1] *unwrap[2]; + v_tally(ii, v); + } + } } // lmp nlist NeighList * list = pair_nnp->list; @@ -466,5 +506,28 @@ void FixDPLR::unpack_reverse_comm(int n, int *list, double *buf) } } +/* ---------------------------------------------------------------------- + return energy added by fix +------------------------------------------------------------------------- */ + +double FixDPLR::compute_scalar(void) +{ + if (efield_force_flag == 0) { + MPI_Allreduce(&efield_fsum[0],&efield_fsum_all[0],4,MPI_DOUBLE,MPI_SUM,world); + efield_force_flag = 1; + } + return efield_fsum_all[0]; +} +/* ---------------------------------------------------------------------- + return total extra force due to fix +------------------------------------------------------------------------- */ +double FixDPLR::compute_vector(int n) +{ + if (efield_force_flag == 0) { + MPI_Allreduce(&efield_fsum[0],&efield_fsum_all[0],4,MPI_DOUBLE,MPI_SUM,world); + efield_force_flag = 1; + } + return efield_fsum_all[n+1]; +} diff --git a/source/lmp/fix_dplr.h b/source/lmp/fix_dplr.h index edc9204874..a0f45185b4 100644 --- a/source/lmp/fix_dplr.h +++ b/source/lmp/fix_dplr.h @@ -32,6 +32,8 @@ namespace LAMMPS_NS { void post_force(int); int pack_reverse_comm(int, int, double *); void unpack_reverse_comm(int, int *, double *); + double compute_scalar(void); + double compute_vector(int); private: PairNNP * pair_nnp; DeepTensor dpt; @@ -45,6 +47,9 @@ namespace LAMMPS_NS { map bk_type_asso; vector dipole_recd; vector dfcorr_buff; + vector efield; + vector efield_fsum, efield_fsum_all; + int efield_force_flag; void get_valid_pairs(vector >& pairs); }; }