Commit 8b1ef1c6 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

import updated pair lj/sf/dipole/sf files from sam genheden 

(cherry picked from commit ae691ab7865b1173d386e003e48df4559291c27e)
parent c3e8cb2f

src/USER-MISC/pair_lj_sf_dipole_sf.cpp

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



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

   Contributing author: Mario Orsi (U Southampton), orsimario@gmail.com

   Contributing authors: Mario Orsi (QMUL), m.orsi@qmul.ac.uk

                         Samuel Genheden (University of Southampton)

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



#include <math.h>
@@ -34,7 +35,6 @@ using namespace LAMMPS_NS; 


PairLJSFDipoleSF::PairLJSFDipoleSF(LAMMPS *lmp) : Pair(lmp)

{

  single_enable = 0;

}



/* ---------------------------------------------------------------------- */
@@ -55,6 +55,7 @@ PairLJSFDipoleSF::~PairLJSFDipoleSF() 
    memory->destroy(lj2);

    memory->destroy(lj3);

    memory->destroy(lj4);

    memory->destroy(scale);

  }

}
@@ -86,6 +87,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) 
  double **torque = atom->torque;

  int *type = atom->type;

  int nlocal = atom->nlocal;

// The global scaling parameters aren't used anymore

  double *special_coul = force->special_coul;

  double *special_lj = force->special_lj;

  int newton_pair = force->newton_pair;
@@ -234,7 +236,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) 


        // total force



        fq = factor_coul*qqrd2e;

        fq = factor_coul*qqrd2e*scale[itype][jtype];

        fx = fq*forcecoulx + delx*forcelj;

        fy = fq*forcecouly + dely*forcelj;

        fz = fq*forcecoulz + delz*forcelj;
@@ -268,7 +270,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag) 
              ecoul += -q[j] * r3inv * pqfac * pidotr;

            if (mu[j][3] > 0.0 && qtmp != 0.0)

              ecoul += qtmp * r3inv * qpfac * pjdotr;

            ecoul *= factor_coul*qqrd2e;

            ecoul *= factor_coul*qqrd2e*scale[itype][jtype];

          } else ecoul = 0.0;



          if (rsq < cut_ljsq[itype][jtype]) {
@@ -315,6 +317,7 @@ void PairLJSFDipoleSF::allocate() 
  memory->create(lj2,n+1,n+1,"pair:lj2");

  memory->create(lj3,n+1,n+1,"pair:lj3");

  memory->create(lj4,n+1,n+1,"pair:lj4");

  memory->create(scale,n+1,n+1,"pair:scale");

}



/* ----------------------------------------------------------------------
@@ -352,7 +355,7 @@ void PairLJSFDipoleSF::settings(int narg, char **arg) 


void PairLJSFDipoleSF::coeff(int narg, char **arg)

{

  if (narg < 4 || narg > 6)

  if (narg < 4 || narg > 7)

    error->all(FLERR,"Incorrect args for pair coefficients");

  if (!allocated) allocate();
@@ -363,10 +366,13 @@ void PairLJSFDipoleSF::coeff(int narg, char **arg) 
  double epsilon_one = force->numeric(FLERR,arg[2]);

  double sigma_one = force->numeric(FLERR,arg[3]);



  double scale_one = 1.0;

  if (narg >= 5) scale_one = force->numeric(FLERR,arg[4]);



  double cut_lj_one = cut_lj_global;

  double cut_coul_one = cut_coul_global;

  if (narg >= 5) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[4]);

  if (narg == 6) cut_coul_one = force->numeric(FLERR,arg[5]);

  if (narg >= 6) cut_coul_one = cut_lj_one = force->numeric(FLERR,arg[5]);

  if (narg == 7) cut_coul_one = force->numeric(FLERR,arg[6]);



  int count = 0;

  for (int i = ilo; i <= ihi; i++) {
@@ -376,6 +382,7 @@ void PairLJSFDipoleSF::coeff(int narg, char **arg) 
      cut_lj[i][j] = cut_lj_one;

      cut_coul[i][j] = cut_coul_one;

      setflag[i][j] = 1;

      scale[i][j] = scale_one;

      count++;

    }

  }
@@ -424,6 +431,7 @@ double PairLJSFDipoleSF::init_one(int i, int j) 
  lj2[j][i] = lj2[i][j];

  lj3[j][i] = lj3[i][j];

  lj4[j][i] = lj4[i][j];

  scale[j][i] = scale[i][j];



  return cut;

}
@@ -445,6 +453,7 @@ void PairLJSFDipoleSF::write_restart(FILE *fp) 
        fwrite(&sigma[i][j],sizeof(double),1,fp);

        fwrite(&cut_lj[i][j],sizeof(double),1,fp);

        fwrite(&cut_coul[i][j],sizeof(double),1,fp);

        fwrite(&scale[i][j],sizeof(double),1,fp);

      }

    }

}
@@ -471,11 +480,13 @@ void PairLJSFDipoleSF::read_restart(FILE *fp) 
          fread(&sigma[i][j],sizeof(double),1,fp);

          fread(&cut_lj[i][j],sizeof(double),1,fp);

          fread(&cut_coul[i][j],sizeof(double),1,fp);

          fread(&scale[i][j],sizeof(double),1,fp);

        }

        MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);

        MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);

        MPI_Bcast(&cut_lj[i][j],1,MPI_DOUBLE,0,world);

        MPI_Bcast(&cut_coul[i][j],1,MPI_DOUBLE,0,world);

        MPI_Bcast(&scale[i][j],1,MPI_DOUBLE,0,world);

      }

    }

}
@@ -506,3 +517,105 @@ void PairLJSFDipoleSF::read_restart_settings(FILE *fp) 
  MPI_Bcast(&cut_coul_global,1,MPI_DOUBLE,0,world);

  MPI_Bcast(&mix_flag,1,MPI_INT,0,world);

}



// PairLJSFDipoleSF: calculation of force is missing (to be implemented)

double PairLJSFDipoleSF::single(int i, int j, int itype, int jtype, double rsq,

				double factor_coul, double factor_lj,

				double &fforce)

{

  double r2inv,r6inv;

  double pdotp,pidotr,pjdotr,pre1,delx,dely,delz;

  double rinv, r3inv,r5inv, rcutlj2inv, rcutcoul2inv,rcutlj6inv;

  double qtmp,xtmp,ytmp,ztmp,bfac,pqfac,qpfac, ecoul, evdwl;



  double **x = atom->x;

  double *q = atom->q;

  double **mu = atom->mu;



  qtmp = q[i];

  xtmp = x[i][0];

  ytmp = x[i][1];

  ztmp = x[i][2];



  r2inv = 1.0/rsq;

  rinv = sqrt(r2inv);

  fforce = 0.0;



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

    delx = xtmp - x[j][0];

    dely = ytmp - x[j][1];

    delz = ztmp - x[j][2];

    // if (qtmp != 0.0 && q[j] != 0.0) {

    //   pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]);

    //   forcecoulx += pre1*delx;

    //   forcecouly += pre1*dely;

    //   forcecoulz += pre1*delz;

    // }

    if (mu[i][3] > 0.0 && mu[j][3] > 0.0) {

      r3inv = r2inv*rinv;

      r5inv = r3inv*r2inv;

      rcutcoul2inv=1.0/cut_coulsq[itype][jtype];

      pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];

      pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;

      pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;

      bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) +

	3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv;

    }

    if (mu[i][3] > 0.0 && q[j] != 0.0) {

      r3inv = r2inv*rinv;

      r5inv = r3inv*r2inv;

      pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;

      rcutcoul2inv=1.0/cut_coulsq[itype][jtype];

      pqfac = 1.0 - 3.0*rsq*rcutcoul2inv +

	2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);

    }

    if (mu[j][3] > 0.0 && qtmp != 0.0) {

      r3inv = r2inv*rinv;

      r5inv = r3inv*r2inv;

      pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;

      rcutcoul2inv=1.0/cut_coulsq[itype][jtype];

      qpfac = 1.0 - 3.0*rsq*rcutcoul2inv +

	2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);

    }

  }

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

    r6inv = r2inv*r2inv*r2inv;

    rcutlj2inv = 1.0 / cut_ljsq[itype][jtype];

    rcutlj6inv = rcutlj2inv * rcutlj2inv * rcutlj2inv;

  }



  double eng = 0.0;

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

    ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype]));

    ecoul *= ecoul;

    ecoul *= qtmp * q[j] * rinv;

    if (mu[i][3] > 0.0 && mu[j][3] > 0.0)

      ecoul += bfac * (r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr);

    if (mu[i][3] > 0.0 && q[j] != 0.0)

      ecoul += -q[j] * r3inv * pqfac * pidotr;

    if (mu[j][3] > 0.0 && qtmp != 0.0)

      ecoul += qtmp * r3inv * qpfac * pjdotr;

    ecoul *= factor_coul*force->qqrd2e*scale[itype][jtype];

    eng += ecoul;

  }

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

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

      rcutlj6inv*(6*lj3[itype][jtype]*rcutlj6inv-3*lj4[itype][jtype])*

      rsq*rcutlj2inv+

      rcutlj6inv*(-7*lj3[itype][jtype]*rcutlj6inv+4*lj4[itype][jtype]);

    eng += evdwl*factor_lj;

  }



  return eng;

}



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



void *PairLJSFDipoleSF::extract(const char *str, int &dim)

{

  dim = 2;

  if (strcmp(str,"epsilon") == 0) return (void *) epsilon;

  if (strcmp(str,"sigma") == 0) return (void *) sigma;

  if (strcmp(str,"scale") == 0) return (void *) scale;

  return NULL;

}

src/USER-MISC/pair_lj_sf_dipole_sf.h

+4 −1
Original line number Diff line number Diff line 
/* -*- c++ -*- ----------------------------------------------------------

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

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov
@@ -37,6 +37,8 @@ class PairLJSFDipoleSF : public Pair { 
  void read_restart(FILE *);

  void write_restart_settings(FILE *);

  void read_restart_settings(FILE *);

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

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



 protected:

  double cut_lj_global,cut_coul_global;
@@ -44,6 +46,7 @@ class PairLJSFDipoleSF : public Pair { 
  double **cut_coul,**cut_coulsq;

  double **epsilon,**sigma;

  double **lj1,**lj2,**lj3,**lj4;

  double **scale;



  void allocate();

};

CRA Git | Maintained and supported by SUSTech CRA and CCSE