From 14091c39531fb98a1b7c10a6bff25546ef711c03 Mon Sep 17 00:00:00 2001 From: Han Wang Date: Sat, 23 May 2020 17:20:50 +0800 Subject: [PATCH 1/2] implement external efield in fix/dplr --- source/lmp/fix_dplr.cpp | 31 ++++++++++++++++++++++++++++++- source/lmp/fix_dplr.h | 1 + 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/source/lmp/fix_dplr.cpp b/source/lmp/fix_dplr.cpp index c0f0bb930c..1119383129 100644 --- a/source/lmp/fix_dplr.cpp +++ b/source/lmp/fix_dplr.cpp @@ -34,7 +34,7 @@ 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) { virial_flag = 1; @@ -54,6 +54,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]) )) { @@ -347,6 +354,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 +374,30 @@ 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]; + 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); + 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; diff --git a/source/lmp/fix_dplr.h b/source/lmp/fix_dplr.h index edc9204874..094b0aa3e1 100644 --- a/source/lmp/fix_dplr.h +++ b/source/lmp/fix_dplr.h @@ -45,6 +45,7 @@ namespace LAMMPS_NS { map bk_type_asso; vector dipole_recd; vector dfcorr_buff; + vector efield; void get_valid_pairs(vector >& pairs); }; } From 2362a2fbec23068a79ba49c14a536c4cc753c1dd Mon Sep 17 00:00:00 2001 From: Han Wang Date: Mon, 10 Aug 2020 09:17:13 +0800 Subject: [PATCH 2/2] add energy calculation in fix/dplr --- source/lmp/fix_dplr.cpp | 36 +++++++++++++++++++++++++++++++++++- source/lmp/fix_dplr.h | 4 ++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/source/lmp/fix_dplr.cpp b/source/lmp/fix_dplr.cpp index 1119383129..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), efield(3, 0.0) + :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; @@ -112,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; @@ -383,11 +388,17 @@ void FixDPLR::post_force(int vflag) 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]; @@ -495,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 094b0aa3e1..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; @@ -46,6 +48,8 @@ namespace LAMMPS_NS { vector dipole_recd; vector dfcorr_buff; vector efield; + vector efield_fsum, efield_fsum_all; + int efield_force_flag; void get_valid_pairs(vector >& pairs); }; }