lammps-sys 0.6.0

/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing authors:
   Xiaohu Hu, CMB/ORNL (hux2@ornl.gov)
   David Hyde-Volpe, Tigran Abramyan, and Robert A. Latour (Clemson University)
   Chris Lorenz (Kings College-London)

   Implementation of the CHARMM CMAP; adds an extra energy term for the
   peptide backbone dihedrals.  The tools/ch2lmp/charmm2lammps.pl
   conversion script, which generates an extra section in the LAMMPS data
   file, is needed in order to generate the info used by this fix style.

   References:
   - MacKerell et al., J. Am. Chem. Soc. 126(2004):698-699.
   - MacKerell et al., J. Comput. Chem. 25(2004):1400-1415.
------------------------------------------------------------------------- */

#include "fix_cmap.h"
#include <mpi.h>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "force.h"
#include "comm.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"

using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;

#define MAXLINE 256
#define LISTDELTA 10000
#define LB_FACTOR 1.5

#define CMAPMAX 6   // max # of CMAP terms stored by one atom
#define CMAPDIM 24  // grid map dimension is 24 x 24
#define CMAPXMIN -360.0
#define CMAPXMIN2 -180.0
#define CMAPDX 15.0 // 360/CMAPDIM

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

FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg),
  crosstermlist(NULL), num_crossterm(NULL), crossterm_type(NULL),
  crossterm_atom1(NULL), crossterm_atom2(NULL), crossterm_atom3(NULL),
  crossterm_atom4(NULL), crossterm_atom5(NULL),
  g_axis(NULL), cmapgrid(NULL), d1cmapgrid(NULL), d2cmapgrid(NULL),
  d12cmapgrid(NULL)
{
  if (narg != 4) error->all(FLERR,"Illegal fix cmap command");

  restart_global = 1;
  restart_peratom = 1;
  peatom_flag = 1;
  virial_flag = 1;
  thermo_virial = 1;
  peratom_freq = 1;
  scalar_flag = 1;
  global_freq = 1;
  extscalar = 1;
  extvector = 1;
  wd_header = 1;
  wd_section = 1;

  MPI_Comm_rank(world,&me);
  MPI_Comm_size(world,&nprocs);

  // allocate memory for CMAP data

  memory->create(g_axis,CMAPDIM,"cmap:g_axis");
  memory->create(cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:grid");
  memory->create(d1cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d1grid");
  memory->create(d2cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d2grid");
  memory->create(d12cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d12grid");

  // read and setup CMAP data

  read_grid_map(arg[3]);

  // perform initial allocation of atom-based arrays
  // register with Atom class

  num_crossterm = NULL;
  crossterm_type = NULL;
  crossterm_atom1 = NULL;
  crossterm_atom2 = NULL;
  crossterm_atom3 = NULL;
  crossterm_atom4 = NULL;
  crossterm_atom5 = NULL;

  nmax_previous = 0;
  grow_arrays(atom->nmax);
  atom->add_callback(0);
  atom->add_callback(1);

  // local list of crossterms

  ncmap = 0;
  maxcrossterm = 0;
  crosstermlist = NULL;
}

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

FixCMAP::~FixCMAP()
{
  // unregister callbacks to this fix from Atom class

  atom->delete_callback(id,0);
  atom->delete_callback(id,1);

  memory->destroy(g_axis);
  memory->destroy(cmapgrid);
  memory->destroy(d1cmapgrid);
  memory->destroy(d2cmapgrid);
  memory->destroy(d12cmapgrid);

  memory->destroy(crosstermlist);

  memory->destroy(num_crossterm);
  memory->destroy(crossterm_type);
  memory->destroy(crossterm_atom1);
  memory->destroy(crossterm_atom2);
  memory->destroy(crossterm_atom3);
  memory->destroy(crossterm_atom4);
  memory->destroy(crossterm_atom5);
}

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

int FixCMAP::setmask()
{
  int mask = 0;
  mask |= PRE_NEIGHBOR;
  mask |= PRE_REVERSE;
  mask |= POST_FORCE;
  mask |= THERMO_ENERGY;
  mask |= POST_FORCE_RESPA;
  mask |= MIN_POST_FORCE;
  return mask;
}

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

void FixCMAP::init()
{
  int i;
  double angle;

  i = 0;
  angle = -180.0;
  while (angle < 180.0) {
    g_axis[i] = angle;
    angle += CMAPDX;
    i++;
  }

  // pre-compute the derivatives of the maps

  for (i = 0; i < 6; i++)
    set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]);

  // define newton_bond here in case restart file was read (not data file)

  newton_bond = force->newton_bond;
}

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

void FixCMAP::setup(int vflag)
{
  pre_neighbor();

  if (strstr(update->integrate_style,"verlet"))
    post_force(vflag);
  else {
    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
    post_force_respa(vflag,nlevels_respa-1,0);
    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
  }
}

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

void FixCMAP::setup_pre_neighbor()
{
  pre_neighbor();
}

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

void FixCMAP::setup_pre_reverse(int eflag, int vflag)
{
  pre_reverse(eflag,vflag);
}

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

void FixCMAP::min_setup(int vflag)
{
  pre_neighbor();
  post_force(vflag);
}

/* ----------------------------------------------------------------------
   store local neighbor list as if newton_bond = OFF, even if actually ON
------------------------------------------------------------------------- */

void FixCMAP::pre_neighbor()
{
  int i,m,atom1,atom2,atom3,atom4,atom5;

  // guesstimate initial length of local crossterm list
  // if ncmap was not set (due to read_restart, no read_data),
  //   then list will grow by LISTDELTA chunks

  if (maxcrossterm == 0) {
    if (nprocs == 1) maxcrossterm = ncmap;
    else maxcrossterm = static_cast<int> (LB_FACTOR*ncmap/nprocs);
    memory->create(crosstermlist,maxcrossterm,6,"cmap:crosstermlist");
  }

  int nlocal = atom->nlocal;

  ncrosstermlist = 0;

  for (i = 0; i < nlocal; i++) {
    for (m = 0; m < num_crossterm[i]; m++) {
      atom1 = atom->map(crossterm_atom1[i][m]);
      atom2 = atom->map(crossterm_atom2[i][m]);
      atom3 = atom->map(crossterm_atom3[i][m]);
      atom4 = atom->map(crossterm_atom4[i][m]);
      atom5 = atom->map(crossterm_atom5[i][m]);

      if (atom1 == -1 || atom2 == -1 || atom3 == -1 ||
          atom4 == -1 || atom5 == -1) {
        char str[128];
        sprintf(str,"CMAP atoms "
                TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT " "
                TAGINT_FORMAT " " TAGINT_FORMAT
                " missing on proc %d at step " BIGINT_FORMAT,
                crossterm_atom1[i][m],crossterm_atom2[i][m],
                crossterm_atom3[i][m],crossterm_atom4[i][m],
                crossterm_atom5[i][m],me,update->ntimestep);
        error->one(FLERR,str);
      }
      atom1 = domain->closest_image(i,atom1);
      atom2 = domain->closest_image(i,atom2);
      atom3 = domain->closest_image(i,atom3);
      atom4 = domain->closest_image(i,atom4);
      atom5 = domain->closest_image(i,atom5);

      if (i <= atom1 && i <= atom2 && i <= atom3 &&
          i <= atom4 && i <= atom5) {
        if (ncrosstermlist == maxcrossterm) {
          maxcrossterm += LISTDELTA;
          memory->grow(crosstermlist,maxcrossterm,6,"cmap:crosstermlist");
        }
        crosstermlist[ncrosstermlist][0] = atom1;
        crosstermlist[ncrosstermlist][1] = atom2;
        crosstermlist[ncrosstermlist][2] = atom3;
        crosstermlist[ncrosstermlist][3] = atom4;
        crosstermlist[ncrosstermlist][4] = atom5;
        crosstermlist[ncrosstermlist][5] = crossterm_type[i][m];
        ncrosstermlist++;
      }
    }
  }
}

/* ----------------------------------------------------------------------
   store eflag, so can use it in post_force to tally per-atom energies
------------------------------------------------------------------------- */

void FixCMAP::pre_reverse(int eflag, int /*vflag*/)
{
  eflag_caller = eflag;
}

/* ----------------------------------------------------------------------
   compute CMAP terms as if newton_bond = OFF, even if actually ON
------------------------------------------------------------------------- */

void FixCMAP::post_force(int vflag)
{
  int n,i1,i2,i3,i4,i5,type,nlist;
  int li1, li2, mli1,mli2,mli11,mli21,t1,li3,li4,mli3,mli4,mli31,mli41;
  int list[5];
  // vectors needed to calculate the cross-term dihedral angles
  double vb21x,vb21y,vb21z,vb32x,vb32y,vb32z,vb34x,vb34y,vb34z;
  double vb23x,vb23y,vb23z;
  double vb43x,vb43y,vb43z,vb45x,vb45y,vb45z,a1x,a1y,a1z,b1x,b1y,b1z;
  double a2x,a2y,a2z,b2x,b2y,b2z,r32,a1sq,b1sq,a2sq,b2sq,dpr21r32,dpr34r32;
  double dpr32r43,dpr45r43,r43,vb12x,vb12y,vb12z,vb54x,vb54y,vb54z;
  // cross-term dihedral angles
  double phi,psi,phi1,psi1;
  double f1[3],f2[3],f3[3],f4[3],f5[3],vcmap[6];
  double gs[4],d1gs[4],d2gs[4],d12gs[4];
  double engfraction;
  // vectors needed for the gradient/force calculation
  double dphidr1x,dphidr1y,dphidr1z,dphidr2x,dphidr2y,dphidr2z;
  double dphidr3x,dphidr3y,dphidr3z,dphidr4x,dphidr4y,dphidr4z;
  double dpsidr1x,dpsidr1y,dpsidr1z,dpsidr2x,dpsidr2y,dpsidr2z;
  double dpsidr3x,dpsidr3y,dpsidr3z,dpsidr4x,dpsidr4y,dpsidr4z;

  // Definition of cross-term dihedrals

  //         phi dihedral
  //   |--------------------|
  //   a1-----a2-----a3-----a4-----a5    cross-term atoms
  //   C      N      CA     C      N     cross-term atom types
  //          |--------------------|
  //               psi dihedral

  double **x = atom->x;
  double **f = atom->f;
  int nlocal = atom->nlocal;

  ecmap = 0.0;
  int eflag = eflag_caller;
  ev_init(eflag,vflag);

  for (n = 0; n < ncrosstermlist; n++) {
    i1 = crosstermlist[n][0];
    i2 = crosstermlist[n][1];
    i3 = crosstermlist[n][2];
    i4 = crosstermlist[n][3];
    i5 = crosstermlist[n][4];

    type = crosstermlist[n][5];
    if (type == 0) continue;

    // calculate bond vectors for both dihedrals

    // phi
    // vb21 = r2 - r1

      vb21x = x[i2][0] - x[i1][0];
      vb21y = x[i2][1] - x[i1][1];
      vb21z = x[i2][2] - x[i1][2];
      vb12x = -1.0*vb21x;
      vb12y = -1.0*vb21y;
      vb12z = -1.0*vb21z;
      vb32x = x[i3][0] - x[i2][0];
      vb32y = x[i3][1] - x[i2][1];
      vb32z = x[i3][2] - x[i2][2];
      vb23x = -1.0*vb32x;
      vb23y = -1.0*vb32y;
      vb23z = -1.0*vb32z;

      vb34x = x[i3][0] - x[i4][0];
      vb34y = x[i3][1] - x[i4][1];
      vb34z = x[i3][2] - x[i4][2];

      // psi
      // bond vectors same as for phi: vb32

      vb43x = -1.0*vb34x;
      vb43y = -1.0*vb34y;
      vb43z = -1.0*vb34z;

      vb45x = x[i4][0] - x[i5][0];
      vb45y = x[i4][1] - x[i5][1];
      vb45z = x[i4][2] - x[i5][2];
      vb54x = -1.0*vb45x;
      vb54y = -1.0*vb45y;
      vb54z = -1.0*vb45z;

      // calculate normal vectors for planes that define the dihedral angles

      a1x = vb12y*vb23z - vb12z*vb23y;
      a1y = vb12z*vb23x - vb12x*vb23z;
      a1z = vb12x*vb23y - vb12y*vb23x;

      b1x = vb43y*vb23z - vb43z*vb23y;
      b1y = vb43z*vb23x - vb43x*vb23z;
      b1z = vb43x*vb23y - vb43y*vb23x;

      a2x = vb23y*vb34z - vb23z*vb34y;
      a2y = vb23z*vb34x - vb23x*vb34z;
      a2z = vb23x*vb34y - vb23y*vb34x;

      b2x = vb45y*vb43z - vb45z*vb43y;
      b2y = vb45z*vb43x - vb45x*vb43z;
      b2z = vb45x*vb43y - vb45y*vb43x;

      // calculate terms used later in calculations

      r32 = sqrt(vb32x*vb32x + vb32y*vb32y + vb32z*vb32z);
      a1sq = a1x*a1x + a1y*a1y + a1z*a1z;
      b1sq = b1x*b1x + b1y*b1y + b1z*b1z;

      r43 = sqrt(vb43x*vb43x + vb43y*vb43y + vb43z*vb43z);
      a2sq = a2x*a2x + a2y*a2y + a2z*a2z;
      b2sq = b2x*b2x + b2y*b2y + b2z*b2z;
      //if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001)
      //  printf("a1sq b1sq a2sq b2sq: %f %f %f %f \n",a1sq,b1sq,a2sq,b2sq);
      if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) continue;
      dpr21r32 = vb21x*vb32x + vb21y*vb32y + vb21z*vb32z;
      dpr34r32 = vb34x*vb32x + vb34y*vb32y + vb34z*vb32z;
      dpr32r43 = vb32x*vb43x + vb32y*vb43y + vb32z*vb43z;
      dpr45r43 = vb45x*vb43x + vb45y*vb43y + vb45z*vb43z;

      // calculate the backbone dihedral angles as VMD and GROMACS

      phi = dihedral_angle_atan2(vb21x,vb21y,vb21z,a1x,a1y,a1z,b1x,b1y,b1z,r32);
      psi = dihedral_angle_atan2(vb32x,vb32y,vb32z,a2x,a2y,a2z,b2x,b2y,b2z,r43);

      if (phi == 180.0) phi= -180.0;
      if (psi == 180.0) psi= -180.0;

      phi1 = phi;
      if (phi1 < 0.0) phi1 += 360.0;
      psi1 = psi;
      if (psi1 < 0.0) psi1 += 360.0;

      // find the neighbor grid point index

      li1 = int(((phi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));
      li2 = int(((psi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));

      li3 = int((phi-CMAPXMIN2)/CMAPDX);
      li4 = int((psi-CMAPXMIN2)/CMAPDX);
      mli3 = li3 % CMAPDIM;
      mli4 = li4 % CMAPDIM;
      mli31 = (li3+1) % CMAPDIM;
      mli41 = (li4+1)  %CMAPDIM;
      mli1 = li1 % CMAPDIM;
      mli2 = li2 % CMAPDIM;
      mli11 = (li1+1) % CMAPDIM;
      mli21 = (li2+1)  %CMAPDIM;
      t1 = type-1;
      if (t1 < 0 || t1 > 5) error->all(FLERR,"Invalid CMAP crossterm_type");

      // determine the values and derivatives for the grid square points

      gs[0] = cmapgrid[t1][mli3][mli4];
      gs[1] = cmapgrid[t1][mli31][mli4];
      gs[2] = cmapgrid[t1][mli31][mli41];
      gs[3] = cmapgrid[t1][mli3][mli41];
      d1gs[0] = d1cmapgrid[t1][mli1][mli2];
      d1gs[1] = d1cmapgrid[t1][mli11][mli2];
      d1gs[2] = d1cmapgrid[t1][mli11][mli21];
      d1gs[3] = d1cmapgrid[t1][mli1][mli21];

      d2gs[0] = d2cmapgrid[t1][mli1][mli2];
      d2gs[1] = d2cmapgrid[t1][mli11][mli2];
      d2gs[2] = d2cmapgrid[t1][mli11][mli21];
      d2gs[3] = d2cmapgrid[t1][mli1][mli21];

      d12gs[0] = d12cmapgrid[t1][mli1][mli2];
      d12gs[1] = d12cmapgrid[t1][mli11][mli2];
      d12gs[2] = d12cmapgrid[t1][mli11][mli21];
      d12gs[3] = d12cmapgrid[t1][mli1][mli21];

      // calculate the cmap energy and the gradient (dE/dphi,dE/dpsi)

      bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs);

      // sum up cmap energy contributions

      engfraction = 0.2 * E;
      if (i1 < nlocal) ecmap += engfraction;
      if (i2 < nlocal) ecmap += engfraction;
      if (i3 < nlocal) ecmap += engfraction;
      if (i4 < nlocal) ecmap += engfraction;
      if (i5 < nlocal) ecmap += engfraction;

      // calculate the derivatives dphi/dr_i

      dphidr1x = 1.0*r32/a1sq*a1x;
      dphidr1y = 1.0*r32/a1sq*a1y;
      dphidr1z = 1.0*r32/a1sq*a1z;

      dphidr2x = -1.0*r32/a1sq*a1x - dpr21r32/a1sq/r32*a1x +
        dpr34r32/b1sq/r32*b1x;
      dphidr2y = -1.0*r32/a1sq*a1y - dpr21r32/a1sq/r32*a1y +
        dpr34r32/b1sq/r32*b1y;
      dphidr2z = -1.0*r32/a1sq*a1z - dpr21r32/a1sq/r32*a1z +
        dpr34r32/b1sq/r32*b1z;

      dphidr3x = dpr34r32/b1sq/r32*b1x - dpr21r32/a1sq/r32*a1x - r32/b1sq*b1x;
      dphidr3y = dpr34r32/b1sq/r32*b1y - dpr21r32/a1sq/r32*a1y - r32/b1sq*b1y;
      dphidr3z = dpr34r32/b1sq/r32*b1z - dpr21r32/a1sq/r32*a1z - r32/b1sq*b1z;

      dphidr4x = r32/b1sq*b1x;
      dphidr4y = r32/b1sq*b1y;
      dphidr4z = r32/b1sq*b1z;

      // calculate the derivatives dpsi/dr_i

      dpsidr1x = 1.0*r43/a2sq*a2x;
      dpsidr1y = 1.0*r43/a2sq*a2y;
      dpsidr1z = 1.0*r43/a2sq*a2z;

      dpsidr2x = r43/a2sq*a2x + dpr32r43/a2sq/r43*a2x - dpr45r43/b2sq/r43*b2x;
      dpsidr2y = r43/a2sq*a2y + dpr32r43/a2sq/r43*a2y - dpr45r43/b2sq/r43*b2y;
      dpsidr2z = r43/a2sq*a2z + dpr32r43/a2sq/r43*a2z - dpr45r43/b2sq/r43*b2z;

      dpsidr3x = dpr45r43/b2sq/r43*b2x - dpr32r43/a2sq/r43*a2x - r43/b2sq*b2x;
      dpsidr3y = dpr45r43/b2sq/r43*b2y - dpr32r43/a2sq/r43*a2y - r43/b2sq*b2y;
      dpsidr3z = dpr45r43/b2sq/r43*b2z - dpr32r43/a2sq/r43*a2z - r43/b2sq*b2z;

      dpsidr4x = r43/b2sq*b2x;
      dpsidr4y = r43/b2sq*b2y;
      dpsidr4z = r43/b2sq*b2z;

      // calculate forces on cross-term atoms: F = -(dE/dPhi)*(dPhi/dr)

      f1[0] = dEdPhi*dphidr1x;
      f1[1] = dEdPhi*dphidr1y;
      f1[2] = dEdPhi*dphidr1z;
      f2[0] = dEdPhi*dphidr2x + dEdPsi*dpsidr1x;
      f2[1] = dEdPhi*dphidr2y + dEdPsi*dpsidr1y;
      f2[2] = dEdPhi*dphidr2z + dEdPsi*dpsidr1z;
      f3[0] = -dEdPhi*dphidr3x - dEdPsi*dpsidr2x;
      f3[1] = -dEdPhi*dphidr3y - dEdPsi*dpsidr2y;
      f3[2] = -dEdPhi*dphidr3z - dEdPsi*dpsidr2z;
      f4[0] = -dEdPhi*dphidr4x - dEdPsi*dpsidr3x;
      f4[1] = -dEdPhi*dphidr4y - dEdPsi*dpsidr3y;
      f4[2] = -dEdPhi*dphidr4z - dEdPsi*dpsidr3z;
      f5[0] = -dEdPsi*dpsidr4x;
      f5[1] = -dEdPsi*dpsidr4y;
      f5[2] = -dEdPsi*dpsidr4z;

      // apply force to each of the 5 atoms

      if (i1 < nlocal) {
        f[i1][0] += f1[0];
        f[i1][1] += f1[1];
        f[i1][2] += f1[2];
      }
      if (i2 < nlocal) {
        f[i2][0] += f2[0];
        f[i2][1] += f2[1];
        f[i2][2] += f2[2];
      }
      if (i3 < nlocal) {
        f[i3][0] += f3[0];
        f[i3][1] += f3[1];
        f[i3][2] += f3[2];
      }
      if (i4 < nlocal) {
        f[i4][0] += f4[0];
        f[i4][1] += f4[1];
        f[i4][2] += f4[2];
      }
      if (i5 < nlocal) {
        f[i5][0] += f5[0];
        f[i5][1] += f5[1];
        f[i5][2] += f5[2];
      }

      // tally energy and/or virial

      if (evflag) {
        nlist = 0;
        if (i1 < nlocal) list[nlist++] = i1;
        if (i2 < nlocal) list[nlist++] = i2;
        if (i3 < nlocal) list[nlist++] = i3;
        if (i4 < nlocal) list[nlist++] = i4;
        if (i5 < nlocal) list[nlist++] = i5;
        vcmap[0] = (vb12x*f1[0])+(vb32x*f3[0])+((vb43x+vb32x)*f4[0])+
          ((vb54x+vb43x+vb32x)*f5[0]);
        vcmap[1] = (vb12y*f1[1])+(vb32y*f3[1])+((vb43y+vb32y)*f4[1])+
          ((vb54y+vb43y+vb32y)*f5[1]);
        vcmap[2] = (vb12z*f1[2])+(vb32z*f3[2])+((vb43z+vb32z)*f4[2])+
          ((vb54z+vb43z+vb32z)*f5[2]);
        vcmap[3] = (vb12x*f1[1])+(vb32x*f3[1])+((vb43x+vb32x)*f4[1])+
          ((vb54x+vb43x+vb32x)*f5[1]);
        vcmap[4] = (vb12x*f1[2])+(vb32x*f3[2])+((vb43x+vb32x)*f4[2])+
          ((vb54x+vb43x+vb32x)*f5[2]);
        vcmap[5] = (vb12y*f1[2])+(vb32y*f3[2])+((vb43y+vb32y)*f4[2])+
          ((vb54y+vb43y+vb32y)*f5[2]);
        ev_tally(nlist,list,5.0,E,vcmap);
        //ev_tally(5,list,nlocal,newton_bond,E,vcmap);
      }
  }
}

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

void FixCMAP::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
  if (ilevel == nlevels_respa-1) post_force(vflag);
}

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

void FixCMAP::min_post_force(int vflag)
{
  post_force(vflag);
}

/* ----------------------------------------------------------------------
   energy of CMAP term
------------------------------------------------------------------------- */

double FixCMAP::compute_scalar()
{
  double all;
  MPI_Allreduce(&ecmap,&all,1,MPI_DOUBLE,MPI_SUM,world);
  return all;
}

// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
// methods to read CMAP potential file, perform interpolation
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------

void FixCMAP::read_grid_map(char *cmapfile)
{
  char linebuf[MAXLINE];
  char *chunk,*line;
  int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter;

  FILE *fp = NULL;
  if (comm->me == 0) {
    fp = force->open_potential(cmapfile);
    if (fp == NULL) {
      char str[128];
      snprintf(str,128,"Cannot open fix cmap file %s",cmapfile);
      error->one(FLERR,str);
    }
  }

  for (int ix1 = 0; ix1 < 6; ix1++)
    for (int ix2 = 0; ix2 < CMAPDIM; ix2++)
      for (int ix3 = 0; ix3 < CMAPDIM; ix3++)
        cmapgrid[ix1][ix2][ix3] = 0.0;

  counter = 0;
  i1 = i2 = i3 = i4 = i5 = i6 = 0;
  j1 = j2 = j3 = j4 = j5 = j6 = 0;

  int done = 0;

  while (!done) {
    // only read on rank 0 and broadcast to all other ranks
    if (comm->me == 0)
      done = (fgets(linebuf,MAXLINE,fp) == NULL);

    MPI_Bcast(&done,1,MPI_INT,0,world);
    if (done) continue;

    MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world);

    // remove leading whitespace
    line = linebuf;
    while (line && (*line == ' ' || *line == '\t' || *line == '\r')) ++line;

    // skip if empty line or comment
    if (!line || *line =='\n' || *line == '\0' || *line == '#') continue;

    // read in the cmap grid point values
    // NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors
    //       will occur if content of the file "cmap.data" is altered
    //
    // Reading order of the maps:
    // 1. Alanine map
    // 2. Alanine before proline map
    // 3. Proline map
    // 4. Two adjacent prolines map
    // 5. Glycine map
    // 6. Glycine before proline map

    chunk = strtok(line, " \r\n");
    while (chunk != NULL) {

      // alanine map

      if (counter < CMAPDIM*CMAPDIM) {
        cmapgrid[0][i1][j1] = atof(chunk);
        chunk = strtok(NULL, " \r\n");
        j1++;
        if (j1 == CMAPDIM) {
          j1 = 0;
          i1++;
        }
        counter++;
      }

      // alanine-proline map

      else if (counter >= CMAPDIM*CMAPDIM &&
               counter < 2*CMAPDIM*CMAPDIM) {
        cmapgrid[1][i2][j2]= atof(chunk);
        chunk = strtok(NULL, " \r\n");
        j2++;
        if (j2 == CMAPDIM) {
          j2 = 0;
          i2++;
        }
        counter++;
      }

      // proline map

      else if (counter >= 2*CMAPDIM*CMAPDIM &&
               counter < 3*CMAPDIM*CMAPDIM) {
        cmapgrid[2][i3][j3] = atof(chunk);
        chunk = strtok(NULL, " \r\n");
        j3++;
        if (j3 == CMAPDIM) {
          j3 = 0;
          i3++;
        }
        counter++;
      }

      // 2 adjacent prolines map

      else if (counter >= 3*CMAPDIM*CMAPDIM &&
               counter < 4*CMAPDIM*CMAPDIM) {
        cmapgrid[3][i4][j4] = atof(chunk);
        chunk = strtok(NULL, " \r\n");
        j4++;
        if (j4 == CMAPDIM) {
          j4 = 0;
          i4++;
        }
        counter++;
      }

      // glycine map

      else if (counter >= 4*CMAPDIM*CMAPDIM &&
               counter < 5*CMAPDIM*CMAPDIM) {
        cmapgrid[4][i5][j5] = atof(chunk);
        chunk = strtok(NULL, " \r\n");
        j5++;
        if (j5 == CMAPDIM) {
          j5 = 0;
          i5++;
        }
        counter++;
      }

      // glycine-proline map

      else if (counter >= 5*CMAPDIM*CMAPDIM &&
               counter < 6*CMAPDIM*CMAPDIM) {
        cmapgrid[5][i6][j6] = atof(chunk);
        chunk = strtok(NULL, " \r\n");
        j6++;
        if (j6 == CMAPDIM) {
          j6 = 0;
          i6++;
        }
        counter++;
      }

      else break;
    }
  }

  if (comm->me == 0) fclose(fp);
}

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

void FixCMAP::spline(double *y, double *ddy, int n)
{
  // create the 2nd dervatives of a taublated function y_i(x_i)
  // at the tabulated points

  int i, j;
  double p, *u;

  memory->create(u,n-1,"cmap:u");

  ddy[0] = u[0] = 0.0;

  for (i = 1; i <= n-2; i++) {
    p = 1.0/(ddy[i-1]+4.0);
    ddy[i] = -p;
    u[i] = ((((6.0*y[i+1])-(12.0*y[i])+(6.0*y[i-1]))/(CMAPDX*CMAPDX))-u[i-1])*p;
  }

  ddy[n-1] = 0.0;

  for (j = n-2; j >= 0; j--)
    ddy[j] = ddy[j]*ddy[j+1] + u[j];

  memory->destroy(u);
}

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

void FixCMAP::spl_interpolate(double x, double *y, double *ddy, double &yo,
                              double &dyo)
{
  // perform a 1D cubic spline interpolation

  int ix;
  double a,b,a1,b1,a2,b2;

  ix = int((x-CMAPXMIN)/CMAPDX-(1./2.));

  a = (CMAPXMIN+(ix*1.0)*CMAPDX-x)/CMAPDX;
  b = (x-CMAPXMIN-(((ix-1)*1.0)*CMAPDX))/CMAPDX;

  a1 = a*a*a-a;
  b1 = b*b*b-b;

  a2 = 3.0*a*a-1.0;
  b2 = 3.0*b*b-1.0;
  yo = a*y[ix]+b*y[ix+1]+(a1*ddy[ix]+b1*ddy[ix+1])*(CMAPDX*CMAPDX)/6.0;
  dyo = (y[ix+1]-y[ix])/CMAPDX-a2/6.0*CMAPDX*ddy[ix]+b2/6.0*CMAPDX*ddy[ix+1];
}

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

void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
                                  double **d12yo)
{
  // precompute the gradient and cross-derivatives of the map grid points.
  // use the bicubic spline to calculate the derivatives

  int i, j, k, ii, jj, xm, p;
  double phi, psi, d1y, d2y, d12y, tyyk,tdyk;
  double *tmp_y, *tmp_dy, *tmp_ddy, **tmap, **tddmap;
  int ix;
  double a,b,a1,b1,a2,b2;

  xm = CMAPDIM/2;
  p = CMAPDIM;

  d1y = 0.;
  d2y = 0.;
  d12y = 0.;

  memory->create(tmp_y,CMAPDIM*2,"cmap:tmp_y");
  memory->create(tmp_dy,CMAPDIM*2,"cmap:tmp_dy");
  memory->create(tmp_ddy,CMAPDIM*2,"cmap:tmp_ddy");
  memory->create(tmap,CMAPDIM*2,CMAPDIM*2,"cmap:tmap");
  memory->create(tddmap,CMAPDIM*2,CMAPDIM*2,"cmap:tddmap");

  // periodically expand the original map
  // use the expanded map for bicubic spline interpolation,
  //   which is used to obtain the derivatives
  // actual interpolation is done with bicubic interpolation

  for (i = 0; i < CMAPDIM*2; i++) {
    ii = ((i+CMAPDIM-xm)%CMAPDIM);
    for (j = 0; j < CMAPDIM*2; j++) {
      jj = ((j+CMAPDIM-xm)%CMAPDIM);
      tmap[i][j] = map[ii][jj];
    }
  }

  for (i = 0; i < CMAPDIM*2; i++)
    spline(tmap[i], tddmap[i], CMAPDIM*2);

  for (i = xm; i < CMAPDIM+xm; i++) {
    phi = (i-xm)*CMAPDX-180.0;
    for (j = xm; j < CMAPDIM+xm; j++) {
      psi = (j-xm)*CMAPDX-180.0;
      ix = int((psi-CMAPXMIN)/CMAPDX);
      a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-psi)/CMAPDX;
      b = (psi-CMAPXMIN-((ix)*1.0)*CMAPDX)/CMAPDX;
      a1 = a*a*a-a;
      b1 = b*b*b-b;
      a2 = 3.0*a*a-1.0;
      b2 = 3.0*b*b-1.0;
      for (k = 0; k < CMAPDIM*2; k++) {
        tyyk = tmp_y[k];
        tdyk = tmp_dy[k];
        tyyk = a*tmap[k][ix]+b*tmap[k][ix+1]+
          (a1*tddmap[k][ix]+b1*tddmap[k][ix+1])*(CMAPDX*CMAPDX)/6.0;
        tdyk = (tmap[k][ix+1]-tmap[k][ix])/CMAPDX-
          (a2/6.0*CMAPDX*tddmap[k][ix])+(b2/6.0*CMAPDX*tddmap[k][ix+1]);
        tmp_y[k] = tyyk;
        tmp_dy[k] = tdyk;
      }

      spline(tmp_y,tmp_ddy,CMAPDIM+xm+xm);
      ix = int((phi-CMAPXMIN)/CMAPDX);
      a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-phi)/CMAPDX;
      b = (phi-CMAPXMIN-(ix*1.0)*CMAPDX)/CMAPDX;
      a1 = a*a*a-a;
      b1 = b*b*b-b;
      a2 = 3.0*a*a-1.0;
      b2 = 3.0*b*b-1.0;
      d1y = (tmp_y[ix+1]-tmp_y[ix])/CMAPDX-
        a2/6.0*CMAPDX*tmp_ddy[ix]+b2/6.0*CMAPDX*tmp_ddy[ix+1];
      spline(tmp_dy,tmp_ddy,CMAPDIM+xm+xm);
      ix = int((phi-CMAPXMIN)/CMAPDX);
      a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-phi)/CMAPDX;
      b = (phi-CMAPXMIN-(ix*1.0)*CMAPDX)/CMAPDX;
      a1 = a*a*a-a;
      b1 = b*b*b-b;
      a2 = 3.0*a*a-1.0;
      b2 = 3.0*b*b-1.0;
      d2y = a*tmp_dy[ix]+b*tmp_dy[ix+1]+
        (a1*tmp_ddy[ix]+b1*tmp_ddy[ix+1])*(CMAPDX*CMAPDX)/6.0;
      d12y = (tmp_dy[ix+1]-tmp_dy[ix])/CMAPDX-
        a2/6.0*CMAPDX*tmp_ddy[ix]+b2/6.0*CMAPDX*tmp_ddy[ix+1];
      d1yo[i%p][j%p] = d1y;
      d2yo[i%p][j%p] = d2y;
      d12yo[i%p][j%p] = d12y;
    }
  }

  memory->destroy(tmp_y);
  memory->destroy(tmp_dy);
  memory->destroy(tmp_ddy);
  memory->destroy(tmap);
  memory->destroy(tddmap);
}

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

double FixCMAP::dihedral_angle_atan2(double fx, double fy, double fz,
                                      double ax, double ay, double az,
                                      double bx, double by, double bz,
                                      double absg)
{
  // calculate the dihedral angle

  double angle = 0.0, arg1, arg2;

  arg1 = absg*(fx*bx+fy*by+fz*bz);
  arg2 = ax*bx+ay*by+az*bz;

  if (arg1 == 0 && arg2 == 0)
    error->all(FLERR,"CMAP: atan2 function cannot take 2 zero arguments");
  else {
    angle = atan2(arg1,arg2);
    angle = angle*180.0/MY_PI;
  }

  return angle;
}

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

void FixCMAP::bc_coeff(double *gs, double *d1gs, double *d2gs, double *d12gs)
{
  // calculate the bicubic interpolation coefficients c_ij

  static int wt[16][16] =
    { {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
      {-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1, 0, 0, 0, 0},
      {2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0},
      {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0},
      {0, 0, 0, 0,-3, 0, 0, 3, 0, 0, 0, 0,-2, 0, 0,-1},
      {0, 0, 0, 0, 2, 0, 0,-2, 0, 0, 0, 0, 1, 0, 0, 1},
      {-3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0},
      {9,-9, 9,-9, 6, 3,-3,-6, 6,-6,-3, 3, 4, 2, 1, 2},
      {-6, 6,-6, 6,-4,-2, 2, 4,-3, 3, 3,-3,-2,-1,-1,-2},
      {2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0},
      {-6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1},
      {4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1}
    };

  int i, j, k, in;
  double xx, x[16];

  for (i = 0; i < 4; i++) {
    x[i] = gs[i];
    x[i+4] = d1gs[i]*CMAPDX;
    x[i+8] = d2gs[i]*CMAPDX;
    x[i+12] = d12gs[i]*CMAPDX*CMAPDX;
  }

  in = 0;
  for (i = 0; i < 4; i++) {
    for (j = 0; j < 4; j++) {
      xx = 0.0;
      for (k = 0; k < 16; k++) xx += wt[in][k]*x[k];
      in++;
      cij[i][j] = xx;
    }
  }
}

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

void FixCMAP::bc_interpol(double x1, double x2, int low1, int low2, double *gs,
                           double *d1gs, double *d2gs, double *d12gs)
{
  // for a given point of interest and its corresponding grid square values,
  //   gradients and cross-derivatives
  // calculate the interpolated value of the point of interest (POI)

  int i;
  double t, u, gs1l, gs2l;

  // set the interpolation coefficients

  bc_coeff(gs,d1gs,d2gs,d12gs);

  gs1l = g_axis[low1];
  gs2l = g_axis[low2];

  t = (x1-gs1l)/CMAPDX;
  u = (x2-gs2l)/CMAPDX;

  E = dEdPhi = dEdPsi = 0.0;

  for (i = 3; i >= 0; i--) {
    E = t*E + ((cij[i][3]*u+cij[i][2])*u+cij[i][1])*u+cij[i][0];
    dEdPhi = u*dEdPhi + (3.0*cij[3][i]*t+2.0*cij[2][i])*t+cij[1][i];
    dEdPsi = t*dEdPsi + (3.0*cij[i][3]*u+2.0*cij[i][2])*u+cij[i][1];
  }

  dEdPhi *= (180.0/MY_PI/CMAPDX);
  dEdPsi *= (180.0/MY_PI/CMAPDX);
}

// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
// methods to read and write data file
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------

void FixCMAP::read_data_header(char *line)
{
  if (strstr(line,"crossterms")) {
    sscanf(line,BIGINT_FORMAT,&ncmap);
  } else error->all(FLERR,"Invalid read data header line for fix cmap");

  // didn't set in constructor b/c this fix could be defined
  // before newton command

  newton_bond = force->newton_bond;
}

/* ----------------------------------------------------------------------
   unpack N lines in buf from section of data file labeled by keyword
   id_offset is applied to atomID fields if multiple data files are read
   store CMAP interactions as if newton_bond = OFF, even if actually ON
------------------------------------------------------------------------- */

void FixCMAP::read_data_section(char *keyword, int n, char *buf,
                                 tagint id_offset)
{
  int m,tmp,itype;
  tagint atom1,atom2,atom3,atom4,atom5;
  char *next;

  next = strchr(buf,'\n');
  *next = '\0';
  int nwords = atom->count_words(buf);
  *next = '\n';

  if (nwords != 7) {
    char str[128];
    snprintf(str,128,"Incorrect %s format in data file",keyword);
    error->all(FLERR,str);
  }

  // loop over lines of CMAP crossterms
  // tokenize the line into values
  // add crossterm to one of my atoms, depending on newton_bond

  for (int i = 0; i < n; i++) {
    next = strchr(buf,'\n');
    *next = '\0';
    sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT
           " " TAGINT_FORMAT " " TAGINT_FORMAT,
           &tmp,&itype,&atom1,&atom2,&atom3,&atom4,&atom5);

    atom1 += id_offset;
    atom2 += id_offset;
    atom3 += id_offset;
    atom4 += id_offset;
    atom5 += id_offset;

    if ((m = atom->map(atom1)) >= 0) {
      if (num_crossterm[m] == CMAPMAX)
        error->one(FLERR,"Too many CMAP crossterms for one atom");
      crossterm_type[m][num_crossterm[m]] = itype;
      crossterm_atom1[m][num_crossterm[m]] = atom1;
      crossterm_atom2[m][num_crossterm[m]] = atom2;
      crossterm_atom3[m][num_crossterm[m]] = atom3;
      crossterm_atom4[m][num_crossterm[m]] = atom4;
      crossterm_atom5[m][num_crossterm[m]] = atom5;
      num_crossterm[m]++;
    }

    if ((m = atom->map(atom2)) >= 0) {
      if (num_crossterm[m] == CMAPMAX)
        error->one(FLERR,"Too many CMAP crossterms for one atom");
      crossterm_type[m][num_crossterm[m]] = itype;
      crossterm_atom1[m][num_crossterm[m]] = atom1;
      crossterm_atom2[m][num_crossterm[m]] = atom2;
      crossterm_atom3[m][num_crossterm[m]] = atom3;
      crossterm_atom4[m][num_crossterm[m]] = atom4;
      crossterm_atom5[m][num_crossterm[m]] = atom5;
      num_crossterm[m]++;
    }

    if ((m = atom->map(atom3)) >= 0) {
      if (num_crossterm[m] == CMAPMAX)
        error->one(FLERR,"Too many CMAP crossterms for one atom");
      crossterm_type[m][num_crossterm[m]] = itype;
      crossterm_atom1[m][num_crossterm[m]] = atom1;
      crossterm_atom2[m][num_crossterm[m]] = atom2;
      crossterm_atom3[m][num_crossterm[m]] = atom3;
      crossterm_atom4[m][num_crossterm[m]] = atom4;
      crossterm_atom5[m][num_crossterm[m]] = atom5;
      num_crossterm[m]++;
    }

    if ((m = atom->map(atom4)) >= 0) {
      if (num_crossterm[m] == CMAPMAX)
        error->one(FLERR,"Too many CMAP crossterms for one atom");
      crossterm_type[m][num_crossterm[m]] = itype;
      crossterm_atom1[m][num_crossterm[m]] = atom1;
      crossterm_atom2[m][num_crossterm[m]] = atom2;
      crossterm_atom3[m][num_crossterm[m]] = atom3;
      crossterm_atom4[m][num_crossterm[m]] = atom4;
      crossterm_atom5[m][num_crossterm[m]] = atom5;
      num_crossterm[m]++;
    }

    if ((m = atom->map(atom5)) >= 0) {
      if (num_crossterm[m] == CMAPMAX)
        error->one(FLERR,"Too many CMAP crossterms for one atom");
      crossterm_type[m][num_crossterm[m]] = itype;
      crossterm_atom1[m][num_crossterm[m]] = atom1;
      crossterm_atom2[m][num_crossterm[m]] = atom2;
      crossterm_atom3[m][num_crossterm[m]] = atom3;
      crossterm_atom4[m][num_crossterm[m]] = atom4;
      crossterm_atom5[m][num_crossterm[m]] = atom5;
      num_crossterm[m]++;
    }

    buf = next + 1;
  }
}

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

bigint FixCMAP::read_data_skip_lines(char * /*keyword*/)
{
  return ncmap;
}

/* ----------------------------------------------------------------------
   write Mth header line to file
   only called by proc 0
------------------------------------------------------------------------- */

void FixCMAP::write_data_header(FILE *fp, int /*mth*/)
{
  fprintf(fp,BIGINT_FORMAT " cmap crossterms\n",ncmap);
}

/* ----------------------------------------------------------------------
   return size I own for Mth data section
   # of data sections = 1 for this fix
   nx = # of crossterms owned by my local atoms
     if newton_bond off, atom only owns crossterm if it is atom3
   ny = columns = type + 5 atom IDs
------------------------------------------------------------------------- */

void FixCMAP::write_data_section_size(int /*mth*/, int &nx, int &ny)
{
  int i,m;

  tagint *tag = atom->tag;
  int nlocal = atom->nlocal;

  nx = 0;
  for (i = 0; i < nlocal; i++)
    for (m = 0; m < num_crossterm[i]; m++)
      if (crossterm_atom3[i][m] == tag[i]) nx++;

  ny = 6;
}

/* ----------------------------------------------------------------------
   pack values for Mth data section into 2d buf
   buf allocated by caller as owned crossterms by 6
------------------------------------------------------------------------- */

void FixCMAP::write_data_section_pack(int /*mth*/, double **buf)
{
  int i,m;

  // 1st column = CMAP type
  // 2nd-6th columns = 5 atom IDs

  tagint *tag = atom->tag;
  int nlocal = atom->nlocal;

  int n = 0;
  for (i = 0; i < nlocal; i++) {
    for (m = 0; m < num_crossterm[i]; m++) {
      if (crossterm_atom3[i][m] != tag[i]) continue;
      buf[n][0] = ubuf(crossterm_type[i][m]).d;
      buf[n][1] = ubuf(crossterm_atom1[i][m]).d;
      buf[n][2] = ubuf(crossterm_atom2[i][m]).d;
      buf[n][3] = ubuf(crossterm_atom3[i][m]).d;
      buf[n][4] = ubuf(crossterm_atom4[i][m]).d;
      buf[n][5] = ubuf(crossterm_atom5[i][m]).d;
      n++;
    }
  }
}

/* ----------------------------------------------------------------------
   write section keyword for Mth data section to file
   use Molecules or Charges if that is only field, else use fix ID
   only called by proc 0
------------------------------------------------------------------------- */

void FixCMAP::write_data_section_keyword(int /*mth*/, FILE *fp)
{
  fprintf(fp,"\nCMAP\n\n");
}

/* ----------------------------------------------------------------------
   write N lines from buf to file
   convert buf fields to int or double depending on styles
   index can be used to prepend global numbering
   only called by proc 0
------------------------------------------------------------------------- */

void FixCMAP::write_data_section(int /*mth*/, FILE *fp,
                                  int n, double **buf, int index)
{
  for (int i = 0; i < n; i++)
    fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
            " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT "\n",
            index+i,(int) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
            (tagint) ubuf(buf[i][2]).i,(tagint) ubuf(buf[i][3]).i,
            (tagint) ubuf(buf[i][4]).i,(tagint) ubuf(buf[i][5]).i);
}

// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
// methods for restart and communication
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------

/* ----------------------------------------------------------------------
   pack entire state of Fix into one write
------------------------------------------------------------------------- */

void FixCMAP::write_restart(FILE *fp)
{
  if (comm->me == 0) {
    int size = sizeof(bigint);
    fwrite(&size,sizeof(int),1,fp);
    fwrite(&ncmap,sizeof(bigint),1,fp);
  }
}

/* ----------------------------------------------------------------------
   use state info from restart file to restart the Fix
------------------------------------------------------------------------- */

void FixCMAP::restart(char *buf)
{
  ncmap = *((bigint *) buf);
}

/* ----------------------------------------------------------------------
   pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */

int FixCMAP::pack_restart(int i, double *buf)
{
  int n = 1;
  for (int m = 0; m < num_crossterm[i]; m++) {
    buf[n++] = ubuf(MAX(crossterm_type[i][m],-crossterm_type[i][m])).d;
    buf[n++] = ubuf(crossterm_atom1[i][m]).d;
    buf[n++] = ubuf(crossterm_atom2[i][m]).d;
    buf[n++] = ubuf(crossterm_atom3[i][m]).d;
    buf[n++] = ubuf(crossterm_atom4[i][m]).d;
    buf[n++] = ubuf(crossterm_atom5[i][m]).d;
  }
  buf[0] = n;

  return n;
}

/* ----------------------------------------------------------------------
   unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */

void FixCMAP::unpack_restart(int nlocal, int nth)
{
  double **extra = atom->extra;

  // skip to Nth set of extra values

   int n = 0;
   for (int i = 0; i < nth; i++) n += static_cast<int> (extra[nlocal][n]);

   int count = static_cast<int> (extra[nlocal][n++]);
   num_crossterm[nlocal] = (count-1)/6;

   for (int m = 0; m < num_crossterm[nlocal]; m++) {
     crossterm_type[nlocal][m] = (int) ubuf(extra[nlocal][n++]).i;
     crossterm_atom1[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
     crossterm_atom2[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
     crossterm_atom3[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
     crossterm_atom4[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
     crossterm_atom5[nlocal][m] = (tagint) ubuf(extra[nlocal][n++]).i;
   }
}

/* ----------------------------------------------------------------------
   maxsize of any atom's restart data
------------------------------------------------------------------------- */

int FixCMAP::maxsize_restart()
{
  return 1 + CMAPMAX*6;
}

/* ----------------------------------------------------------------------
   size of atom nlocal's restart data
------------------------------------------------------------------------- */

int FixCMAP::size_restart(int nlocal)
{
  return 1 + num_crossterm[nlocal]*6;
}

/* ----------------------------------------------------------------------
   allocate atom-based array
------------------------------------------------------------------------- */

void FixCMAP::grow_arrays(int nmax)
{
  num_crossterm = memory->grow(num_crossterm,nmax,"cmap:num_crossterm");
  crossterm_type = memory->grow(crossterm_type,nmax,CMAPMAX,
                                "cmap:crossterm_type");
  crossterm_atom1 = memory->grow(crossterm_atom1,nmax,CMAPMAX,
                                 "cmap:crossterm_atom1");
  crossterm_atom2 = memory->grow(crossterm_atom2,nmax,CMAPMAX,
                                 "cmap:crossterm_atom2");
  crossterm_atom3 = memory->grow(crossterm_atom3,nmax,CMAPMAX,
                                 "cmap:crossterm_atom3");
  crossterm_atom4 = memory->grow(crossterm_atom4,nmax,CMAPMAX,
                                 "cmap:crossterm_atom4");
  crossterm_atom5 = memory->grow(crossterm_atom5,nmax,CMAPMAX,
                                 "cmap:crossterm_atom5");

  // must initialize num_crossterm to 0 for added atoms
  // may never be set for some atoms when data file is read

  for (int i = nmax_previous; i < nmax; i++) num_crossterm[i] = 0;
  nmax_previous = nmax;
}

/* ----------------------------------------------------------------------
   copy values within local atom-based array
------------------------------------------------------------------------- */

void FixCMAP::copy_arrays(int i, int j, int /*delflag*/)
{
  num_crossterm[j] = num_crossterm[i];

  for (int k = 0; k < num_crossterm[j]; k++){
    crossterm_type[j][k] = crossterm_type[i][k];
    crossterm_atom1[j][k] = crossterm_atom1[i][k];
    crossterm_atom2[j][k] = crossterm_atom2[i][k];
    crossterm_atom3[j][k] = crossterm_atom3[i][k];
    crossterm_atom4[j][k] = crossterm_atom4[i][k];
    crossterm_atom5[j][k] = crossterm_atom5[i][k];
  }
}

/* ----------------------------------------------------------------------
   initialize one atom's array values, called when atom is created
------------------------------------------------------------------------- */

void FixCMAP::set_arrays(int i)
{
  num_crossterm[i] = 0;
}

/* ----------------------------------------------------------------------
   pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */

int FixCMAP::pack_exchange(int i, double *buf)
{
  int n = 0;
  buf[n++] = ubuf(num_crossterm[i]).d;
  for (int m = 0; m < num_crossterm[i]; m++) {
    buf[n++] = ubuf(crossterm_type[i][m]).d;
    buf[n++] = ubuf(crossterm_atom1[i][m]).d;
    buf[n++] = ubuf(crossterm_atom2[i][m]).d;
    buf[n++] = ubuf(crossterm_atom3[i][m]).d;
    buf[n++] = ubuf(crossterm_atom4[i][m]).d;
    buf[n++] = ubuf(crossterm_atom5[i][m]).d;
  }
  return n;
}

/* ----------------------------------------------------------------------
   unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */

int FixCMAP::unpack_exchange(int nlocal, double *buf)
{
  int n = 0;
  num_crossterm[nlocal] = (int) ubuf(buf[n++]).i;
  for (int m = 0; m < num_crossterm[nlocal]; m++) {
    crossterm_type[nlocal][m] = (int) ubuf(buf[n++]).i;
    crossterm_atom1[nlocal][m] = (tagint) ubuf(buf[n++]).i;
    crossterm_atom2[nlocal][m] = (tagint) ubuf(buf[n++]).i;
    crossterm_atom3[nlocal][m] = (tagint) ubuf(buf[n++]).i;
    crossterm_atom4[nlocal][m] = (tagint) ubuf(buf[n++]).i;
    crossterm_atom5[nlocal][m] = (tagint) ubuf(buf[n++]).i;
  }
  return n;
}

/* ----------------------------------------------------------------------
   memory usage of local atom-based arrays
------------------------------------------------------------------------- */

double FixCMAP::memory_usage()
{
  int nmax = atom->nmax;
  double bytes = nmax * sizeof(int);        // num_crossterm
  bytes += nmax*CMAPMAX * sizeof(int);      // crossterm_type
  bytes += 5*nmax*CMAPMAX * sizeof(int);    // crossterm_atom 12345
  bytes += maxcrossterm*6 * sizeof(int);    // crosstermlist
  return bytes;
}