Merge pull request #169 from ludwig-cf/fix-wall-colloid

Re-instate wall colloid potential

CHANGES.md

@@ -1,6 +1,12 @@
 
 ### Changes
 
+version 0.16.0
+
+- Reinstated the boundary (wall) - colloid soft sphere potential.
+  See https://ludwig.epcc.ed.ac.uk/inputs/colloid.html
+  Thanks to Rishish Mishra for spotting this problem.
+
 version 0.15.0
 
 - Active stress implementation is updated to conform to the documented

src/colloids_rt.c

@@ -7,7 +7,7 @@
  *  Edinburgh Soft Matter and Statistical Physics Group and
  *  Edinburgh Parallel Computing Centre
  *
- *  (c) 2014-2021 The University of Edinburgh
+ *  (c) 2014-2022 The University of Edinburgh
  *
  *  Contributing authors:
  *  Kevin Stratford ([email protected])
@@ -32,6 +32,7 @@
 #include "pair_yukawa.h"
 #include "bond_fene.h"
 #include "angle_cosine.h"
+#include "wall_ss_cut.h"
 
 #include "colloids_halo.h"
 #include "colloids_init.h"
@@ -47,9 +48,11 @@ int pair_yukawa_init(pe_t * pe, cs_t * cs, rt_t * rt, interact_t * inter);
 int pair_lj_cut_init(pe_t * pe, cs_t * cs, rt_t * rt, interact_t * inter);
 int bond_fene_init(pe_t * pe, cs_t * cs, rt_t * rt, interact_t * interact);
 int angle_cosine_init(pe_t * pe, cs_t * cs, rt_t * rt, interact_t * interact);
-
 int pair_ss_cut_ij_init(pe_t * pe, cs_t * cs, rt_t * rt, interact_t * intrct);
 
+int wall_ss_cut_init(pe_t * pe, cs_t * cs, rt_t * rt, wall_t * wall,
+		     interact_t * inter);
+
 int colloids_rt_dynamics(cs_t * cs, colloids_info_t * cinfo, wall_t * wall,
 			 map_t * map, const lb_model_t * model);
 int colloids_rt_gravity(pe_t * pe, rt_t * rt, colloids_info_t * cinfo);
@@ -146,6 +149,8 @@ int colloids_init_rt(pe_t * pe, rt_t * rt, cs_t * cs, colloids_info_t ** pinfo,
 
   pair_ss_cut_ij_init(pe, cs, rt, *interact);
 
+  wall_ss_cut_init(pe, cs, rt, wall, *interact);
+
   colloids_rt_cell_list_checks(pe, cs, pinfo, *interact);
   colloids_init_halo_range_check(pe, cs, *pinfo);
   if (nc > 1) interact_range_check(*interact, *pinfo);
@@ -756,18 +761,13 @@ int pair_ss_cut_init(pe_t * pe, cs_t * cs, rt_t * rt, interact_t * inter) {
   double epsilon ;
   double sigma;
   int nu;
-  double kt;
   double cutoff;
 
-  physics_t * phys = NULL;
   pair_ss_cut_t * pair = NULL;
 
   assert(pe);
   assert(rt);
 
-  physics_ref(&phys);
-  physics_kt(phys, &kt);
-
   rt_int_parameter(rt, "soft_sphere_on", &on);
 
   if (on) {
@@ -1071,3 +1071,46 @@ int colloids_init_halo_range_check(pe_t * pe, cs_t * cs,
 
   return ifail;
 }
+
+/*****************************************************************************
+ *
+ *  wall_ss_cut_init
+ *
+ *****************************************************************************/
+
+int wall_ss_cut_init(pe_t * pe, cs_t * cs, rt_t * rt, wall_t * wall,
+		     interact_t * interact) {
+
+  int have_wall_ss_cut = 0;
+
+  assert(pe);
+  assert(cs);
+  assert(rt);
+
+  have_wall_ss_cut = rt_switch(rt, "wall_ss_cut_on");
+
+  if (have_wall_ss_cut) {
+
+    wall_ss_cut_t * wall_ss_cut = NULL;
+    wall_ss_cut_options_t opts = {};
+
+    rt_key_required(rt, "wall_ss_cut_epsilon", RT_FATAL);
+    rt_key_required(rt, "wall_ss_cut_sigma", RT_FATAL);
+    rt_key_required(rt, "wall_ss_cut_nu", RT_FATAL);
+    rt_key_required(rt, "wall_ss_cut_hc", RT_FATAL);
+
+    rt_double_parameter(rt, "wall_ss_cut_epsilon", &opts.epsilon);
+    rt_double_parameter(rt, "wall_ss_cut_sigma", &opts.sigma);
+    rt_double_parameter(rt, "wall_ss_cut_nu", &opts.nu);
+    rt_double_parameter(rt, "wall_ss_cut_hc", &opts.hc);
+
+    if (opts.nu <= 0) pe_fatal(pe, "Please ensure wall_ss_cut_nu is +ve\n");
+    if (opts.hc <= 0) pe_fatal(pe, "Please ensure wall_ss_cut_hc is +ve\n");
+
+    wall_ss_cut_create(pe, cs, wall, &opts, &wall_ss_cut);
+    wall_ss_cut_register(wall_ss_cut, interact);
+    wall_ss_cut_info(wall_ss_cut);
+  }
+
+  return 0;
+}

src/wall_ss_cut.c

@@ -23,7 +23,7 @@
  *  Edinburgh Soft Matter and Statistical Physics Group and
  *  Edinburgh Parallel Computing Centre
  *
- *  (c) 2010-2016 The University of Edinburgh
+ *  (c) 2010-2022 The University of Edinburgh
  *  Juho Lintuvuori ([email protected])
  *
  *****************************************************************************/
@@ -47,7 +47,7 @@ struct wall_ss_cut_s {
   double hc;             /* cut-off */
   double vlocal;         /* local contribution to energy */
   double hminlocal;      /* local nearest separation */
-  double rminlocal;      /* local min centre-centre separation */
+  double rminlocal;      /* local min wall-centre separation */
 };
 
 /*****************************************************************************
@@ -57,12 +57,15 @@ struct wall_ss_cut_s {
  *****************************************************************************/
 
 int wall_ss_cut_create(pe_t * pe, cs_t * cs, wall_t * wall,
+		       const wall_ss_cut_options_t * opts,
 		       wall_ss_cut_t ** pobj) {
 
   wall_ss_cut_t * obj = NULL;
 
   assert(pe);
   assert(cs);
+  assert(wall);
+  assert(opts);
   assert(pobj);
 
   obj = (wall_ss_cut_t *) calloc(1, sizeof(wall_ss_cut_t));
@@ -73,6 +76,11 @@ int wall_ss_cut_create(pe_t * pe, cs_t * cs, wall_t * wall,
   obj->cs = cs;
   obj->wall = wall;
 
+  obj->epsilon = opts->epsilon;
+  obj->sigma = opts->sigma;
+  obj->nu = opts->nu;
+  obj->hc = opts->hc;
+
   *pobj = obj;
 
   return 0;
@@ -93,24 +101,6 @@ int wall_ss_cut_free(wall_ss_cut_t * obj) {
   return 0;
 }
 
-/*****************************************************************************
- *
- *  wall_ss_cut_param_set
- *
- *****************************************************************************/
-
-int wall_ss_cut_param_set(wall_ss_cut_t * obj, double epsilon, double sigma,
-			  int nu, double hc) {
-  assert(obj);
-
-  obj->epsilon = epsilon;
-  obj->sigma = sigma;
-  obj->nu = nu;
-  obj->hc = hc;
-
-  return 0;
-}
-
 /*****************************************************************************
  *
  *  wall_ss_cut_info
@@ -126,11 +116,12 @@ int wall_ss_cut_info(wall_ss_cut_t * obj) {
   physics_kt(phys, &kt);
 
   pe_info(obj->pe, "\n");
-  pe_info(obj->pe, "Soft sphere for wall potential\n");
+  pe_info(obj->pe, "Wall-colloid soft-sphere potential\n");
+  pe_info(obj->pe, "----------------------------------\n");
   pe_info(obj->pe, "epsilon:                  %14.7e\n", obj->epsilon);
   pe_info(obj->pe, "sigma:                    %14.7e\n", obj->sigma);
   pe_info(obj->pe, "exponent nu:              %14.7e\n", obj->nu);
-  pe_info(obj->pe, "cut off (surface-surface) %14.7e\n", obj->hc);
+  pe_info(obj->pe, "cut off hc (wall-surface) %14.7e\n", obj->hc);
   if (kt > 0.0) {
     pe_info(obj->pe, "epsilon / kT              %14.7e\n", obj->epsilon/kt);
   }
@@ -166,23 +157,11 @@ int wall_ss_cut_compute(colloids_info_t * cinfo, void * obj) {
 
   wall_ss_cut_t * self = (wall_ss_cut_t *) obj;
 
-  int ic1, jc1, kc1;
-  int ncell[3];
-  int ia;
-  int iswall[3];
-
-  double r;                             /* centre-centre sepration */
-  double h;                             /* surface-surface separation */
-  double rh;                            /* reciprocal h */
-  double rsigma;                        /* reciproal sigma */
-  double vcut;                          /* potential at cut off */
-  double dvcut;                         /* derivative at cut off */
-  double forcewall[3];                  /* force on the wall for accounting purposes */
+  double forcewall[3] = {};             /* force on the wall for accounting */
   double lmin[3];
   double ltot[3];
-  double f;
 
-  colloid_t * pc1;
+  colloid_t * pc = NULL;
 
   assert(cinfo);
   assert(self);
@@ -194,65 +173,49 @@ int wall_ss_cut_compute(colloids_info_t * cinfo, void * obj) {
   self->hminlocal = dmax(ltot[X], dmax(ltot[Y], ltot[Z]));
   self->rminlocal = self->hminlocal;
 
-  rsigma = 1.0/self->sigma;
-  vcut = self->epsilon*pow(self->sigma/self->hc, self->nu);
-  dvcut = -self->epsilon*self->nu*rsigma*pow(self->sigma/self->hc, self->nu+1);
+  colloids_info_local_head(cinfo, &pc);
 
-  forcewall[X] = 0.0;
-  forcewall[Y] = 0.0;
-  forcewall[Z] = 0.0;
+  for (; pc; pc = pc->nextlocal) {
 
-  wall_present_dim(self->wall, iswall);
-  colloids_info_ncell(cinfo, ncell);
+    for (int ia = 0; ia < 3; ia++) {
 
-  for (ic1 = 1; ic1 <= ncell[X]; ic1++) {
-    for (jc1 = 1; jc1 <= ncell[Y]; jc1++) {
-      for (kc1 = 1; kc1 <= ncell[Z]; kc1++) {
-
-        colloids_info_cell_list_head(cinfo, ic1, jc1, kc1, &pc1);
-        for (; pc1; pc1 = pc1->next) {
+      double f = 0.0;
+      double r = 0.0;                       /* wall-centre sepration */
+      double h = 0.0;                       /* wall-surface separation */
 
-	  for (ia = 0; ia < 3; ia++) {
-	    if (iswall[ia]) {
-
-	      f = 0.0;
-	      /* lower wall */
-	      r = pc1->s.r[ia] - lmin[ia];
-	      if (r < self->rminlocal) self->rminlocal = r;
+      if (self->wall->param->isboundary[ia] == 0) continue;
 
-	      h = r - pc1->s.ah;
-	      assert(h > 0.0);
-	      if (h < self->hminlocal) self->hminlocal = h;
-
-	      if (h < self->hc) {
-		rh = 1.0/h;
-		self->vlocal += self->epsilon*pow(rh*self->sigma, self->nu)
-		  - vcut - (h - self->hc)*dvcut;
-		f = -(-self->epsilon*self->nu*rsigma
-		      *pow(rh*self->sigma, self->nu+1) - dvcut);
-	      }
+      /* lower wall */
+      r = pc->s.r[ia] - lmin[ia];
+      h = r - pc->s.ah;
+
+      if (h < self->hminlocal) self->hminlocal = h;
+      if (r < self->rminlocal) self->rminlocal = r;
 
-	      /*upper wall*/
-	      r = lmin[ia] + ltot[ia] - pc1->s.r[ia];
-	      if (r < self->rminlocal) self->rminlocal = r;
+      if (h < self->hc) {
+	double v = 0.0;
+	double fl = 0.0;
+	wall_ss_cut_single(self, h, &fl, &v);
+	self->vlocal += v;
+	f += fl;
+      }
 
-	      h = r - pc1->s.ah;
-	      assert(h > 0.0);
-	      if (h < self->hminlocal) self->hminlocal = h;
+      /* upper wall */
+      r = lmin[ia] + ltot[ia] - pc->s.r[ia];
+      h = r - pc->s.ah;
+
+      if (r < self->rminlocal) self->rminlocal = r;
+      if (h < self->hminlocal) self->hminlocal = h;
 
-	      if (h < self->hc) {
-		rh = 1.0/h;
-		self->vlocal += self->epsilon*pow(rh*self->sigma, self->nu)
-		  - vcut - (h - self->hc)*dvcut;
-		f -= -(-self->epsilon*self->nu*rsigma
-		       *pow(rh*self->sigma, self->nu+1) - dvcut);
-	      }
-	      pc1->force[ia] += f;
-	      forcewall[ia] -= f;
-	    }
-	  }
-	}
+      if (h < self->hc) {
+	double v = 0.0;
+	double fu = 0.0;
+	wall_ss_cut_single(self, h, &fu, &v);
+	self->vlocal += v;
+	f -= fu;            /* upper wall gives -ve (repulsive) force */
       }
+      pc->force[ia] += f;
+      forcewall[ia] -= f;
     }
   }
 

src/wall_ss_cut.h

@@ -5,7 +5,7 @@
  *  Edinburgh Soft Matter and Statistical Physics Group and
  *  Edinburgh Parallel Computing Centre
  *
- *  (c) The University of Edinburgh (2014)
+ *  (c) 2014-2022 The University of Edinburgh
  *  Juho Lintuvuori ([email protected])
  *
  *****************************************************************************/
@@ -19,17 +19,24 @@
 #include "colloids.h"
 #include "interaction.h"
 
+typedef struct wall_ss_cut_options_s wall_ss_cut_options_t;
 typedef struct wall_ss_cut_s wall_ss_cut_t;
 
+struct wall_ss_cut_options_s {
+  double epsilon;       /* Energy scale */
+  double sigma;         /* Sigma (length) */
+  double nu;            /* exponent */
+  double hc;            /* Surface-surface cut off */
+};
+
 int wall_ss_cut_create(pe_t * pe, cs_t * cs, wall_t * wall,
+		       const wall_ss_cut_options_t * opts,
 		       wall_ss_cut_t ** pobj);
 int wall_ss_cut_free(wall_ss_cut_t * obj);
 int wall_ss_cut_info(wall_ss_cut_t * obj);
 int wall_ss_cut_register(wall_ss_cut_t * obj, interact_t * parent);
 int wall_ss_cut_compute(colloids_info_t * cinfo, void * self);
 int wall_ss_cut_stats(void * self, double * stats);
 int wall_ss_cut_single(wall_ss_cut_t * obj, double h, double * f, double * v);
-int wall_ss_cut_param_set(wall_ss_cut_t * obj, double epsilon, double sigma,
-			  int nu, double hc);
 
 #endif