Commit 67d45d04 authored by sjplimp's avatar sjplimp
Browse files 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@945 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent 15abc50f

src/USER-EWALDN/pair_buck_coul.cpp

+93 −46
Original line number Diff line number Diff line
@@ -238,44 +238,6 @@ void PairBuckCoul::coeff(int narg, char **arg) 
  if (count == 0) error->all("Incorrect args for pair coefficients");

}



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

   init for one type pair i,j and corresponding j,i

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



double PairBuckCoul::init_one(int i, int j)

{

  if (setflag[i][j] == 0) error->all("All pair coeffs are not set");



  cut_buck[i][j] = cut_buck_read[i][j];

  buck_a[i][j] = buck_a_read[i][j];

  buck_c[i][j] = buck_c_read[i][j];

  buck_rho[i][j] = buck_rho_read[i][j];



  double cut = MAX(cut_buck[i][j], cut_coul);

  cutsq[i][j] = cut*cut;

  cut_bucksq[i][j] = cut_buck[i][j] * cut_buck[i][j];



  buck1[i][j] = buck_a[i][j]/buck_rho[i][j];

  buck2[i][j] = 6.0*buck_c[i][j];

  rhoinv[i][j] = 1.0/buck_rho[i][j];

     

  if (offset_flag) {

    double rexp = exp(-cut_buck[i][j]/buck_rho[i][j]);

    offset[i][j] = buck_a[i][j]*rexp - buck_c[i][j]/pow(cut_buck[i][j],6.0);

  } else offset[i][j] = 0.0;



  cutsq[j][i] = cutsq[i][j];

  cut_bucksq[j][i] = cut_bucksq[i][j];

  buck_a[j][i] = buck_a[i][j];

  buck_c[j][i] = buck_c[i][j];

  rhoinv[j][i] = rhoinv[i][j];

  buck1[j][i] = buck1[i][j];

  buck2[j][i] = buck2[i][j];

  offset[j][i] = offset[i][j];



  return cut;

}



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

   init specific to this pair style

------------------------------------------------------------------------- */
@@ -286,12 +248,46 @@ void PairBuckCoul::init_style() 


  // require an atom style with charge defined



  //if (atom->charge_allow == 0)

    //error->all("Must use charged atom style with this pair style");

  if (!atom->q_flag && (ewald_order&(1<<1)))

    error->all(

	"Invoking coulombic in pair style lj/coul requires atom attribute q");



  // request regular or rRESPA neighbor lists



  int irequest;



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

    int respa = 0;

    if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;

    if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;



    if (respa == 0) irequest = neighbor->request(this);

    else if (respa == 1) {

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 1;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respainner = 1;

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 3;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respaouter = 1;

    } else {

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 1;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respainner = 1;

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 2;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respamiddle = 1;

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 3;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respaouter = 1;

    }



  } else irequest = neighbor->request(this);



  cut_coulsq = cut_coul * cut_coul;



  // set & error check interior rRESPA cutoffs
@@ -330,6 +326,57 @@ void PairBuckCoul::init_style() 
  if (ncoultablebits) init_tables();

}



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

   neighbor callback to inform pair style of neighbor list to use

   regular or rRESPA

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



void PairBuckCoul::init_list(int id, NeighList *ptr)

{

  if (id == 0) list = ptr;

  else if (id == 1) listinner = ptr;

  else if (id == 2) listmiddle = ptr;

  else if (id == 3) listouter = ptr;

}



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

   init for one type pair i,j and corresponding j,i

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



double PairBuckCoul::init_one(int i, int j)

{

  if (setflag[i][j] == 0) error->all("All pair coeffs are not set");



  cut_buck[i][j] = cut_buck_read[i][j];

  buck_a[i][j] = buck_a_read[i][j];

  buck_c[i][j] = buck_c_read[i][j];

  buck_rho[i][j] = buck_rho_read[i][j];



  double cut = MAX(cut_buck[i][j], cut_coul);

  cutsq[i][j] = cut*cut;

  cut_bucksq[i][j] = cut_buck[i][j] * cut_buck[i][j];



  buck1[i][j] = buck_a[i][j]/buck_rho[i][j];

  buck2[i][j] = 6.0*buck_c[i][j];

  rhoinv[i][j] = 1.0/buck_rho[i][j];

     

  if (offset_flag) {

    double rexp = exp(-cut_buck[i][j]/buck_rho[i][j]);

    offset[i][j] = buck_a[i][j]*rexp - buck_c[i][j]/pow(cut_buck[i][j],6.0);

  } else offset[i][j] = 0.0;



  cutsq[j][i] = cutsq[i][j];

  cut_bucksq[j][i] = cut_bucksq[i][j];

  buck_a[j][i] = buck_a[i][j];

  buck_c[j][i] = buck_c[i][j];

  rhoinv[j][i] = rhoinv[i][j];

  buck1[j][i] = buck1[i][j];

  buck2[j][i] = buck2[i][j];

  offset[j][i] = offset[i][j];



  return cut;

}



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

  proc 0 writes to restart file

------------------------------------------------------------------------- */
@@ -601,7 +648,7 @@ void PairBuckCoul::compute_inner() 
  double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi, *offseti;

  vector xi, d;



  ineighn = (ineigh = list->ilist)+list->inum;

  ineighn = (ineigh = listinner->ilist)+listinner->inum;



  for (; ineigh<ineighn; ++ineigh) {			// loop over my atoms

    i = *ineigh; fi = f0+3*i;
@@ -610,7 +657,7 @@ void PairBuckCoul::compute_inner() 
    offseti = offset[typei = type[i]];

    cut_bucksqi = cut_bucksq[typei];

    buck1i = buck1[typei]; buck2i = buck2[typei]; rhoinvi = rhoinv[typei];

    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];

    jneighn = (jneigh = listinner->firstneigh[i])+listinner->numneigh[i];



    for (; jneigh<jneighn; ++jneigh) {			// loop over neighbors

      if ((j = *jneigh) < nall) ni = -1;
@@ -692,7 +739,7 @@ void PairBuckCoul::compute_middle() 
  double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi, *offseti;

  vector xi, d;



  ineighn = (ineigh = list->ilist)+list->inum;

  ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum;



  for (; ineigh<ineighn; ++ineigh) {			// loop over my atoms

    i = *ineigh; fi = f0+3*i;
@@ -701,7 +748,7 @@ void PairBuckCoul::compute_middle() 
    offseti = offset[typei = type[i]];

    cut_bucksqi = cut_bucksq[typei];

    buck1i = buck1[typei]; buck2i = buck2[typei]; rhoinvi = rhoinv[typei];

    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];

    jneighn = (jneigh = listmiddle->firstneigh[i])+listmiddle->numneigh[i];



    for (; jneigh<jneighn; ++jneigh) {			// loop over neighbors

      if ((j = *jneigh) < nall) ni = -1;
@@ -790,7 +837,7 @@ void PairBuckCoul::compute_outer(int eflag, int vflag) 
  eng_vdwl = eng_coul = 0.0;				// reset energy&virial

  if (vflag) memset(virial, 0, sizeof(shape));



  ineighn = (ineigh = list->ilist)+list->inum;

  ineighn = (ineigh = listouter->ilist)+listouter->inum;



  for (; ineigh<ineighn; ++ineigh) {			// loop over my atoms

    i = *ineigh; fi = f0+3*i;
@@ -800,7 +847,7 @@ void PairBuckCoul::compute_outer(int eflag, int vflag) 
    buckai = buck_a[typei]; buckci = buck_c[typei]; rhoinvi = rhoinv[typei];

    cutsqi = cutsq[typei]; cut_bucksqi = cut_bucksq[typei];

    memcpy(xi, x0+(i+(i<<1)), sizeof(vector));

    jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i];

    jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i];



    for (; jneigh<jneighn; ++jneigh) {			// loop over neighbors

      if ((j = *jneigh) < nall) ni = -1;

src/USER-EWALDN/pair_buck_coul.h

+7 −6
Original line number Diff line number Diff line
@@ -28,16 +28,17 @@ class PairBuckCoul : public Pair { 


  virtual void settings(int, char **);

  void coeff(int, char **);

  void init_style();

  void init_list(int, class NeighList *);

  double init_one(int, int);

  virtual void init_style();

  void write_restart(FILE *);

  void read_restart(FILE *);

  

  virtual void write_restart_settings(FILE *);

  virtual void read_restart_settings(FILE *);

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

  virtual void *extract_ptr(char *);

  virtual void extract_long(double *);

  void write_restart_settings(FILE *);

  void read_restart_settings(FILE *);

  void single(int, int, int, int, double, double, double, int, One &);

  void *extract_ptr(char *);

  void extract_long(double *);



  void compute_inner();

  void compute_middle();

src/USER-EWALDN/pair_lj_coul.cpp

+95 −48
Original line number Diff line number Diff line
@@ -233,52 +233,6 @@ void PairLJCoul::coeff(int narg, char **arg) 
  if (count == 0) error->all("Incorrect args for pair coefficients");

}



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

   init for one type pair i,j and corresponding j,i

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



double PairLJCoul::init_one(int i, int j)

{

  if ((ewald_order&(1<<6))||(setflag[i][j] == 0)) {

    epsilon[i][j] = mix_energy(epsilon_read[i][i],epsilon_read[j][j],

			       sigma_read[i][i],sigma_read[j][j]);

    sigma[i][j] = mix_distance(sigma_read[i][i],sigma_read[j][j]);

    if (ewald_order&(1<<6))

      cut_lj[i][j] = cut_lj_global;

    else

      cut_lj[i][j] = mix_distance(cut_lj_read[i][i],cut_lj_read[j][j]);

  }

  else {

    sigma[i][j] = sigma_read[i][j];

    epsilon[i][j] = epsilon_read[i][j];

    cut_lj[i][j] = cut_lj_read[i][j];

  }



  double cut = MAX(cut_lj[i][j], cut_coul);

  cutsq[i][j] = cut*cut;

  cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];



  lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);

  lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);

  lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);

  lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);

     

  if (offset_flag) {

    double ratio = sigma[i][j] / cut_lj[i][j];

    offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));

  } else offset[i][j] = 0.0;



  cutsq[j][i] = cutsq[i][j];

  cut_ljsq[j][i] = cut_ljsq[i][j];

  lj1[j][i] = lj1[i][j];

  lj2[j][i] = lj2[i][j];

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

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

  offset[j][i] = offset[i][j];



  return cut;

}



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

   init specific to this pair style

------------------------------------------------------------------------- */
@@ -291,12 +245,46 @@ void PairLJCoul::init_style() 


  // require an atom style with charge defined



  //if (atom->charge_allow == 0)

    //error->all("Must use charged atom style with this pair style");

  if (!atom->q_flag && (ewald_order&(1<<1)))

    error->all(

	"Invoking coulombic in pair style lj/coul requires atom attribute q");



  // request regular or rRESPA neighbor lists



  int irequest;



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

    int respa = 0;

    if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;

    if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;



    if (respa == 0) irequest = neighbor->request(this);

    else if (respa == 1) {

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 1;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respainner = 1;

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 3;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respaouter = 1;

    } else {

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 1;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respainner = 1;

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 2;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respamiddle = 1;

      irequest = neighbor->request(this);

      neighbor->requests[irequest]->id = 3;

      neighbor->requests[irequest]->half = 0;

      neighbor->requests[irequest]->respaouter = 1;

    }



  } else irequest = neighbor->request(this);



  cut_coulsq = cut_coul * cut_coul;



  // set & error check interior rRESPA cutoffs
@@ -332,6 +320,65 @@ void PairLJCoul::init_style() 
  if (ncoultablebits) init_tables();

}



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

   neighbor callback to inform pair style of neighbor list to use

   regular or rRESPA

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



void PairLJCoul::init_list(int id, NeighList *ptr)

{

  if (id == 0) list = ptr;

  else if (id == 1) listinner = ptr;

  else if (id == 2) listmiddle = ptr;

  else if (id == 3) listouter = ptr;

}



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

   init for one type pair i,j and corresponding j,i

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



double PairLJCoul::init_one(int i, int j)

{

  if ((ewald_order&(1<<6))||(setflag[i][j] == 0)) {

    epsilon[i][j] = mix_energy(epsilon_read[i][i],epsilon_read[j][j],

			       sigma_read[i][i],sigma_read[j][j]);

    sigma[i][j] = mix_distance(sigma_read[i][i],sigma_read[j][j]);

    if (ewald_order&(1<<6))

      cut_lj[i][j] = cut_lj_global;

    else

      cut_lj[i][j] = mix_distance(cut_lj_read[i][i],cut_lj_read[j][j]);

  }

  else {

    sigma[i][j] = sigma_read[i][j];

    epsilon[i][j] = epsilon_read[i][j];

    cut_lj[i][j] = cut_lj_read[i][j];

  }



  double cut = MAX(cut_lj[i][j], cut_coul);

  cutsq[i][j] = cut*cut;

  cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j];



  lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0);

  lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0);

  lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0);

  lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0);

     

  if (offset_flag) {

    double ratio = sigma[i][j] / cut_lj[i][j];

    offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0));

  } else offset[i][j] = 0.0;



  cutsq[j][i] = cutsq[i][j];

  cut_ljsq[j][i] = cut_ljsq[i][j];

  lj1[j][i] = lj1[i][j];

  lj2[j][i] = lj2[i][j];

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

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

  offset[j][i] = offset[i][j];



  return cut;

}



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

  proc 0 writes to restart file

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

src/USER-EWALDN/pair_lj_coul.h

+7 −6
Original line number Diff line number Diff line
@@ -27,16 +27,17 @@ class PairLJCoul : public Pair { 
  virtual void compute(int, int);

  virtual void settings(int, char **);

  void coeff(int, char **);

  void init_style();

  void init_list(int, class NeighList *);

  double init_one(int, int);

  virtual void init_style();

  void write_restart(FILE *);

  void read_restart(FILE *);

  

  virtual void write_restart_settings(FILE *);

  virtual void read_restart_settings(FILE *);

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

  virtual void *extract_ptr(char *);

  virtual void extract_long(double *);

  void write_restart_settings(FILE *);

  void read_restart_settings(FILE *);

  void single(int, int, int, int, double, double, double, int, One &);

  void *extract_ptr(char *);

  void extract_long(double *);



  void compute_inner();

  void compute_middle();

CRA Git | Maintained and supported by SUSTech CRA and CCSE