mods to fix numdiff

Related commands

""""""""""""""""

:doc:`fix phonon <fix_phonon>`

:doc:`fix phonon <fix_phonon>`, :doc:`fix numdiff <fix_numdiff>`, 

:doc:`compute hma <compute_hma>` uses an analytic formulation of the

Hessian provided by a pair_style's Pair::single\_hessian() function,

Description

"""""""""""

Calculate forces through finite difference of energy versus position.

The resulting per-atom forces can be used by :doc:`dump custom <dump>`.

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.

----------

**Restart, fix\_modify, output, run start/stop, minimize info:**

No information about this fix is written to :doc:`binary restart files <restart>`.  None of the :doc:`fix_modify <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 <Howto_output>`.

No information about this fix is written to :doc:`binary restart files

<restart>`.  None of the :doc:`fix_modify <fix_modify>` options are

relevant to this fix.

This fix produces a per-atom array which can be accessed by

various :doc:`output commands <Howto_output>`.  .  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.

This fix produces a per-atom array which can be accessed by various

:doc:`output commands <Howto_output>`, 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 <run>` command.  This fix is invoked during :doc:`energy minimization <minimize>`.

the :doc:`run <run>` command.  This fix is invoked during :doc:`energy

minimization <minimize>`.

Restrictions

""""""""""""

   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 <algorithm>

#include <mpi.h>

#include <cstring>

#include <memory.h>

#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;

  numdiff_forces = NULL;

  peratom_freq = nevery;

  size_peratom_cols = 3;

  respa_level_support = 1;

  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;

}

/* ---------------------------------------------------------------------- */

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);

{

  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);

}

/* ---------------------------------------------------------------------- */

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;

  int flag,allflag;

  double denominator = 0.5 / delta;

  int *mask = atom->mask;

  int ntotal = static_cast<tagint> (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 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;

      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;

    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) {

    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;

}

/* ----------------------------------------------------------------------

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);

  }

}

/* ----------------------------------------------------------------------

   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];

  if (nbytes) memset(&forces[0][0],0,3*nbytes);

}

/* ----------------------------------------------------------------------

   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;

}

   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

  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)

  char *id_pe;

  class Compute *pe;

  void update_energy(int vflag);

  void force_clear(double **forces);

  // virtual void openfile(const char* filename);

  double **numdiff_forces;          // finite diff forces

  double **temp_f;                  // original forces

  double **temp_x;                  // original coords

 private:

  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);

  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;

};

}