diff --git a/doc/src/manifolds.txt b/doc/src/manifolds.txt index eb3bd6d486..c9bb1ce57f 100644 --- a/doc/src/manifolds.txt +++ b/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 diff --git a/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp b/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp index 4f6b62590d..67f298d8df 100644 --- a/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp +++ b/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; diff --git a/src/USER-MANIFOLD/manifold.h b/src/USER-MANIFOLD/manifold.h index d0ffa214ac..b89e765a6e 100644 --- a/src/USER-MANIFOLD/manifold.h +++ b/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 diff --git a/src/USER-MANIFOLD/manifold_factory.cpp b/src/USER-MANIFOLD/manifold_factory.cpp index 57c6e8305f..f98180ddb6 100644 --- a/src/USER-MANIFOLD/manifold_factory.cpp +++ b/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 ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); + make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); - make_manifold_if ( &man, mname, lmp, narg, arg ); + make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); make_manifold_if ( &man, mname, lmp, narg, arg ); diff --git a/src/USER-MANIFOLD/manifold_gaussian_bump.cpp b/src/USER-MANIFOLD/manifold_gaussian_bump.cpp new file mode 100644 index 0000000000..200edd2fb0 --- /dev/null +++ b/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) ); +} diff --git a/src/USER-MANIFOLD/manifold_gaussian_bump.h b/src/USER-MANIFOLD/manifold_gaussian_bump.h new file mode 100644 index 0000000000..f3401a4a33 --- /dev/null +++ b/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