Merge pull request #519 from Pakketeretet2/USER-MANIFOLD-gaussian-bump

Some extensions/cleanup in USER-MANIFOLD

doc/src/manifolds.txt

@@ -24,8 +24,9 @@ to the relevant fixes.
 {manifold} @ {parameters} @ {equation} @ {description}
 cylinder @ R @ x^2 + y^2 - R^2 = 0 @ Cylinder along z-axis, axis going through (0,0,0)
 cylinder_dent @ R l a @ x^2 + y^2 - r(z)^2 = 0, r(x) = R if | z | > l, r(z) = R - a*(1 + cos(z/l))/2 otherwise @ A cylinder with a dent around z = 0
-dumbbell @ a A B c @ -( x^2 + y^2 ) * (a^2 - z^2/c^2) * ( 1 + (A*sin(B*z^2))^4) = 0 @ A dumbbell @
+dumbbell @ a A B c @ -( x^2 + y^2 ) + (a^2 - z^2/c^2) * ( 1 + (A*sin(B*z^2))^4) = 0 @ A dumbbell
 ellipsoid @ a  b c @ (x/a)^2 + (y/b)^2 + (z/c)^2 = 0 @ An ellipsoid
+gaussian_bump @ A l rc1 rc2 @ if( x < rc1) -z + A * exp( -x^2 / (2 l^2) ); else if( x < rc2 ) -z + a + b*x + c*x^2 + d*x^3; else z @ A Gaussian bump at x = y = 0, smoothly tapered to a flat plane z = 0.
 plane @ a b c x0 y0 z0 @ a*(x-x0) + b*(y-y0) + c*(z-z0) = 0 @ A plane with normal (a,b,c) going through point (x0,y0,z0)
 plane_wiggle @ a w @ z - a*sin(w*x) = 0 @ A plane with a sinusoidal modulation on z along x.
 sphere @ R @ x^2 + y^2 + z^2 - R^2 = 0 @ A sphere of radius R

src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp

@@ -115,6 +115,13 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
     error->all(FLERR, "Error creating manifold arg arrays");
   }
 
+  // Check if you have enough args:
+  if( 6 + nvars > narg ){
+    char msg[2048];
+    sprintf(msg, "Not enough args for manifold %s, %d expected but got %d\n",
+            ptr_m->id(), nvars, narg - 6);
+    error->all(FLERR, msg);
+  }
   // Loop over manifold args:
   for( int i = 0; i < nvars; ++i ){
     int len = 0, offset = 0;

src/USER-MANIFOLD/manifold.h

@@ -24,7 +24,7 @@
    testing purposes) and a wave-y plane.
    See the README file for more info.
 
-   Stefan Paquay, s.paquay@tue.nl
+   Stefan Paquay, stefanpaquay@gmail.com
    Applied Physics/Theory of Polymers and Soft Matter,
    Eindhoven University of Technology (TU/e), The Netherlands
 

src/USER-MANIFOLD/manifold_factory.cpp

@@ -37,6 +37,7 @@
 #include "manifold_cylinder_dent.h"
 #include "manifold_dumbbell.h"
 #include "manifold_ellipsoid.h"
+#include "manifold_gaussian_bump.h"
 #include "manifold_plane.h"
 #include "manifold_plane_wiggle.h"
 #include "manifold_sphere.h"
@@ -58,12 +59,13 @@ manifold* LAMMPS_NS::user_manifold::create_manifold(const char *mname,
   make_manifold_if<manifold_cylinder_dent>  ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_dumbbell>       ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_ellipsoid>      ( &man, mname, lmp, narg, arg );
+  make_manifold_if<manifold_gaussian_bump>  ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_plane>          ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_plane_wiggle>   ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_sphere>         ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_supersphere>    ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_spine>          ( &man, mname, lmp, narg, arg );
-  make_manifold_if<manifold_spine_two>    ( &man, mname, lmp, narg, arg );
+  make_manifold_if<manifold_spine_two>      ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_thylakoid>      ( &man, mname, lmp, narg, arg );
   make_manifold_if<manifold_torus>          ( &man, mname, lmp, narg, arg );
 

src/USER-MANIFOLD/manifold_gaussian_bump.cpp

@@ -0,0 +1,138 @@
+#include "manifold_gaussian_bump.h"
+
+using namespace LAMMPS_NS;
+using namespace user_manifold;
+
+// The constraint is z = f(r = sqrt( x^2 + y^2) )
+// f(x) = A*exp( -x^2 / 2 l^2 )       if x < rc1
+//      = a + b*x + c*x**2 + d*x**3   if rc1 <= x < rc2
+//      = 0                           if x >= rc2
+//
+double manifold_gaussian_bump::g( const double *x )
+{
+	double xf[3];
+	xf[0] = x[0];
+	xf[1] = x[1];
+	xf[2] = 0.0;
+	
+	double x2 = dot(xf,xf);
+	if( x2 < rc12 ){
+		return x[2] - gaussian_bump_x2( x2 );
+	}else if( x2 < rc22 ){
+		double rr = sqrt(x2);
+		double xi = rr - rc1;
+		xi *= inv_dr;
+		double xi2 = x2 * inv_dr*inv_dr;
+		double xi3 = xi*xi2;
+		return x[2] - ( aa + bb*xi + cc*xi2 + dd*xi3 );
+		
+	}else{
+		return x[2];
+	}
+}
+
+void   manifold_gaussian_bump::n( const double *x, double *nn )
+{
+	double xf[3];
+	xf[0] = x[0];
+	xf[1] = x[1];
+	xf[2] = 0.0;
+	nn[2] = 1.0;
+	
+	double x2 = dot(xf,xf);
+	
+	if( x2 < rc12 ){
+		double factor = gaussian_bump_x2(x2);
+		factor /= (ll*ll);
+		nn[0] = factor * x[0];
+		nn[1] = factor * x[1];
+	}else if( x2 < rc22 ){
+		double rr = sqrt(x2);
+		double xi = rr - rc1;
+		xi *= inv_dr;
+		double xi2 = x2 * inv_dr*inv_dr;
+		double der = bb + 2*cc*xi + 3*dd*xi2;
+		
+		nn[0] = -der * x[0] / rr;
+		nn[1] = -der * x[1] / rr;
+	}else{
+		nn[0] = nn[1] = 0.0;
+	}
+}
+
+double manifold_gaussian_bump::g_and_n( const double *x, double *nn )
+{
+	double xf[3];
+	xf[0] = x[0];
+	xf[1] = x[1];
+	xf[2] = 0.0;
+	nn[2] = 1.0;
+	
+	double x2 = dot(xf,xf);
+	if( x2 < rc12 ){
+		double gb = gaussian_bump_x2(x2);
+		double factor = gb / (ll*ll);
+		nn[0] = factor * x[0];
+		nn[1] = factor * x[1];
+		
+		return x[2] - gb;
+	}else if( x2 < rc22 ){
+		
+		double rr = sqrt(x2);
+		double xi = rr - rc1;
+		xi *= inv_dr;
+		double xi2 = x2 * inv_dr*inv_dr;
+		double xi3 = xi*xi2;
+
+		double der = bb + 2*cc*xi + 3*dd*xi2;
+		
+		nn[0] = -der * x[0] / rr;
+		nn[1] = -der * x[1] / rr;
+
+		
+		return x[2] - ( aa + bb*xi + cc*xi2 + dd*xi3 );
+		
+	}else{
+		nn[0] = nn[1] = 0.0;
+		return x[2];
+	}
+	
+}
+
+
+void manifold_gaussian_bump::post_param_init()
+{
+	// Read in the params:
+	AA  = params[0];
+	ll  = params[1];
+	rc1 = params[2];
+	rc2 = params[3];
+
+	ll2 = 2.0*ll*ll;
+
+	f_at_rc  = gaussian_bump_x2 ( rc12 );
+	fp_at_rc = gaussian_bump_der( rc12 );
+
+	rc12 = rc1*rc1;
+	rc22 = rc2*rc2;
+	dr = rc2 - rc1;
+	inv_dr = 1.0 / dr;
+}
+
+
+double manifold_gaussian_bump::gaussian_bump( double x )
+{
+	double x2 = x*x;
+	return gaussian_bump_x2( x2 );
+}
+
+double manifold_gaussian_bump::gaussian_bump_x2( double x2 )
+{
+	return AA*exp( -x2 / ll2 );
+}
+
+double manifold_gaussian_bump::gaussian_bump_der( double x )
+{
+	double x2 = x*x;
+	return gaussian_bump_x2( x2 )*( -x/(ll*ll) );
+}

src/USER-MANIFOLD/manifold_gaussian_bump.h

@@ -0,0 +1,82 @@
+/* ----------------------------------------------------------------------
+   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.
+   -----------------------------------------------------------------------
+
+   This file is a part of the USER-MANIFOLD package.
+
+   This package allows LAMMPS to perform MD simulations of particles
+   constrained on a manifold (i.e., a 2D subspace of the 3D simulation
+   box). It achieves this using the RATTLE constraint algorithm applied
+   to single-particle constraint functions g(xi,yi,zi) = 0 and their
+   derivative (i.e. the normal of the manifold) n = grad(g).
+
+   It is very easy to add your own manifolds to the current zoo
+   (we now have sphere, a dendritic spine approximation, a 2D plane (for
+   testing purposes) and a wave-y plane.
+   See the README file for more info.
+
+   Stefan Paquay, spaquay@brandeis.edu
+   Brandeis University, Waltham, MA, USA.
+
+   This package was mainly developed at
+   Applied Physics/Theory of Polymers and Soft Matter,
+   Eindhoven University of Technology (TU/e), The Netherlands
+
+   Thanks to Remy Kusters at TU/e for testing.
+
+   This software is distributed under the GNU General Public License.
+
+------------------------------------------------------------------------- */
+
+#ifndef LMP_MANIFOLD_GAUSSIAN_BUMP_H
+#define LMP_MANIFOLD_GAUSSIAN_BUMP_H
+
+#include "manifold.h"
+
+namespace LAMMPS_NS {
+
+namespace user_manifold {
+
+  // A Gaussian bump with a smoothed decay to flat 2D.
+  class manifold_gaussian_bump : public manifold {
+   public:
+    enum { NPARAMS = 4 };
+    manifold_gaussian_bump(class LAMMPS*, int, char **) : manifold(lmp) {}
+    virtual ~manifold_gaussian_bump(){}
+
+    virtual double g( const double * );
+    virtual void   n( const double *, double * );
+
+    // Variant of g that computes n at the same time.
+    virtual double g_and_n( const double *x, double *nn );
+
+    static const char* type(){ return "gaussian_bump"; }
+    virtual const char *id() { return type(); }
+
+    virtual int nparams(){ return NPARAMS; }
+    virtual void post_param_init();
+   private:
+    // Some private constants:
+    double aa, bb, cc, dd, AA, ll, ll2, f_at_rc, fp_at_rc;
+    double rc1, rc2, rc12, rc22, dr, inv_dr;
+
+    double gaussian_bump    ( double );
+    double gaussian_bump_x2 ( double );
+    double gaussian_bump_der( double );
+
+  };
+}
+
+}
+
+
+#endif // LMP_MANIFOLD_GAUSSIAN_BUMP_H