Commit 735ec9de authored by Steven J. Plimpton's avatar Steven J. Plimpton
Browse files

more general MD c/s protocol

parent 33d5fe45
Loading
Loading
Loading
Loading
+34 −9
Original line number Diff line number Diff line
@@ -28,9 +28,9 @@ using namespace LAMMPS_NS;
using namespace CSLIB_NS;
using namespace FixConst;

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

/* ---------------------------------------------------------------------- */
@@ -47,7 +47,7 @@ FixClientMD::FixClientMD(LAMMPS *lmp, int narg, char **arg) :

  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");
  else units = OTHER;

  scalar_flag = 1;
  global_freq = 1;
@@ -55,18 +55,19 @@ FixClientMD::FixClientMD(LAMMPS *lmp, int narg, char **arg) :
  virial_flag = 1;
  thermo_virial = 1;

  inv_nprocs = 1.0 / comm->nprocs;
  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;
  // unit conversion factors for REAL
  // otherwise not needed

  fconvert = econvert = vconvert = 1.0;
  if (units == REAL) {
    // NOTE: still need to set these
    fconvert = 1.0;
    econvert = 1.0;
    vconvert = 1.0;
@@ -120,10 +121,12 @@ void FixClientMD::setup(int vflag)
{
  CSlib *cs = (CSlib *) lmp->cslib;

  // SETUP send at beginning of each run
  // required fields: DIM, PERIODICTY, ORIGIN, BOX, NATOMS, NTYPES, TYPES, COORDS
  // optional fields: others in enum above

  int nfields = 8;
  if (units == OTHER) nfields++;
  if (atom->q_flag) nfields++;

  cs->send(SETUP,nfields);
@@ -142,12 +145,20 @@ void FixClientMD::setup(int vflag)
  pack_coords();
  cs->pack_parallel(COORDS,4,atom->nlocal,atom->tag,3,xpbc);

  if (units == OTHER) cs->pack_string(UNITS,update->unit_style);

  if (atom->q_flag)
    cs->pack_parallel(CHARGE,4,atom->nlocal,atom->tag,1,atom->q);

  // receive initial forces, energy, virial

  receive_fev(vflag);

  if (server_error) {
    str = char[64];
    sprintf(str,"Fix client/md received server error %d\n",server_error);
    error->all(FLERR,str);
  }
}

/* ---------------------------------------------------------------------- */
@@ -168,6 +179,7 @@ void FixClientMD::post_force(int vflag)
  if (vflag) v_setup(vflag);
  else evflag = 0;

  // STEP send every step
  // required fields: COORDS
  // optional fields: ORIGIN, BOX

@@ -189,9 +201,15 @@ void FixClientMD::post_force(int vflag)
    cs->pack(BOX,4,9,&box[0][0]);
  }

  // recv forces, energy, virial
  // receive forces, energy, virial

  receive_fev(vflag);

  if (server_error) {
    str = char[64];
    sprintf(str,"Fix client/md received server error %d\n",server_error);
    error->all(FLERR,str);
  }
}

/* ---------------------------------------------------------------------- */
@@ -202,7 +220,7 @@ void FixClientMD::min_post_force(int vflag)
}

/* ----------------------------------------------------------------------
   potential energy from QM code
   potential energy from server
------------------------------------------------------------------------- */

double FixClientMD::compute_scalar()
@@ -255,7 +273,9 @@ void FixClientMD::pack_box()
}

/* ----------------------------------------------------------------------
   receive message from server with forces, energy, virial
   receive message from server
   required fields: FORCES, ENERGY, VIRIAL
   optional field: ERROR
------------------------------------------------------------------------- */

void FixClientMD::receive_fev(int vflag)
@@ -291,4 +311,9 @@ void FixClientMD::receive_fev(int vflag)
    for (int i = 0; i < 6; i++)
      virial[i] = inv_nprocs * vconvert * v[i];
  }

  // error return

  server_error = 0;
  if (nfield == 4) server_error = cs->unpack_int(ERROR);
}
+62 −24
Original line number Diff line number Diff line
@@ -33,21 +33,11 @@
using namespace LAMMPS_NS;
using namespace CSLIB_NS;

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

// NOTE: features that could be added to this interface
// allow client to set periodicity vs shrink-wrap
//   currently just assume server is same as client
// test that triclinic boxes actually work
// send new box size/shape every step, for NPT client
// unit check between client/server with unit conversion if needed
// option for client to send other per-atom quantities, e.g. rmass
// more precise request of energy/virial (global or peratom) by client
//   maybe Verlet should have a single(eflag,vflag) method to more easily comply

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

ServerMD::ServerMD(LAMMPS *lmp) : Pointers(lmp)
@@ -61,7 +51,27 @@ ServerMD::ServerMD(LAMMPS *lmp) : Pointers(lmp)

  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");
  else units = OTHER;

  // unit conversion factors for REAL
  // otherwise not needed

  fconvert = econvert = vconvert = 1.0;
  if (units == REAL) {
    // NOTE: still need to set these
    fconvert = 1.0;
    econvert = 1.0;
    vconvert = 1.0;
  }

  fcopy = NULL;
}

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

ServerMD::~ServerMD();
{
  memory->destroy(fcopy);
}

/* ---------------------------------------------------------------------- */
@@ -89,8 +99,9 @@ void ServerMD::loop()
    msgID = cs->recv(nfield,fieldID,fieldtype,fieldlen);
    if (msgID < 0) break;

    // SETUP call at beginning of each run
    // required fields: NATOMS, NTYPES, BOXLO, BOXHI, TYPES, COORDS
    // SETUP receive at beginning of each run
    // required fields: DIM, PERIODICTY, ORIGIN, BOX, 
    //                  NATOMS, NTYPES, TYPES, COORDS
    // optional fields: others in enum above

    if (msgID == SETUP) {
@@ -103,6 +114,7 @@ void ServerMD::loop()
      double *box = NULL;
      int *types = NULL;
      double *coords = NULL;
      char *units = NULL;
      double *charge = NULL;

      for (int ifield = 0; ifield < nfield; ifield++) {
@@ -132,6 +144,9 @@ void ServerMD::loop()
          types = (int *) cs->unpack(TYPES);
        } else if (fieldID[ifield] == COORDS) {
          coords = (double *) cs->unpack(COORDS);

        } else if (fieldID[ifield] == UNITS) {
          units = (char *) cs->unpack(UNITS);
        } else if (fieldID[ifield] == CHARGE) {
          charge = (double *) cs->unpack(CHARGE);
        } else error->all(FLERR,"Server md setup field unknown");
@@ -141,6 +156,9 @@ void ServerMD::loop()
          natoms < 0 || ntypes < 0 || !types || !coords)
        error->all(FLERR,"Required server md setup field not received");

      if (units && strcmp(units,update->unit_style) != 0)
        error->all(FLERR,"Server md does not match client units");

      if (charge && atom->q_flag == 0)
        error->all(FLERR,"Server md does not define atom attribute q");

@@ -183,6 +201,13 @@ void ServerMD::loop()
      atom->map_set();
      atom->natoms = natoms;

      // allocate fcopy if needed

      if (units == REAL) {
        memory->destroy(fcopy);
        memory->create(fcopy,atom->nlocal,3,"server/md:fcopy");
      }

      // perform system setup() for dynamics
      // also OK for minimization, since client runs minimizer
      // return forces, energy, virial to client
@@ -195,7 +220,7 @@ void ServerMD::loop()

      send_fev(msgID);

    // STEP call at each timestep of run or minimization
    // STEP receive at each timestep of run or minimization
    // required fields: COORDS
    // optional fields: BOXLO, BOXHI, BOXTILT

@@ -309,13 +334,13 @@ void ServerMD::box_change(double *origin, double *box)
  domain->xy = box[3];
  domain->xz = box[6];
  domain->yz = box[7];


}

/* ----------------------------------------------------------------------
   send return message with forces, energy, pressure tensor
   pressure tensor should be just pair style virial
   return message with forces, energy, pressure tensor
   pressure tensor should be just pair and KSpace contributions
   required fields: FORCES, ENERGY, VIRIAL
   optional field: ERROR (not ever sending)
------------------------------------------------------------------------- */

void ServerMD::send_fev(int msgID)
@@ -325,12 +350,25 @@ void ServerMD::send_fev(int msgID)
  cs->send(msgID,3);
  
  double *forces = NULL;
  if (atom->nlocal) forces = &atom->f[0][0];
  if (atom->nlocal) {
    if (units != REAL) forces = &atom->f[0][0];
    else {
      double **f = atom->f;
      int nlocal = atom->nlocal;
      for (i = 0; i < nlocal; i++) {
        fcopy[i][0] = fconvert*f[i][0];
        fcopy[i][1] = fconvert*f[i][1];
        fcopy[i][2] = fconvert*f[i][2];
      }
      forces = &fcopy[0][0];
    }
  }
  cs->pack_parallel(FORCES,4,atom->nlocal,atom->tag,3,forces);

  double eng = force->pair->eng_vdwl + force->pair->eng_coul;
  double engall;
  MPI_Allreduce(&eng,&engall,1,MPI_DOUBLE,MPI_SUM,world);
  engall *= econvert;
  cs->pack_double(ENERGY,engall);
  
  double v[6],vall[6];
@@ -340,7 +378,7 @@ void ServerMD::send_fev(int msgID)
    for (int i = 0; i < 6; i++)
      v[i] += force->kspace->virial[i];

  for (int i = 0; i < 6; i++) v[i] *= vconvert;
  MPI_Allreduce(&v,&vall,6,MPI_DOUBLE,MPI_SUM,world);
  cs->pack(VIRIAL,4,6,vall);
}
+4 −0
Original line number Diff line number Diff line
@@ -21,9 +21,13 @@ namespace LAMMPS_NS {
class ServerMD : protected Pointers {
 public:
  ServerMD(class LAMMPS *);
  ~ServerMD();
  void loop();

 private:
  int units;
  double **fcopy;

  void box_change(double *, double *, double *);
  void send_fev(int);
};