Commit 6956dc1a authored by Agilio Padua's avatar Agilio Padua
Browse files

Added safeguard to thole pair styles

parent 7434267f

src/USER-DRUDE/pair_lj_cut_thole_long.cpp

+11 −97
Original line number Diff line number Diff line
@@ -45,6 +45,10 @@ using namespace MathConst; 
#define B4       -5.80844129e-3

#define B5        1.14652755e-1



#define EPSILON 1.0e-20

#define EPS_EWALD 1.0e-6

#define EPS_EWALD_SQR 1.0e-12



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



PairLJCutTholeLong::PairLJCutTholeLong(LAMMPS *lmp) : Pair(lmp)
@@ -149,6 +153,7 @@ void PairLJCutTholeLong::compute(int eflag, int vflag) 
      jtype = type[j];



      if (rsq < cutsq[itype][jtype]) {

         rsq += EPSILON; // Add Epsilon for case: r = 0; DC-DP 1-1 interaction must be removed by special bond;

         r2inv = 1.0/rsq;



        if (rsq < cut_coulsq) {
@@ -156,14 +161,15 @@ void PairLJCutTholeLong::compute(int eflag, int vflag) 
          r = sqrt(rsq);



          if (!ncoultablebits || rsq <= tabinnersq) {

            grij = g_ewald * r;

            grij = g_ewald * (r + EPS_EWALD);

            expm2 = exp(-grij*grij);

            t = 1.0 / (1.0 + EWALD_P*grij);

            u = 1. - t;

            erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;

            prefactor = qqrd2e * qi*qj/r;

            prefactor = qqrd2e * qi*qj/(r + EPS_EWALD);

            forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

            if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;

            r2inv = 1.0/(rsq + EPS_EWALD_SQR);

          } else {

            union_int_float_t rsq_lookup;

            rsq_lookup.f = rsq;
@@ -579,100 +585,8 @@ void PairLJCutTholeLong::write_data_all(FILE *fp) 


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



double PairLJCutTholeLong::single(int i, int j, int itype, int jtype,

                                 double rsq, double factor_coul,

                                 double factor_lj, double &fforce)

{

  double r2inv,r6inv,r,grij,expm2,t,erfc,prefactor,u;

  double fraction,table,forcecoul,forcelj,phicoul,philj;

  int itable;

  double factor_f,factor_e;

  double dqi,dqj,dcoul,asr,exp_asr;

  int di, dj, di_closest;



  int *drudetype = fix_drude->drudetype;

  tagint *drudeid = fix_drude->drudeid;

  int *type = atom->type;



  r2inv = 1.0/rsq;

  if (rsq < cut_coulsq) {

    r = sqrt(rsq);

    if (!ncoultablebits || rsq <= tabinnersq) {

      grij = g_ewald * r;

      expm2 = exp(-grij*grij);

      t = 1.0 / (1.0 + EWALD_P*grij);

      u = 1. - t;

      erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;

      prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;

      forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

      if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;

    } else {

      union_int_float_t rsq_lookup_single;

      rsq_lookup_single.f = rsq;

      itable = rsq_lookup_single.i & ncoulmask;

      itable >>= ncoulshiftbits;

      fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];

      table = ftable[itable] + fraction*dftable[itable];

      forcecoul = atom->q[i]*atom->q[j] * table;

      if (factor_coul < 1.0) {

        table = ctable[itable] + fraction*dctable[itable];

        prefactor = atom->q[i]*atom->q[j] * table;

        forcecoul -= (1.0-factor_coul)*prefactor;

      }

    }

    if (drudetype[type[i]] != NOPOL_TYPE && drudetype[type[j]] != NOPOL_TYPE) {

      di = atom->map(drudeid[i]);

      di_closest = domain->closest_image(i, di);

      if (j != di_closest){

        if (drudetype[i] == CORE_TYPE) dqi = -atom->q[di];

        else if (drudetype[i] == DRUDE_TYPE) dqi = atom->q[i];

        else dqi = 0.0;

        if (drudetype[j] == CORE_TYPE) {

          dj = atom->map(drudeid[j]);

          dqj = -atom->q[dj];

        } else if (drudetype[j] == DRUDE_TYPE) dqj = atom->q[j];

        else dqj = 0.0;

        asr = ascreen[itype][jtype] * r;

        exp_asr = exp(-asr);

        dcoul = force->qqrd2e * dqi * dqj / r;

        factor_f = 0.5*(2. + (exp_asr * (-2. - asr * (2. + asr))))

            - factor_coul;

        forcecoul += factor_f * dcoul;

        factor_e = 0.5*(2. - (exp_asr * (2. + asr))) - factor_coul;

      }

    }

  } else forcecoul = 0.0;



  if (rsq < cut_ljsq[itype][jtype]) {

    r6inv = r2inv*r2inv*r2inv;

    forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);

  } else forcelj = 0.0;



  fforce = (forcecoul + factor_lj*forcelj) * r2inv;



  double eng = 0.0;

  if (rsq < cut_coulsq) {

    if (!ncoultablebits || rsq <= tabinnersq)

      phicoul = prefactor*erfc;

    else {

      table = etable[itable] + fraction*detable[itable];

      phicoul = atom->q[i]*atom->q[j] * table;

    }

    if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;

    if (drudetype[type[i]] != NOPOL_TYPE && drudetype[type[j]] != NOPOL_TYPE &&

        di_closest != j)

      phicoul += factor_e * dcoul;

    eng += phicoul;

  }



  if (rsq < cut_ljsq[itype][jtype]) {

    philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -

      offset[itype][jtype];

    eng += factor_lj*philj;

  }



  return eng;

}

// No point in having single() since it  has no information about topology or Drude particles.

// Charges qi and qj are defined by the user (or 1.0 by defaut)



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

src/USER-DRUDE/pair_lj_cut_thole_long.h

+0 −2
Original line number Diff line number Diff line
@@ -40,8 +40,6 @@ class PairLJCutTholeLong : public Pair { 
  virtual void read_restart_settings(FILE *);

  void write_data(FILE *);

  void write_data_all(FILE *);

  virtual double single(int, int, int, int, double, double, double, double &);



  virtual void *extract(const char *, int &);



 protected:

src/USER-DRUDE/pair_thole.cpp

+6 −24
Original line number Diff line number Diff line
@@ -367,30 +367,12 @@ double PairThole::single(int i, int j, int itype, int jtype, 
{

  double r2inv,rinv,r,phicoul;

  double qi,qj,factor_f,factor_e,dcoul,asr,exp_asr;

  int di, dj;

   

  int *drudetype = fix_drude->drudetype;

  tagint *drudeid = fix_drude->drudeid;

  int *type = atom->type;

  // single() has no information about topology or Drude particles.

  // Charges qi and qj are defined by the user (or 1.0 by defaut)

  

  // only on core-drude pair, but not on the same pair

  if (drudetype[type[i]] == NOPOL_TYPE || drudetype[type[j]] == NOPOL_TYPE ||

      j == i)

    return 0.0;



  // get dq of the core via the drude charge

  if (drudetype[type[i]] == DRUDE_TYPE)

  qi = atom->q[i];

  else {

    di = domain->closest_image(i, atom->map(drudeid[i]));

    qi = -atom->q[di];

  }

  if (drudetype[type[j]] == DRUDE_TYPE)

  qj = atom->q[j];

  else {

    dj = domain->closest_image(j, atom->map(drudeid[j]));

    qj = -atom->q[dj];

  }



  r2inv = 1.0/rsq;

  fforce = phicoul = 0.0;

src/USER-OMP/pair_lj_cut_thole_long_omp.cpp

+8 −2
Original line number Diff line number Diff line
@@ -45,6 +45,10 @@ using namespace MathConst; 
#define B4       -5.80844129e-3

#define B5        1.14652755e-1



#define EPSILON 1.0e-20

#define EPS_EWALD 1.0e-6

#define EPS_EWALD_SQR 1.0e-12



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



PairLJCutTholeLongOMP::PairLJCutTholeLongOMP(LAMMPS *lmp) :
@@ -168,6 +172,7 @@ void PairLJCutTholeLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
      jtype = type[j];



      if (rsq < cutsqi[jtype]) {

        rsq += EPSILON; // Add Epsilon for case: r = 0; DC-DP 1-1 interaction must be removed by special bond;

        r2inv = 1.0/rsq;



        if (rsq < cut_coulsq) {
@@ -175,14 +180,15 @@ void PairLJCutTholeLongOMP::eval(int iifrom, int iito, ThrData * const thr) 
          r = sqrt(rsq);



          if (!ncoultablebits || rsq <= tabinnersq) {

            grij = g_ewald * r;

            grij = g_ewald * (r + EPS_EWALD);

            expm2 = exp(-grij*grij);

            t = 1.0 / (1.0 + EWALD_P*grij);

            u = 1. - t;

            erfc = t * (1.+u*(B0+u*(B1+u*(B2+u*(B3+u*(B4+u*B5)))))) * expm2;

            prefactor = qqrd2e * qi*qj/r;

            prefactor = qqrd2e * qi*qj/(r + EPS_EWALD);

            forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);

            if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;

            r2inv = 1.0/(rsq + EPS_EWALD_SQR);

          } else {

            union_int_float_t rsq_lookup;

            rsq_lookup.f = rsq;

src/fix_propel_self.cpp

0 → 100644
+268 −0
Original line number Diff line number Diff line 
/* ----------------------------------------------------------------------

   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.

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



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

   Contributed by Stefan Paquay @ Brandeis University



   Thanks to Liesbeth Janssen @ Eindhoven University for useful discussions!

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



#include "fix_propel_self.h"



#include <cstdio>

#include <cstring>

#include <string>



#include "atom.h"

#include "atom_vec_ellipsoid.h"

#include "citeme.h"

#include "comm.h"

#include "error.h"

#include "force.h"

#include "group.h"

#include "math.h"

#include "math_const.h"

#include "math_extra.h"

#include "math_vector.h"

#include "modify.h"

#include "random_mars.h"

#include "respa.h"

#include "update.h"



using namespace LAMMPS_NS;

using namespace FixConst;

using namespace MathConst;



#define PRINT_DEBUG_OUTPUT 0



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



FixPropelSelf::FixPropelSelf( LAMMPS *lmp, int narg, char **argv )

  : Fix(lmp, narg, argv), magnitude(0.0),

    mode(VELOCITY), n_types_filter(0), apply_to_type(NULL)

{

  if (narg < 5) error->all(FLERR, "Illegal fix propel/self command");



  // The fix is to support the following cases:

  // 1. Simple atoms, in which case the force points along the velocity

  // 2. Aspherical particles with an orientation.

  // The first argument (mode) is used to differentiate between these.



  // args: fix ID all propel/self mode magnitude

  // Optional args are



  const char *mode_str = argv[3];



  if (strcmp(mode_str, "velocity") == 0) {

    mode = VELOCITY;



  } else if (strcmp(mode_str, "quat") == 0) {



    // This mode should only be supported if the atom style has

    // a quaternion (and if all atoms in the group have it)



    if (!atoms_have_quaternion()) {

      error->all(FLERR, "All fix atoms need to be extended particles");

    }

    mode = QUATERNION;



  } else {

    char msg[2048];

    sprintf(msg, "Illegal mode \"%s\" for fix propel/self", mode_str);

    error->all(FLERR, msg);

  }



  magnitude = force->numeric( FLERR, argv[4] );



  // Handle rest of args:



  int iarg = 5;

  while (iarg < narg) {



    if (strcmp(argv[iarg],"types") == 0) {



      apply_to_type = new int[atom->ntypes+1];

      memset(apply_to_type, 0, atom->ntypes * sizeof(int));



      // consume all following numerical arguments as types



      iarg++;

      int flag=0;

      while (iarg < narg) {

        if (isdigit(argv[iarg][0])) {

          int thistype = force->inumeric(FLERR,argv[iarg]);

          if ((thistype < 1) || (thistype > atom->ntypes))

            error->all(FLERR,"Illegal atom type to types keyword");

          apply_to_type[thistype] = 1;

          flag = 1;

          iarg++;

        } else break;

      }

      if (!flag)

        error->all(FLERR,"'types' keyword requires at least one type");

      else

        n_types_filter = 1;



    } else {

      error->all(FLERR,"Illegal fix propel/self command.");

    }

  }

}



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



FixPropelSelf::~FixPropelSelf()

{

  delete[] apply_to_type;

}

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



int FixPropelSelf::setmask()

{

  int mask = 0;

  mask |= POST_FORCE;



  return mask;

}



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



double FixPropelSelf::memory_usage()

{

  // magnitude + thermostat_orient + mode + n_types_filter + apply_to_type

  double bytes = sizeof(double) + 3*sizeof(int) + sizeof(int*);

  bytes += sizeof(int)*atom->ntypes*n_types_filter;



  return bytes;

}



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



void FixPropelSelf::post_force(int vflag )

{

  switch(mode) {

  case QUATERNION:

    if (n_types_filter) post_force_quaternion<1>(vflag);

    else                post_force_quaternion<0>(vflag);

    break;

  case VELOCITY:

    if (n_types_filter) post_force_velocity<1>(vflag);

    else                post_force_velocity<0>(vflag);

    break;

  default:

    ;

  }

}



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



template <int filter_by_type>

void FixPropelSelf::post_force_quaternion(int /* vflag */ )

{

  double **f = atom->f;

  AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->avec);



  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  int *type = atom->type;

  int* ellipsoid = atom->ellipsoid;



  AtomVecEllipsoid::Bonus *bonus = av->bonus;



  // Add the active force to the atom force:



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

    if( mask[i] & groupbit ){

      if (filter_by_type && !apply_to_type[type[i]]) {

        continue;

      }



      double f_act[3] = { 1.0, 0.0, 0.0 };

      double f_rot[3];



      double *quat  = bonus[ellipsoid[i]].quat;



      double Q[3][3];

      MathExtra::quat_to_mat( quat, Q );

      MathExtra::matvec( Q, f_act, f_rot );



      f[i][0] += magnitude * f_rot[0];

      f[i][1] += magnitude * f_rot[1];

      f[i][2] += magnitude * f_rot[2];

    }

  }

}



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



template <int filter_by_type>

void FixPropelSelf::post_force_velocity(int /*vflag*/)

{

  double **f = atom->f;

  double **v = atom->v;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  int *type = atom->type;



  // Add the active force to the atom force:



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

    if( mask[i] & groupbit ){

      if (filter_by_type && !apply_to_type[type[i]]) {

        continue;

      }



      const double *vi = v[i];

      double f_act[3] = { vi[0], vi[1], vi[2] };

      double nv2 = vi[0]*vi[0] + vi[1]*vi[1] + vi[2]*vi[2];

      double fnorm = 0.0;

      const double TOL = 1e-14;



      if (nv2 > TOL) {



        // Without this check you can run into numerical

        // issues because fnorm will blow up.



        fnorm = magnitude / sqrt(nv2);

      }



      f[i][0] += fnorm * f_act[0];

      f[i][1] += fnorm * f_act[1];

      f[i][2] += fnorm * f_act[2];

    }

  }

}



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



int FixPropelSelf::atoms_have_quaternion()

{

  if (!atom->ellipsoid_flag) {

    error->all(FLERR, "Mode 'quat' requires atom style ellipsoid");

    return 0;

  }



  int *mask = atom->mask;

  int flag=0,flagall=0;



  // Make sure all atoms have ellipsoid data:



  for (int i = 0; i < atom->nlocal; ++i)

    if (mask[i] & groupbit)

      if (atom->ellipsoid[i] < 0) ++flag;



  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);

  if (flagall > 0) return 0;



  return 1;

}

CRA Git | Maintained and supported by SUSTech CRA and CCSE