Commit a3ca177d authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

various small changes and reformatting 

- add memory usage estimator
- test against varying number of atoms
- test against non-consecutive atom tags
- test for 32-bit overflow in number of atoms
- test for 32-bit overflow in timestep
- reduce tail correction error to warning
- more LAMMPS style formatting of the source code
- remove trailing whitespace
- avoid leaking memory from allocated arrays for masses/charges/tags
parent a3c0fe77

src/USER-PLUMED/fix_plumed.cpp

+108 −48
Original line number Diff line number Diff line
@@ -12,7 +12,8 @@ 
------------------------------------------------------------------------- */



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

   Contributing author: Gareth Tribello (Queens U, Belfast)

   Contributing authors: Gareth Tribello (Queens U, Belfast)

                         Pablo Piaggi (EPFL)

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



#include <cmath>
@@ -56,11 +57,13 @@ FixPlumed::FixPlumed(LAMMPS *lmp, int narg, char **arg) : 
  if (!atom->tag_enable)

    error->all(FLERR,"Fix plumed requires atom tags");



  if (atom->tag_consecutive() == 0)

    error->all(FLERR,"Fix plumed requires consecutive atom IDs");



  if (igroup != 0 && comm->me == 0)

    error->warning(FLERR,"Fix group for fix plumed is not 'all'. "

                   "Group will be ignored.");





  p=new PLMD::Plumed;



  // Check API version
@@ -106,29 +109,53 @@ FixPlumed::FixPlumed(LAMMPS *lmp, int narg, char **arg) : 
    // LAMMPS units lj

    p->cmd("setNaturalUnits");

  } else {

    // Conversion factor from LAMMPS energy units to kJ mol-1 (units of PLUMED) 



    // Conversion factor from LAMMPS energy units to kJ/mol (units of PLUMED)



    double energyUnits=1.0;



    // LAMMPS units real :: kcal/mol;

    if (strcmp(update->unit_style,"real") == 0) energyUnits=4.184;



    if (strcmp(update->unit_style,"real") == 0) {

      energyUnits=4.184;



      // LAMMPS units metal :: eV;

    else if (strcmp(update->unit_style,"metal") == 0) energyUnits=96.48530749925792;



    } else if (strcmp(update->unit_style,"metal") == 0) {

      energyUnits=96.48530749925792;



      // LAMMPS units si :: Joule;

    else if (strcmp(update->unit_style,"si") == 0) energyUnits=0.001;



    } else if (strcmp(update->unit_style,"si") == 0) {

      energyUnits=0.001;



      // LAMMPS units cgs :: erg;

    else if (strcmp(update->unit_style,"cgs") == 0) energyUnits=6.0221418e13;



    } else if (strcmp(update->unit_style,"cgs") == 0) {

      energyUnits=6.0221418e13;



      // LAMMPS units electron :: Hartree;

    else if (strcmp(update->unit_style,"electron") == 0) energyUnits=2625.5257;

    else error->all(FLERR,"Odd LAMMPS units, plumed cannot work with that");



    } else if (strcmp(update->unit_style,"electron") == 0) {

      energyUnits=2625.5257;



    } else error->all(FLERR,"Fix plumed cannot handle your choice of units");



    // Conversion factor from LAMMPS length units to nm (units of PLUMED)



    double lengthUnits=0.1/force->angstrom;



    // Conversion factor from LAMMPS time unit to ps (units of PLUMED)



    double timeUnits=0.001/force->femtosecond;



    p->cmd("setMDEnergyUnits",&energyUnits);

    p->cmd("setMDLengthUnits",&lengthUnits);

    p->cmd("setMDTimeUnits",&timeUnits);

  }



  // Read fix parameters:



  int next=0;

  for (int i=3;i<narg;++i) {

    if (!strcmp(arg[i],"outfile")) {
@@ -161,7 +188,10 @@ FixPlumed::FixPlumed(LAMMPS *lmp, int narg, char **arg) : 


  p->cmd("setMDEngine","LAMMPS");



  int natoms=int(atom->natoms);

  if (atom->natoms > MAXSMALLINT)

    error->all(FLERR,"Fix plumed can only handle up to 2.1 billion atoms");



  natoms=int(atom->natoms);

  p->cmd("setNatoms",&natoms);



  double dt=update->dt;
@@ -232,6 +262,9 @@ FixPlumed::~FixPlumed() 
  delete p;

  modify->delete_compute(id_pe);

  modify->delete_compute(id_press);

  delete[] masses;

  delete[] charges;

  delete[] gatindex;

}



int FixPlumed::setmask()
@@ -249,8 +282,10 @@ void FixPlumed::init() 
{

  if (strcmp(update->integrate_style,"respa") == 0)

    nlevels_respa = ((Respa *) update->integrate)->nlevels;



  // This avoids nan pressure if compute_pressure is called

  // in a setup method



  for (int i=0;i<6;i++) virial[i] = 0.;

}
@@ -283,21 +318,29 @@ void FixPlumed::min_setup(int vflag) 


void FixPlumed::post_force(int /* vflag */)

{

// Check tag is enabled

  if ( !atom->tag_enable ) error->all(FLERR,"to run PLUMED you must have tag_enable==1");



  int update_gatindex=0;



  if (natoms != int(atom->natoms))

    error->all(FLERR,"Fix plumed does not support simulations with varying "

               "numbers of atoms");



  // Try to find out if the domain decomposition has been updated:



  if (nlocal != atom->nlocal) {



    if (charges) delete [] charges;

    if (masses) delete [] masses;

    if (gatindex) delete [] gatindex;



    nlocal=atom->nlocal;

    gatindex=new int [nlocal];

    masses=new double [nlocal];

    charges=new double [nlocal];

    update_gatindex=1;



  } else {



    for (int i=0;i<nlocal;i++) {

      if (gatindex[i]!=atom->tag[i]-1) {

        update_gatindex=1;
@@ -309,6 +352,7 @@ void FixPlumed::post_force(int /* vflag */) 


  // In case it has been updated, rebuild the local mass/charges array

  // and tell plumed about the change:



  if (update_gatindex) {

    for (int i=0;i<nlocal;i++) gatindex[i]=atom->tag[i]-1;

    // Get masses
@@ -339,12 +383,16 @@ void FixPlumed::post_force(int /* vflag */) 
  box[2][1]=domain->h[3];

  box[2][0]=domain->h[4];

  box[1][0]=domain->h[5];



  // Make initial of virial of this fix zero

  // The following line is very important, otherwise the compute pressure will include 

  // The following line is very important, otherwise

  // the compute pressure will include

  for (int i=0;i<6;++i) virial[i] = 0.;



  // local variable with timestep:

  int step=update->ntimestep;

  if (update->ntimestep > MAXSMALLINT)

    error->all(FLERR,"Fix plumed can only handle up to 2.1 billion timesteps");

  int step=int(update->ntimestep);



  // pass all pointers to plumed:

  p->cmd("setStep",&step);
@@ -352,8 +400,9 @@ void FixPlumed::post_force(int /* vflag */) 
  p->cmd("setBox",&box[0][0]);

  p->cmd("setForces",&atom->f[0][0]);

  p->cmd("setMasses",&masses[0]);

  if (atom->q) p->cmd("setCharges",&charges[0]);

  p->cmd("setCharges",&charges[0]);

  p->cmd("getBias",&bias);



  // Pass virial to plumed

  // If energy is needed virial_plmd is equal to Lammps' virial

  // If energy is not needed virial_plmd is initialized to zero
@@ -369,30 +418,37 @@ void FixPlumed::post_force(int /* vflag */) 
  double *virial_lmp;

  if (plumedNeedsEnergy) {

    // Error if tail corrections are included

    if (force->pair && force->pair->tail_flag)

      error->all(FLERR,"Tail corrections to the pair potential included."

                 " The energy cannot be biased in this case."

    if (force->pair && force->pair->tail_flag && comm->me == 0)

      error->warning(FLERR,"Tail corrections to the pair potential included."

                     " The energy cannot be biased correctly in this case."

                     " Remove the tail corrections by removing the"

                     " command: pair_modify tail yes");



    // compute the potential energy

    double pot_energy = 0.;

    c_pe->compute_scalar();

    pot_energy = c_pe->scalar;



    // Divide energy by number of processes

    // Plumed wants it this way

    int nprocs;

    MPI_Comm_size(world,&nprocs);

    pot_energy /= nprocs;

    p->cmd("setEnergy",&pot_energy);



    // Compute pressure due to the virial (no kinetic energy term!)

    c_press->compute_vector();

    virial_lmp = c_press->vector;



    // Check if pressure is finite

    if (!std::isfinite(virial_lmp[0]) || !std::isfinite(virial_lmp[1]) || !std::isfinite(virial_lmp[2])

        || !std::isfinite(virial_lmp[3]) || !std::isfinite(virial_lmp[4]) || !std::isfinite(virial_lmp[5]))

    if (!std::isfinite(virial_lmp[0]) || !std::isfinite(virial_lmp[1])

        || !std::isfinite(virial_lmp[2]) || !std::isfinite(virial_lmp[3])

        || !std::isfinite(virial_lmp[4]) || !std::isfinite(virial_lmp[5]))

      error->all(FLERR,"Non-numeric virial - Plumed cannot work with that");



    // Convert pressure to virial per number of MPI processes

    // From now on all virials are divided by the number of MPI processes



    double nktv2p = force->nktv2p;

    double inv_volume;

    if (domain->dimension == 3) {
@@ -415,12 +471,12 @@ void FixPlumed::post_force(int /* vflag */) 
  // do the real calculation:

  p->cmd("performCalc");



// retransform virial to lammps representation and assign it to this fix's virial.

// Plumed is giving back the full virial and therefore we have to subtract the

// initial virial i.e. virial_lmp. 

// The vector virial contains only the contribution added by plumed

// The calculation of the pressure will be done by a compute pressure and will

// include this contribution. 

  // retransform virial to lammps representation and assign it to this

  // fix's virial.  Plumed is giving back the full virial and therefore

  // we have to subtract the initial virial i.e. virial_lmp.

  // The vector virial contains only the contribution added by plumed.

  // The calculation of the pressure will be done by a compute pressure

  // and will include this contribution.

  virial[0] = -plmd_virial[0][0]-virial_lmp[0];

  virial[1] = -plmd_virial[1][1]-virial_lmp[1];

  virial[2] = -plmd_virial[2][2]-virial_lmp[2];
@@ -448,7 +504,7 @@ void FixPlumed::min_post_force(int vflag) 


void FixPlumed::reset_dt()

{

  error->all(FLERR,"cannot reset_dt within a fix plumed action"); 

  error->all(FLERR,"Cannot change the time step when fix plumed is active");

}



double FixPlumed::compute_scalar()
@@ -497,3 +553,7 @@ int FixPlumed::modify_param(int narg, char **arg) 
  return 0;

}



double FixPlumed::memory_usage()

{

  return double((8+8+4)*atom->nlocal);

}

src/USER-PLUMED/fix_plumed.h

+2 −0
Original line number Diff line number Diff line
@@ -43,11 +43,13 @@ class FixPlumed : public Fix { 
  double compute_scalar();

  void reset_dt();

  int modify_param(int narg, char **arg);

  double memory_usage();



 private:



  PLMD::Plumed *p;         // pointer to plumed object

  int nlocal;              // number of atoms local to this process

  int natoms;              // total number of atoms

  int *gatindex;           // array of atom indexes local to this process

  double *masses;          // array of masses for local atoms

  double *charges;         // array of charges for local atoms

CRA Git | Maintained and supported by SUSTech CRA and CCSE