Commit 8930a88e authored by Steven J. Plimpton's avatar Steven J. Plimpton
Browse files

enrich the AIMD interface

parent 03b880a3
Loading
Loading
Loading
Loading
+60 −37
Original line number Diff line number Diff line
@@ -28,8 +28,9 @@ using namespace LAMMPS_NS;
using namespace CSLIB_NS;
using namespace FixConst;

enum{REAL,METAL}
enum{SETUP=1,STEP};
enum{UNITS=1,DIM,NATOMS,NTYPES,BOXLO,BOXHI,BOXTILT,TYPES,COORDS,CHARGE};
enum{DIM=1,PERIODICITY,ORIGIN,BOX,NATOMS,NTYPES,TYPES,COORDS,CHARGE};
enum{FORCES=1,ENERGY,VIRIAL};

/* ---------------------------------------------------------------------- */
@@ -44,14 +45,32 @@ FixClientMD::FixClientMD(LAMMPS *lmp, int narg, char **arg) :
  if (sizeof(tagint) != 4) 
    error->all(FLERR,"Fix client/md requires 4-byte atom IDs");

  if (strcmp(update->unit_style,"real") == 0) units = REAL;
  else if (strcmp(update->unit_style,"metal") == 0) units = METAL;
  else error->all(FLERR,"Units must be real or metal for fix client/md");

  scalar_flag = 1;
  global_freq = 1;
  extscalar = 1;
  virial_flag = 1;
  thermo_virial = 1;

  if (domain->dimension == 2)
    box[0][2] = box[1][2] = box[2][0] = box[2][1] = box[2][2] = 0.0;

  maxatom = 0;
  xpbc = NULL;

  // unit conversion factors

  double inv_nprocs = 1.0 / comm->nprocs;

  fconvert = econvert = vconvert = 1.0;
  if (units == REAL) {
    fconvert = 1.0;
    econvert = 1.0;
    vconvert = 1.0;
  }
}

/* ---------------------------------------------------------------------- */
@@ -101,34 +120,28 @@ void FixClientMD::setup(int vflag)
{
  CSlib *cs = (CSlib *) lmp->cslib;

  // required fields: NATOMS, NTYPES, BOXLO, BOXHI, TYPES, COORDS
  // required fields: DIM, PERIODICTY, ORIGIN, BOX, NATOMS, NTYPES, TYPES, COORDS
  // optional fields: others in enum above

  int nfields = 6;
  if (domain->dimension == 2) nfields++;
  if (domain->triclinic) nfields++;
  int nfields = 8;
  if (atom->q_flag) nfields++;

  cs->send(SETUP,nfields);

  cs->pack_int(DIM,domain->dimension);
  cs->pack(PERIODICITY,1,3,domain->periodicity);

  pack_box();
  cs->pack(ORIGIN,4,3,domain->boxlo);
  cs->pack(BOX,4,9,&box[0][0]);

  cs->pack_int(NATOMS,atom->natoms);
  cs->pack_int(NTYPES,atom->ntypes);
  cs->pack(BOXLO,4,3,domain->boxlo);
  cs->pack(BOXHI,4,3,domain->boxhi);

  cs->pack_parallel(TYPES,1,atom->nlocal,atom->tag,1,atom->type);
  pack_coords();
  cs->pack_parallel(COORDS,4,atom->nlocal,atom->tag,3,xpbc);

  if (domain->dimension == 2) cs->pack_int(DIM,domain->dimension);
  if (domain->triclinic) {
    double boxtilt[3];
    boxtilt[0] = domain->xy;
    if (domain->dimension == 3) {
      boxtilt[1] = domain->xz;
      boxtilt[2] = domain->yz;
    } else boxtilt[1] = boxtilt[2] = 0.0;
    cs->pack(BOXTILT,4,3,boxtilt);
  }
  if (atom->q_flag)
    cs->pack_parallel(CHARGE,4,atom->nlocal,atom->tag,1,atom->q);

@@ -156,7 +169,7 @@ void FixClientMD::post_force(int vflag)
  else evflag = 0;

  // required fields: COORDS
  // optional fields: BOXLO, BOXHI, BOXTILT
  // optional fields: ORIGIN, BOX

  // send coords

@@ -164,7 +177,6 @@ void FixClientMD::post_force(int vflag)

  int nfields = 1;
  if (domain->box_change) nfields += 2;
  if (domain->box_change && domain->triclinic) nfields++;;

  cs->send(STEP,nfields);

@@ -172,17 +184,9 @@ void FixClientMD::post_force(int vflag)
  cs->pack_parallel(COORDS,4,atom->nlocal,atom->tag,3,xpbc);

  if (domain->box_change) {
    cs->pack(BOXLO,4,3,domain->boxlo);
    cs->pack(BOXHI,4,3,domain->boxhi);
    if (domain->triclinic) {
      double boxtilt[3];
      boxtilt[0] = domain->xy;
      if (domain->dimension == 3) {
        boxtilt[1] = domain->xz;
        boxtilt[2] = domain->yz;
      } else boxtilt[1] = boxtilt[2] = 0.0;
      cs->pack(BOXTILT,4,3,boxtilt);
    }
    pack_box();
    cs->pack(ORIGIN,4,3,domain->boxlo);
    cs->pack(BOX,4,9,&box[0][0]);
  }

  // recv forces, energy, virial
@@ -230,6 +234,26 @@ void FixClientMD::pack_coords()
  }
}

/* ----------------------------------------------------------------------
   pack box info into box = 3 edge vectors of simulation box
------------------------------------------------------------------------- */

void FixClientMD::pack_box()
{
  double *boxlo = domain->boxlo;
  double *boxhi = domain->boxhi;

  box[0][0] = boxhi[0] - boxlo[0];
  box[0][1] = 0.0;
  box[0][2] = 0.0;
  box[1][0] = domain->xy;
  box[1][1] = boxhi[1] - boxlo[1];
  box[1][2] = 0.0;
  box[2][0] = domain->xz;
  box[2][1] = domain->yz;
  box[2][2] = boxhi[2] - boxlo[2];
}

/* ----------------------------------------------------------------------
   receive message from server with forces, energy, virial
------------------------------------------------------------------------- */
@@ -254,18 +278,17 @@ void FixClientMD::receive_fev(int vflag)
    m = atom->map(id);
    if (m < 0 || m >= nlocal) j += 3;
    else {
      f[m][0] += forces[j++];
      f[m][1] += forces[j++];
      f[m][2] += forces[j++];
      f[m][0] += fconvert * forces[j++];
      f[m][1] += fconvert * forces[j++];
      f[m][2] += fconvert * forces[j++];
    }
  }

  eng = cs->unpack_double(ENERGY);
  eng = econvert * cs->unpack_double(ENERGY);
  
  if (vflag) {
    double *v = (double *) cs->unpack(VIRIAL);
    double invnprocs = 1.0 / comm->nprocs;
    for (int i = 0; i < 6; i++)
      virial[i] = invnprocs*v[i];
      virial[i] = inv_nprocs * vconvert * v[i];
  }
}
+5 −1
Original line number Diff line number Diff line
@@ -38,10 +38,14 @@ class FixClientMD : public Fix {

 private:
  void *cslib;
  int maxatom;
  int maxatom,units;
  double box[3][3];
  double eng;
  double *xpbc;

  double inv_nprocs;
  double fconvert,econvert,vconvert;

  void pack_coords();
  void receive_fev(int );
};
+69 −56
Original line number Diff line number Diff line
@@ -33,8 +33,9 @@
using namespace LAMMPS_NS;
using namespace CSLIB_NS;

enum{REAL,METAL}
enum{SETUP=1,STEP};
enum{UNITS=1,DIM,NATOMS,NTYPES,BOXLO,BOXHI,BOXTILT,TYPES,COORDS,CHARGE};
enum{DIM=1,PERIODICITY,ORIGIN,BOX,NATOMS,NTYPES,TYPES,COORDS,CHARGE};
enum{FORCES=1,ENERGY,VIRIAL};

// NOTE: features that could be added to this interface
@@ -54,9 +55,13 @@ ServerMD::ServerMD(LAMMPS *lmp) : Pointers(lmp)
  if (domain->box_exist == 0)
    error->all(FLERR,"Server command before simulation box is defined");

  if (!atom->map_style) error->all(FLERR,"Server md mode requires atom map");
  if (atom->tag_enable == 0) error->all(FLERR,"Server md mode requires atom IDs");
  if (!atom->map_style) error->all(FLERR,"Server md requires atom map");
  if (atom->tag_enable == 0) error->all(FLERR,"Server md requires atom IDs");
  if (sizeof(tagint) != 4) error->all(FLERR,"Server md requires 4-byte atom IDs");

  if (strcmp(update->unit_style,"real") == 0) units = REAL;
  else if (strcmp(update->unit_style,"metal") == 0) units = METAL;
  else error->all(FLERR,"Units must be real or metal for server md");
}

/* ---------------------------------------------------------------------- */
@@ -90,24 +95,31 @@ void ServerMD::loop()

    if (msgID == SETUP) {

      int dim = 0;
      int *periodicity = NULL;
      int natoms = -1;
      int ntypes = -1;
      double *boxlo = NULL;
      double *boxhi = NULL;
      double *boxtilt = NULL;
      double *origin = NULL;
      double *box = NULL;
      int *types = NULL;
      double *coords = NULL;
      double *charge = NULL;

      for (int ifield = 0; ifield < nfield; ifield++) {
        if (fieldID[ifield] == UNITS) {
          char *units = cs->unpack_string(UNITS);
          if (strcmp(units,update->unit_style) != 0)
            error->all(FLERR,"Server md units mis-match with client");
        } else if (fieldID[ifield] == DIM) {
          int dim = cs->unpack_int(DIM);
        if (fieldID[ifield] == DIM) {
          dim = cs->unpack_int(DIM);
          if (dim != domain->dimension)
            error->all(FLERR,"Server md dimension mis-match with client");
        } else if (fieldID[ifield] == PERIODICITY) {
          periodicity = (int *) cs->unpack(PERIODICITY);
          if (periodicity[0] != domain->periodicity[0] ||
              periodicity[1] != domain->periodicity[1]
              periodicity[2] != domain->periodicity[2])
            error->all(FLERR,"Server md periodicity mis-match with client");
        } else if (fieldID[ifield] == ORIGIN) {
          origin = (double *) cs->unpack(ORIGIN);
        } else if (fieldID[ifield] == BOX) {
          box = (double *) cs->unpack(BOX);
        } else if (fieldID[ifield] == NATOMS) {
          natoms = cs->unpack_int(NATOMS);
          if (3*natoms > INT_MAX)
@@ -116,12 +128,6 @@ void ServerMD::loop()
          ntypes = cs->unpack_int(NTYPES);
          if (ntypes != atom->ntypes)
            error->all(FLERR,"Server md ntypes mis-match with client");
        } else if (fieldID[ifield] == BOXLO) {
          boxlo = (double *) cs->unpack(BOXLO);
        } else if (fieldID[ifield] == BOXHI) {
          boxhi = (double *) cs->unpack(BOXHI);
        } else if (fieldID[ifield] == BOXTILT) {
          boxtilt = (double *) cs->unpack(BOXTILT);
        } else if (fieldID[ifield] == TYPES) {
          types = (int *) cs->unpack(TYPES);
        } else if (fieldID[ifield] == COORDS) {
@@ -131,7 +137,8 @@ void ServerMD::loop()
        } else error->all(FLERR,"Server md setup field unknown");
      }

      if (natoms < 0 || ntypes < 0 || !boxlo || !boxhi || !types || !coords)
      if (dim == 0 || !periodicity || !origin || !box || 
          natoms < 0 || ntypes < 0 || !types || !coords)
        error->all(FLERR,"Required server md setup field not received");

      if (charge && atom->q_flag == 0)
@@ -140,7 +147,11 @@ void ServerMD::loop()
      // reset box, global and local
      // reset proc decomposition
 
      box_change(boxlo,boxhi,boxtilt);
      if ((box[3] != 0.0 || box[6] != 0.0 || box[7] != 0.0) && 
          domain->triclinic == 0)
        error->all(FLERR,"Server md is not initialized for a triclinic box");

      box_change(origin,box);

      domain->set_initial_box();
      domain->set_global_box();
@@ -166,7 +177,7 @@ void ServerMD::loop()
      int ntotal;
      MPI_Allreduce(&atom->nlocal,&ntotal,1,MPI_INT,MPI_SUM,world);
      if (ntotal != natoms) 
        error->all(FLERR,"Server md atom cound does not match client");
        error->all(FLERR,"Server md atom count does not match client");

      atom->map_init();
      atom->map_set();
@@ -190,29 +201,34 @@ void ServerMD::loop()

    } else if (msgID == STEP) {

      double *boxlo = NULL;
      double *boxhi = NULL;
      double *boxtilt = NULL;
      double *coords = NULL;
      double *origin = NULL;
      double *box = NULL;

      for (int ifield = 0; ifield < nfield; ifield++) {
        if (fieldID[ifield] == BOXLO) {
          boxlo = (double *) cs->unpack(BOXLO);
        } else if (fieldID[ifield] == BOXHI) {
          boxhi = (double *) cs->unpack(BOXHI);
        } else if (fieldID[ifield] == BOXTILT) {
          boxtilt = (double *) cs->unpack(BOXTILT);
        } else if (fieldID[ifield] == COORDS) {
        if (fieldID[ifield] == COORDS) {
          coords = (double *) cs->unpack(COORDS);
        } else if (fieldID[ifield] == ORIGIN) {
          origin = (double *) cs->unpack(ORIGIN);
        } else if (fieldID[ifield] == BOX) {
          box = (double *) cs->unpack(BOX);
        } else error->all(FLERR,"Server md step field unknown");
      }

      if (!coords)
        error->all(FLERR,"Required server md step field not received");

      // change box size/shape, only if both box lo/hi received
      // change box size/shape, only if origin and box received
      // reset global/local box like FixDeform at end_of_step()

      if (boxlo && boxhi) box_change(boxlo,boxhi,boxtilt);
      if (origin && box) {
        if ((box[3] != 0.0 || box[6] != 0.0 || box[7] != 0.0) && 
            domain->triclinic == 0)
          error->all(FLERR,"Server md is not initialized for a triclinic box");
        box_change(origin,box);
        domain->set_global_box();
        domain->set_local_box();
      }

      // assign received coords to owned atoms
      // closest_image() insures xyz matches current server PBC
@@ -231,9 +247,11 @@ void ServerMD::loop()
        }
      }

      // if no need to reneighbor, just ghost comm
      // else setup_minimal(1) which includes reneigh
      // setup_minimal computes forces for 0 or 1
      // if no need to reneighbor:
      //   ghost comm
      //   setup_minimal(0) which just computes forces
      // else:
      //   setup_minimal(1) for pbc, reset_box, reneigh, forces

      int nflag = neighbor->decide();
      if (nflag == 0) {
@@ -278,26 +296,21 @@ void ServerMD::loop()
   box change due to received message
------------------------------------------------------------------------- */

void ServerMD::box_change(double *boxlo, double *boxhi, double *boxtilt)
void ServerMD::box_change(double *origin, double *box)
{
  domain->boxlo[0] = boxlo[0];
  domain->boxhi[0] = boxhi[0];
  domain->boxlo[1] = boxlo[1];
  domain->boxhi[1] = boxhi[1];
  if (domain->dimension == 3) {
    domain->boxlo[2] = boxlo[2];
    domain->boxhi[2] = boxhi[2];
  }
  domain->boxlo[0] = origin[0];
  domain->boxlo[1] = origin[1];
  domain->boxlo[2] = origin[2];

  domain->boxhi[0] = origin[0] + box[0];
  domain->boxhi[1] = origin[1] + box[4];
  domain->boxhi[2] = origin[2] + box[8];

  domain->xy = box[3];
  domain->xz = box[6];
  domain->yz = box[7];


  if (boxtilt) {
    if (!domain->triclinic) 
      error->all(FLERR,"Server md not setup for triclinic box");
    domain->xy = boxtilt[0];
    if (domain->dimension == 3) {
      domain->xz = boxtilt[1];
      domain->yz = boxtilt[2];
    }
  }
}

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