diff --git a/doc/src/dynamical_matrix.rst b/doc/src/dynamical_matrix.rst index 623ab4b4f1..514eb94ce4 100644 --- a/doc/src/dynamical_matrix.rst +++ b/doc/src/dynamical_matrix.rst @@ -75,7 +75,7 @@ See the :doc:`Build package ` doc page for more info. Related commands """""""""""""""" -:doc:`fix phonon ` +:doc:`fix phonon `, :doc:`fix numdiff `, :doc:`compute hma ` uses an analytic formulation of the Hessian provided by a pair_style's Pair::single\_hessian() function, diff --git a/doc/src/fix_num_diff.rst b/doc/src/fix_num_diff.rst deleted file mode 100644 index 189b26e3f7..0000000000 --- a/doc/src/fix_num_diff.rst +++ /dev/null @@ -1,78 +0,0 @@ -.. index:: fix numdiff - -fix numdiff command -==================== - -Syntax -"""""" - - -.. parsed-literal:: - - fix ID group-ID numdiff Nevery Delta - -* ID, group-ID are documented in :doc:`fix ` command -* numdiff = style name of this fix command -* Nevery = calculate force by finite difference every this many timesteps -* Delta = finite difference displacement length (distance units) - - -Examples -"""""""" - - -.. parsed-literal:: - - fix 1 all numdiff 1 0.0001 - fix 1 all numdiff 10 1e-6 - fix 1 all numdiff 100 0.01 - -Description -""""""""""" - -Calculate forces through finite difference of energy versus position. -The resulting per-atom forces can be used by :doc:`dump custom `. - -The group specified with the command means only atoms within the group -have their averages computed. Results are set to 0.0 for atoms not in -the group. - - ----------- - - -The *Nevery*\ argument specifies on what timesteps the force will -be used calculated by finite difference. - -The *Delta*\ argument specifies the positional displacement each -atom will undergo. - - ----------- - - -**Restart, fix\_modify, output, run start/stop, minimize info:** - -No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options -are relevant to this fix. No global scalar or vector quantities are -stored by this fix for access by various :doc:`output commands `. - -This fix produces a per-atom array which can be accessed by -various :doc:`output commands `. . This fix produces -a per-atom array of the forces calculated by finite difference. The -per-atom values should only be accessed on timesteps that are multiples -of *Nfreq* since that is when the finite difference forces are calculated. - -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is invoked during :doc:`energy minimization `. - -Restrictions -"""""""""""" - none - -Related commands -"""""""""""""""" - -:doc:`dynamical_matrix `, - -**Default:** none diff --git a/doc/src/fix_numdiff.rst b/doc/src/fix_numdiff.rst new file mode 100644 index 0000000000..20d95cdf0a --- /dev/null +++ b/doc/src/fix_numdiff.rst @@ -0,0 +1,97 @@ +.. index:: fix numdiff + +fix numdiff command +==================== + +Syntax +"""""" + + +.. parsed-literal:: + + fix ID group-ID numdiff Nevery Delta + +* ID, group-ID are documented in :doc:`fix ` command +* numdiff = style name of this fix command +* Nevery = calculate force by finite difference every this many timesteps +* Delta = finite difference displacement length (distance units) + + +Examples +"""""""" + + +.. parsed-literal:: + + fix 1 all numdiff 1 0.0001 + fix 1 all numdiff 10 1e-6 + fix 1 all numdiff 100 0.01 + +Description +""""""""""" + +Calculate forces through finite difference calculations of energy +versus position. These forces can be compared to analytic forces +computed by pair styles, bond styles, etc. E.g. for debugging +purposes. + +The group specified with the command means only atoms within the group +have their averages computed. Results are set to 0.0 for atoms not in +the group. + +This fix performs a loop over all atoms (in the group). For each atom +and each component of force it adds *Delta* to the position, and +computes the new energy of the entire system. It then subtracts +*Delta* from the original position and again computes the new energy +of the system. It then restores the original position. That +component of force is calculated as the difference in energy divided +by two times *Delta*. + +.. note:: + +The cost of each energy evaluation is essentially the cost of an MD +timestep. This invoking this fix once has a cost of 2N timesteps, +where N is the total number of atoms in the system (assuming all atoms +are included in the group). So this fix can be very expensive to use +for large systems. + +---------- + + +The *Nevery*\ argument specifies on what timesteps the force will +be used calculated by finite difference. + +The *Delta*\ argument specifies the positional displacement each +atom will undergo. + + +---------- + + +**Restart, fix\_modify, output, run start/stop, minimize info:** + +No information about this fix is written to :doc:`binary restart files +`. None of the :doc:`fix_modify ` options are +relevant to this fix. + +This fix produces a per-atom array which can be accessed by various +:doc:`output commands `, which stores the components of +the force on each atom as calculated by finite difference. The +per-atom values can only be accessed on timesteps that are multiples +of *Nfreq* since that is when the finite difference forces are +calculated. + +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. This fix is invoked during :doc:`energy +minimization `. + +Restrictions +"""""""""""" + none + +Related commands +"""""""""""""""" + +:doc:`dynamical_matrix `, + +**Default:** none diff --git a/src/fix_num_diff.cpp b/src/fix_numdiff.cpp similarity index 52% rename from src/fix_num_diff.cpp rename to src/fix_numdiff.cpp index e713072c91..ce66ca49f8 100644 --- a/src/fix_num_diff.cpp +++ b/src/fix_numdiff.cpp @@ -9,94 +9,69 @@ the GNU General Public License. See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ - Fix_num_diff was created by Charles Sievers (UC Davis) +/* ---------------------------------------------------------------------- + Contributing author: Charles Sievers (UC Davis) ------------------------------------------------------------------------- */ #include "fix_num_diff.h" -#include #include -#include +#include #include "atom.h" +#include "domain.h" #include "update.h" #include "modify.h" -#include "comm.h" #include "compute.h" -#include "domain.h" -#include "region.h" #include "respa.h" +#include "force.h" #include "pair.h" #include "bond.h" #include "angle.h" #include "dihedral.h" #include "improper.h" #include "kspace.h" -#include "input.h" -#include "neighbor.h" -#include "timer.h" -#include "variable.h" #include "memory.h" #include "error.h" -#include "force.h" -#include "group.h" using namespace LAMMPS_NS; using namespace FixConst; -//enum{}; - /* ---------------------------------------------------------------------- */ FixNumDiff::FixNumDiff(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), numdiff_forces(NULL), temp_f(NULL), temp_x(NULL), id_pe(NULL) + Fix(lmp, narg, arg), numdiff_forces(NULL), temp_f(NULL), + temp_x(NULL), id_pe(NULL) { - respa_level_support = 1; - ilevel_respa = 0; - eflag = 1; - energy = 0.0; - size_peratom_cols = 3; - if (narg < 5) error->all(FLERR,"Illegal fix numdiff command"); - nevery = force->inumeric(FLERR,arg[3]); - del = force->numeric(FLERR,arg[4]); - - if (nevery <= 0) - error->all(FLERR,"Illegal fix numdiff command"); - peratom_flag = 1; - if (force->pair && force->pair->compute_flag) pair_compute_flag = 1; - else pair_compute_flag = 0; - if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1; - else kspace_compute_flag = 0; + peratom_freq = nevery; + size_peratom_cols = 3; + respa_level_support = 1; - numdiff_forces = NULL; - - maxatom1 = 0; + nevery = force->inumeric(FLERR,arg[3]); + delta = force->numeric(FLERR,arg[4]); + if (nevery <= 0 || delta <= 0.0) + error->all(FLERR,"Illegal fix numdiff command"); int n = strlen(id) + 6; id_pe = new char[n]; strcpy(id_pe,id); strcat(id_pe,"_pe"); - char **newarg = new char*[10]; + char **newarg = new char*[3]; newarg[0] = id_pe; newarg[1] = (char *) "all"; newarg[2] = (char *) "pe"; - newarg[3] = (char *) "pair"; - newarg[4] = (char *) "bond"; - newarg[5] = (char *) "angle"; - newarg[6] = (char *) "dihedral"; - newarg[7] = (char *) "improper"; - newarg[8] = (char *) "kspace"; - newarg[9] = (char *) "fix"; modify->add_compute(3,newarg); delete [] newarg; - nmax = 0; - memory->create(numdiff_forces,atom->natoms,3,"numdiff:numdiff_force"); - force_clear(numdiff_forces); - array_atom = numdiff_forces; + maxatom = 0; + numdiff_forces = NULL; + temp_x = NULL; + temp_f = NULL; + array_atom = NULL; } /* ---------------------------------------------------------------------- */ @@ -128,13 +103,22 @@ int FixNumDiff::setmask() void FixNumDiff::init() { - // check variables + // require consecutive atom IDs + + if (!atom->tag_enable || !atom->tag_consecutive()) + error->all(FLERR,"Fix numdiff requires consecutive atom IDs"); + + // check for PE compute int icompute = modify->find_compute(id_pe); - if (icompute < 0) - error->all(FLERR,"Compute ID for fix numdiff does not exist"); + if (icompute < 0) error->all(FLERR,"Compute ID for fix numdiff does not exist"); pe = modify->compute[icompute]; + if (force->pair && force->pair->compute_flag) pair_compute_flag = 1; + else pair_compute_flag = 0; + if (force->kspace && force->kspace->compute_flag) kspace_compute_flag = 1; + else kspace_compute_flag = 0; + if (strstr(update->integrate_style,"respa")) { ilevel_respa = ((Respa *) update->integrate)->nlevels-1; if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa); @@ -147,13 +131,7 @@ void FixNumDiff::post_force(int vflag) { if (update->ntimestep % nevery) return; - // energy and virial setup - calculate_forces(vflag); - - // if (lmp->kokkos) - // atom->sync_modify(Host, (unsigned int) (F_MASK | MASK_MASK), - // (unsigned int) F_MASK); } /* ---------------------------------------------------------------------- */ @@ -176,123 +154,130 @@ void FixNumDiff::min_post_force(int vflag) void FixNumDiff::calculate_forces(int vflag) { - int local_idx; // local index - bigint natoms = atom->natoms; - double **f = atom->f; - double **x = atom->x; - double temp = 1.0 / 2 / del; - int nlocal = atom->nlocal; - int *mask = atom->mask; - int flag = 0; - int allflag = 0; + int i,j,ilocal; + double energy; - if (atom->nmax > maxatom1) { + if (atom->nmax > maxatom) { memory->destroy(numdiff_forces); memory->destroy(temp_f); memory->destroy(temp_x); - maxatom1 = atom->nmax; - memory->create(numdiff_forces,maxatom1,3,"numdiff:numdiff_force"); - memory->create(temp_f,maxatom1,3,"numdiff:temp_f"); - memory->create(temp_x,maxatom1,3,"numdiff:temp_x"); + maxatom = atom->nmax; + memory->create(numdiff_forces,maxatom,3,"numdiff:numdiff_force"); + memory->create(temp_f,maxatom,3,"numdiff:temp_f"); + memory->create(temp_x,maxatom,3,"numdiff:temp_x"); array_atom = numdiff_forces; } - for (bigint i = 0; i < natoms; i++) - for (int j = 0; j < 3; j++) - temp_f[i][j] = f[i][j]; + // store copy of current coords and forces for owned and ghost atoms - for (bigint i = 0; i < natoms; i++) - for (int j = 0; j < 3; j++) + double **x = atom->x; + double **f = atom->f; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + for (i = 0; i < nall; i++) + for (j = 0; j < 3; j++) { temp_x[i][j] = x[i][j]; + temp_f[i][j] = f[i][j]; + } + + // initialize numerical forces to zero - //initialize forces to all zeros force_clear(numdiff_forces); - for (bigint i=1; i<=natoms; i++){ - local_idx = atom->map(i); - if (mask[local_idx] & groupbit) flag = 1; - else flag = 0; - MPI_Allreduce(&flag, &allflag,1, MPI_INT, MPI_SUM,world); - if (!allflag) continue; - for (int alpha=0; alpha<3; alpha++){ - displace_atom(local_idx, alpha, 1); - update_energy(vflag); - if (local_idx >= 0 && local_idx < nlocal) - numdiff_forces[local_idx][alpha] -= energy; + int flag,allflag; + double denominator = 0.5 / delta; - displace_atom(local_idx,alpha,-2); - update_energy(vflag); - if (local_idx >= 0 && local_idx < nlocal) { - numdiff_forces[local_idx][alpha] += energy; - numdiff_forces[local_idx][alpha] *= temp; + int *mask = atom->mask; + int ntotal = static_cast (atom->natoms); + int dimension = domain->dimension; + + for (tagint m = 1; m <= ntotal; m++) { + ilocal = atom->map(m); + flag = 0; + if ((ilocal >= 0 && ilocal < nlocal) && (mask[ilocal] & groupbit)) flag = 1; + MPI_Allreduce(&flag,&allflag,1,MPI_INT,MPI_SUM,world); + if (!allflag) continue; + + for (int idim = 0; idim < dimension; idim++) { + displace_atom(ilocal,idim,1); + energy = update_energy(vflag); + if (ilocal >= 0 && ilocal < nlocal) + numdiff_forces[ilocal][idim] -= energy; + + displace_atom(ilocal,idim,-2); + energy = update_energy(vflag); + if (ilocal >= 0 && ilocal < nlocal) { + numdiff_forces[ilocal][idim] += energy; + numdiff_forces[ilocal][idim] *= denominator; } - displace_atom(local_idx,alpha,1); + + // NOTE: will this introduce round-off in subsequent per-atom forces? + // maybe better to restore original coord here? + // in which case, replace displace_atom with something like + // reset_atom_coord(ilocal,idim,newcoord) + // in all 3 invocations in this method? + // then don't need to restore original coords at end, just forces + + displace_atom(ilocal,idim,1); } } - for (bigint i = 0; i < natoms; i++) - for (int j = 0; j < 3; j++) - f[i][j] = temp_f[i][j]; + // restore original coords and forces for owned and ghost atoms - for (bigint i = 0; i < natoms; i++) - for (int j = 0; j < 3; j++) + for (i = 0; i < nall; i++) + for (j = 0; j < 3; j++) { x[i][j] = temp_x[i][j]; - + f[i][j] = temp_f[i][j]; + } } /* ---------------------------------------------------------------------- - Displace atoms - ---------------------------------------------------------------------- */ + displace coord of all owned and ghost copies of ilocal +---------------------------------------------------------------------- */ -void FixNumDiff::displace_atom(int local_idx, int direction, int magnitude) +void FixNumDiff::displace_atom(int ilocal, int idim, int magnitude) { - if (local_idx < 0) return; + if (ilocal < 0) return; double **x = atom->x; int *sametag = atom->sametag; - int j = local_idx; - x[local_idx][direction] += del*magnitude; + int j = ilocal; + x[ilocal][idim] += delta*magnitude; - while (sametag[j] >= 0){ + while (sametag[j] >= 0) { j = sametag[j]; - x[j][direction] += del*magnitude; + x[j][idim] += delta*magnitude; } } /* ---------------------------------------------------------------------- evaluate potential energy and forces - may migrate atoms due to reneighboring - return new energy, which should include nextra_global dof - return negative gradient stored in atom->f - return negative gradient for nextra_global dof in fextra + same logic as in Verlet ------------------------------------------------------------------------- */ -void FixNumDiff::update_energy(int vflag) +double FixNumDiff::update_energy(int vflag) { force_clear(atom->f); - if (pair_compute_flag) { - force->pair->compute(eflag,vflag); - timer->stamp(Timer::PAIR); - } + int eflag = 1; + int vvflag = 0; // NOTE: + // I think we just want to disable virial comp here ? + // else we will be changing what the pair style stores ? + + if (pair_compute_flag) force->pair->compute(eflag,vvflag); + if (atom->molecular) { - if (force->bond) { - force->bond->compute(eflag,vflag); - } if (force->angle) { - force->angle->compute(eflag,vflag); - } if (force->dihedral) { - force->dihedral->compute(eflag,vflag); - } if (force->improper) { - force->improper->compute(eflag,vflag); - } - timer->stamp(Timer::BOND); - } - if (kspace_compute_flag) { - force->kspace->compute(eflag,vflag); - timer->stamp(Timer::KSPACE); + if (force->bond) force->bond->compute(eflag,vvflag); + if (force->angle) force->angle->compute(eflag,vvflag); + if (force->dihedral) force->dihedral->compute(eflag,vvflag); + if (force->improper) force->improper->compute(eflag,vvflag); } - energy = pe->compute_scalar(); + if (kspace_compute_flag) force->kspace->compute(eflag,vvflag); + + double energy = pe->compute_scalar(); + return energy; } /* ---------------------------------------------------------------------- @@ -301,57 +286,18 @@ void FixNumDiff::update_energy(int vflag) void FixNumDiff::force_clear(double **forces) { - size_t nbytes = sizeof(double) * atom->nlocal; if (force->newton) nbytes += sizeof(double) * atom->nghost; - - if (nbytes) { - memset(&forces[0][0],0,3*nbytes); - } + if (nbytes) memset(&forces[0][0],0,3*nbytes); } /* ---------------------------------------------------------------------- - pack values in local atom-based array for exchange with another proc -------------------------------------------------------------------------- */ - -int FixNumDiff::pack_exchange(int i, double *buf) -{ - int n = 0; - buf[n++] = numdiff_forces[i][0]; - buf[n++] = numdiff_forces[i][1]; - buf[n++] = numdiff_forces[i][2]; - return n; -} - -/* ---------------------------------------------------------------------- - unpack values in local atom-based array from exchange with another proc -------------------------------------------------------------------------- */ - -int FixNumDiff::unpack_exchange(int nlocal, double *buf) -{ - int n = 0; - numdiff_forces[nlocal][0] = buf[n++]; - numdiff_forces[nlocal][1] = buf[n++]; - numdiff_forces[nlocal][2] = buf[n++]; - return n; -} - -/* ---------------------------------------------------------------------- - return force calculated by numerical difference -------------------------------------------------------------------------- */ - -double FixNumDiff::compute_array(int i, int j) -{ - return numdiff_forces[i][j]; -} - -/* ---------------------------------------------------------------------- - memory usage of local atom-based array + memory usage of local atom-based arrays ------------------------------------------------------------------------- */ double FixNumDiff::memory_usage() { bigint bytes = 0.0; - bytes += 3 * atom->natoms * 3 * sizeof(double); // temp_f, temp_x, numdiff_f + bytes += 3 * maxatom*3 * sizeof(double); return bytes; } diff --git a/src/fix_num_diff.h b/src/fix_numdiff.h similarity index 69% rename from src/fix_num_diff.h rename to src/fix_numdiff.h index 29f7aa8681..1497139a01 100644 --- a/src/fix_num_diff.h +++ b/src/fix_numdiff.h @@ -9,8 +9,6 @@ the GNU General Public License. See the README file in the top-level LAMMPS directory. - - Fix_num_diff was created by Charles Sievers (UC Davis) ------------------------------------------------------------------------- */ #ifdef FIX_CLASS @@ -35,47 +33,28 @@ class FixNumDiff : public Fix { void post_force(int); void post_force_respa(int, int, int); void min_post_force(int); - double compute_array(int, int); double memory_usage(); - int pack_exchange(int i, double *buf); - int unpack_exchange(int nlocal, double *buf); -protected: - int eflag; // flags for energy/virial computation - int external_force_clear; // clear forces locally or externally - - int triclinic; // 0 if domain is orthog, 1 if triclinic - int pairflag; +private: + double delta; + int maxatom; + int ilevel_respa; int pair_compute_flag; // 0 if pair->compute is skipped int kspace_compute_flag; // 0 if kspace->compute is skipped - double **numdiff_forces; // local forces from numerical difference (this might be usefull for debugging) - - void update_energy(int vflag); - void force_clear(double **forces); - // virtual void openfile(const char* filename); - - private: - void create_groupmap(); - void displace_atom(int local_idx, int direction, int magnitude); - void calculate_forces(int vflag); - void compute_energy(); - - int ilevel_respa; - double del; - int nmax; - - int scaleflag; - int me; - double **temp_f; - double **temp_x; - double energy; - int maxatom1; - char *id_pe; class Compute *pe; + double **numdiff_forces; // finite diff forces + double **temp_f; // original forces + double **temp_x; // original coords + + double update_energy(int vflag); + void force_clear(double **forces); + void create_groupmap(); + void displace_atom(int local_idx, int direction, int magnitude); + void calculate_forces(int vflag); }; }