diff --git a/src/Makefile b/src/Makefile index 892810e28e..16d3e5326c 100755 --- a/src/Makefile +++ b/src/Makefile @@ -48,7 +48,7 @@ PACKAGE = asphere body class2 colloid compress coreshell dipole gpu \ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars user-cuda \ user-diffraction user-dpd user-drude user-eff user-fep user-h5md \ - user-intel user-lb user-mgpt \ + user-intel user-lb user-manifold user-mgpt \ user-misc user-molfile user-omp user-phonon user-qmmm user-qtb \ user-quip user-reaxc user-smd user-smtbq user-sph user-tally \ user-vtk diff --git a/src/USER-MANIFOLD/README b/src/USER-MANIFOLD/README new file mode 100644 index 0000000000..3656391641 --- /dev/null +++ b/src/USER-MANIFOLD/README @@ -0,0 +1,84 @@ ++==============================================================================+ +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). + +Stefan Paquay, s.paquay@tue.nl +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. ++==============================================================================+ + +At the moment we have a few manifolds available, extending them is very easy: +To add a new manifold, do the following in the "USER-MANIFOLD" directory: + +0. Create a new pair of source/header files, and name them "manifold_*.cpp/h", + where you should replace '*' with some (descriptive) name. +1. In the header file, add the following: + a. Include guards. + b. Make sure you #include "manifold.h" + c. In namespace LAMMPS_NS, add a new class that inherits publicly from manifold + and protectedly from Pointers. + d. The constructor has to take ( LAMMPS*, int, char ** ) as arguments, + and should initialise Pointers. + e. The header file has to contain somewhere the macro ManifoldStyle with as + first argument the name of the manifold and as second argument the name + of the class implementing this manifold. The macro expands into some code + that registers the manifold during static initialisation, before main is + entered. + +2. In the source file, make sure you implement the following (of course, + you can also implement these in the header): ++====================================+=========================================+ +| Function signature | Purpose | ++====================================+=========================================+ +| destructor | Free space (can be empty) | +| constructor | Has to take (LAMMPS *lmp, int, char **) | +| | as arguments and has to initialise | +| | Pointers with lmp. | +| double g( double *x ) | A function that is 0 if the 3-vector | +| | x is on the manifold. | +| void n( double *x, double *n ) | Stores the normal of position x in the | +| | 3-vector n. | +| static const char *type() | Returns text that identifies the | +| | manifold in LAMMPS input scripts. | +| const char *id() | Should return whatever type() returns. | +| static int expected_argc() | Returns the number of arguments needed | +| | for the construction/initialisation of | +| | your manifold. Example: Sphere only | +| | needs a radius, so it returns 1. The | +| | spine needs 5 parameters, so it | +| | returns 5. | +| int nparams() | Should return same as expected_argc() | ++====================================+=========================================+ + +If the above instructions seem a bit intimidating, you can get by just fine +by copying an existing manifold and modifying it. See e.g. manifold_sphere for +a relatively simple manifold. + +With those things in place, the install script should be able to add your +manifold to LAMMPS without any extra work, so just running +make yes-user-manifold +make +should (re)compile LAMMPS with the manifolds added. + ++==============================================================================+ +Obviously, you need some parameters to represent the surface, such as the radius +in the case of a sphere. These should be passed to the nve/manifold/rattle fix, +with the following syntax: +fix ID group-ID nve/manifold/rattle tol maxit manifold_style args + +tol = tolerance to which RATTLE tries to satisfy the constraints +maxit = maximum number of iterations RATTLE uses each time step +manifold_style = manifold style, should equal what type() of the desired + manifold returns +args = parameters for the manifold, order is manifold-dependent. + can be equal style variables ++==============================================================================+ diff --git a/src/USER-MANIFOLD/fix_manifoldforce.cpp b/src/USER-MANIFOLD/fix_manifoldforce.cpp new file mode 100644 index 0000000000..da794e4a48 --- /dev/null +++ b/src/USER-MANIFOLD/fix_manifoldforce.cpp @@ -0,0 +1,185 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "string.h" +#include "stdlib.h" +#include "atom.h" +#include "update.h" +#include "respa.h" +#include "error.h" +#include "force.h" + +#include "manifold.h" +#include "fix_manifoldforce.h" // For stuff +#include "manifold_factory.h" // For constructing manifold + + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace user_manifold; + + +// Helper functions for parameters/equal style variables in input script +inline bool was_var( const char *arg ) +{ + return strstr( arg, "v_" ) == arg; +} + +inline bool str_eq( const char *str1, const char *str2 ) +{ + return strcmp(str1,str2) == 0; +} + +/* ---------------------------------------------------------------------- */ + +FixManifoldForce::FixManifoldForce(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + int me = -1; + MPI_Comm_rank(world,&me); + + + // Check the min-style: + int good_minner = str_eq(update->minimize_style,"hftn") | + str_eq(update->minimize_style,"quickmin"); + if( !good_minner){ + error->warning(FLERR,"Minimizing with fix manifoldforce without hftn or quickmin is fishy"); + } + + + // Command is given as + // fix manifoldforce manifold_name manifold_args + if( narg < 5 ){ + error->all(FLERR,"Illegal fix manifoldforce! No manifold given"); + } + const char *m_name = arg[3]; + ptr_m = create_manifold(m_name,lmp,narg,arg); + + // Construct manifold from factory: + if( !ptr_m ){ + char msg[2048]; + snprintf(msg,2048,"Manifold pointer for manifold '%s' was NULL for some reason", arg[3]); + error->all(FLERR,msg); + } + + + // After constructing the manifold, you can safely make + // room for the parameters + nvars = ptr_m->nparams(); + if( narg < nvars+4 ){ + char msg[2048]; + sprintf(msg,"Manifold %s needs at least %d argument(s)!", + m_name, nvars); + error->all(FLERR,msg); + } + + *(ptr_m->get_params()) = new double[nvars]; + if( ptr_m->get_params() == NULL ){ + error->all(FLERR,"Parameter pointer was NULL!"); + } + + // This part here stores the names/text of each argument, + // determines which params are equal-style variables, + // and sets the values of those arguments that were _not_ + // equal style vars (so that they are not overwritten each time step). + + double *params = *(ptr_m->get_params()); + for( int i = 0; i < nvars; ++i ){ + if( was_var( arg[i+4] ) ) + error->all(FLERR,"Equal-style variables not allowed with fix manifoldforce"); + + // Use force->numeric to trigger an error if arg is not a number. + params[i] = force->numeric(FLERR,arg[i+4]); + } + + + // Perform any further initialisation for the manifold that depends on params: + ptr_m->post_param_init(); +} + +/* ---------------------------------------------------------------------- */ + +int FixManifoldForce::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + mask |= MIN_POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixManifoldForce::setup(int vflag) +{ + if (strstr(update->integrate_style,"verlet")) + post_force(vflag); + else { + int nlevels_respa = ((Respa *) update->integrate)->nlevels; + for (int ilevel = 0; ilevel < nlevels_respa; ilevel++) { + ((Respa *) update->integrate)->copy_flevel_f(ilevel); + post_force_respa(vflag,ilevel,0); + ((Respa *) update->integrate)->copy_f_flevel(ilevel); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixManifoldForce::min_setup(int vflag) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixManifoldForce::post_force(int vflag) +{ + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double n[3]; + double invn2; + double dot; + for (int i = 0; i < nlocal; i++){ + if (mask[i] & groupbit) { + // Determine normal of particle: + ptr_m->n(x[i],n); + + invn2 = 1.0 / ( n[0]*n[0] + n[1]*n[1] + n[2]*n[2] ); + dot = f[i][0]*n[0] + f[i][1]*n[1] + f[i][2]*n[2]; + + f[i][0] -= dot*n[0] * invn2; + f[i][1] -= dot*n[1] * invn2; + f[i][2] -= dot*n[2] * invn2; + + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixManifoldForce::post_force_respa(int vflag, int ilevel, int iloop) +{ + post_force(vflag); +} + +/* ---------------------------------------------------------------------- */ + +void FixManifoldForce::min_post_force(int vflag) +{ + post_force(vflag); +} diff --git a/src/USER-MANIFOLD/fix_manifoldforce.h b/src/USER-MANIFOLD/fix_manifoldforce.h new file mode 100644 index 0000000000..89bd6bd378 --- /dev/null +++ b/src/USER-MANIFOLD/fix_manifoldforce.h @@ -0,0 +1,81 @@ +/* -*- c++ -*- ---------------------------------------------------------- + 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. + + Copyright (2013-2015) Stefan Paquay, Eindhoven University of Technology. + License: GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + This file is part of the user-manifold package written by + Stefan Paquay at the Eindhoven University of Technology. + This module makes it possible to do MD with particles constrained + to pretty arbitrary manifolds characterised by some constraint function + g(x,y,z) = 0 and its normal grad(g). The number of manifolds available + right now is limited but can be extended straightforwardly by making + a new class that inherits from manifold and implements all pure virtual + methods. + + This fix subtracts force components that point out of the manifold, + useful for minimisation on surfaces. + + +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(manifoldforce,FixManifoldForce) + +#else + +#ifndef LMP_FIX_MANIFOLDFORCE_H +#define LMP_FIX_MANIFOLDFORCE_H + +#include "fix.h" +#include "manifold.h" + +namespace LAMMPS_NS { + +class FixManifoldForce : public Fix { + public: + FixManifoldForce(class LAMMPS *, int, char **); + int setmask(); + void setup(int); + void min_setup(int); + void post_force(int); + void post_force_respa(int, int, int); + void min_post_force(int); + + + private: + user_manifold::manifold *ptr_m; + + // Stuff to store the parameters in. + int nvars; // # of args after manifold name. +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +*/ diff --git a/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp b/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp new file mode 100644 index 0000000000..d2d73de6bd --- /dev/null +++ b/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp @@ -0,0 +1,641 @@ +/* ---------------------------------------------------------------------- + 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. + + Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology. + License: GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + This file is part of the user-manifold package written by + Stefan Paquay at the Eindhoven University of Technology. + This module makes it possible to do MD with particles constrained + to pretty arbitrary manifolds characterised by some constraint function + g(x,y,z) = 0 and its normal grad(g). The number of manifolds available + right now is limited but can be extended straightforwardly by making + a new class that inherits from manifold and implements all pure virtual + methods. + + Thanks to Remy Kusters for beta-testing! + +------------------------------------------------------------------------- */ + + +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "atom.h" +#include "force.h" +#include "update.h" +#include "respa.h" +#include "error.h" +#include "group.h" +#include "math.h" +#include "input.h" +#include "variable.h" +#include "citeme.h" +#include "memory.h" +#include "comm.h" + + +#include "fix_nve_manifold_rattle.h" +#include "manifold_factory.h" +#include "manifold.h" + + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace user_manifold; + + +enum { CONST, EQUAL }; // For treating the variables. + + +static const char* cite_fix_nve_manifold_rattle = + "fix nve/manifold/rattle command:\n\n" + "@article{paquay-2016,\n" + " author = {Paquay, Stefan and Kusters, Remy}, \n" + " eprint = {arXiv:1411.3019}, \n" + " title = {A method for molecular dynamics on curved surfaces}, \n" + " url = {http://arxiv.org/abs/1411.3019}, \n" + " year = 2016, \n" + " journal = {Arxiv preprint, \\url{http://arxiv.org/abs/1411.3019}}, \n" + " note = {To be published in Biophys. J.}, \n" + "}\n\n"; + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg, + int error_on_unknown_keyword ) + : Fix(lmp,narg,arg) +{ + if( lmp->citeme) lmp->citeme->add(cite_fix_nve_manifold_rattle); + if( narg < 6 ) error->all(FLERR, "Illegal fix nve/manifold/rattle command"); + + // Set all bits/settings: + time_integrate = 1; + dynamic_group_allow = 1; + size_vector = 0; + dof_flag = 1; + + nevery = 0; + dtv = dtf = 0; + + tolerance = force->numeric( FLERR, arg[3] ); + max_iter = force->numeric( FLERR, arg[4] ); + + ptr_m = create_manifold(arg[5], lmp, narg, arg); + if( !ptr_m ){ + error->all(FLERR,"Error creating manifold pointer"); + } + + nvars = ptr_m->nparams(); + tstrs = new char*[nvars]; + tvars = new int[nvars]; + tstyle = new int[nvars]; + is_var = new int[nvars]; + + if( !tstrs || !tvars || !tstyle || !is_var ){ + error->all(FLERR, "Error creating manifold arg arrays"); + } + + // Loop over manifold args: + for( int i = 0; i < nvars; ++i ){ + int len = 0, offset = 0; + if( was_var( arg[i+6] ) ){ + len = strlen(arg[i+6]) - 1; // -1 because -2 for v_, +1 for \0. + is_var[i] = 1; + offset = 2; + }else{ + force->numeric(FLERR,arg[i+6]); // Check if legal number. + len = strlen( arg[i+6] ) + 1; // +1 for \0. + is_var[i] = 0; + } + tstrs[i] = new char[len]; + if( tstrs[i] == NULL ) error->all(FLERR,"Error allocating space for args."); + strcpy( tstrs[i], arg[i+6] + offset ); + } + + *ptr_m->get_params() = new double[nvars]; + if( !(*ptr_m->get_params()) ) error->all(FLERR,"Failed to allocate params!"); + for( int i = 0; i < nvars; ++i ){ + // If param i was variable type, it will be set later... + (*ptr_m->get_params())[i] = is_var[i] ? 0.0 : force->numeric( FLERR, arg[i+6] ); + } + ptr_m->post_param_init(); + + + // Loop over rest of args: + int argi = 6 + nvars; + while( argi < narg ){ + if( strcmp(arg[argi], "every") == 0 ){ + nevery = force->inumeric(FLERR,arg[argi+1]); + argi += 2; + }else if( error_on_unknown_keyword ){ + char msg[2048]; + sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]); + error->all(FLERR, msg); + }else{ + argi += 1; + } + } + +} + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +FixNVEManifoldRattle::~FixNVEManifoldRattle() +{ + if( tstrs ){ + for( int i = 0; i < nvars; ++i ){ + delete [] tstrs[i]; + } + delete [] tstrs; + } + + if( tvars ) delete [] tvars; + if( tstyle ) delete [] tstyle; + if( is_var ) delete [] is_var; +} + + + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::reset_dt() +{ + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; + + + +} + +void FixNVEManifoldRattle::print_stats( const char *header ) +{ + double n = stats.natoms; + if( n > 0 ){ + stats.x_iters_per_atom += stats.x_iters / n; + stats.v_iters_per_atom += stats.v_iters / n; + } + + double x_iters = 0, v_iters = 0; + bigint ntimestep = update->ntimestep; + int me = -1; + + + MPI_Comm_rank(world,&me); + MPI_Allreduce(&stats.x_iters_per_atom,&x_iters,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&stats.v_iters_per_atom,&v_iters,1,MPI_DOUBLE,MPI_SUM,world); + + // Set iters back to zero: + stats.x_iters_per_atom = stats.x_iters = 0; + stats.v_iters_per_atom = stats.v_iters = 0; + + + if( me == 0 ){ + double inv_tdiff = 1.0/( static_cast(ntimestep) - stats.last_out ); + stats.last_out = ntimestep; + + fprintf(screen, "%s stats for time step " BIGINT_FORMAT " on %d atoms:\n", + header, ntimestep, stats.natoms); + fprintf(screen, " iters/atom: x = %f, v = %f, dofs removed %d", + x_iters * inv_tdiff, v_iters * inv_tdiff, stats.dofs_removed); + fprintf(screen,"\n"); + } +} + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +int FixNVEManifoldRattle::was_var( const char *str ) +{ + if( strlen(str) > 2 ){ + return (str[0] == 'v') && (str[1] == '_'); + }else{ + return 0; + } +} + + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +int FixNVEManifoldRattle::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + if( nevery > 0 ) mask |= END_OF_STEP; + + return mask; +} + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::init() +{ + // Makes sure the manifold params are set initially. + + update_var_params(); + reset_dt(); +} + + +void FixNVEManifoldRattle::update_var_params() +{ + if( nevery > 0 ){ + stats.x_iters = 0; + stats.v_iters = 0; + stats.natoms = 0; + stats.x_iters_per_atom = 0.0; + stats.v_iters_per_atom = 0.0; + } + + double **ptr_params = ptr_m->get_params(); + for( int i = 0; i < nvars; ++i ){ + if( is_var[i] ){ + tvars[i] = input->variable->find(tstrs[i]); + if( tvars[i] < 0 ){ + error->all(FLERR, + "Variable name for fix nve/manifold/rattle does not exist"); + } + if( input->variable->equalstyle(tvars[i]) ){ + tstyle[i] = EQUAL; + double new_val = input->variable->compute_equal(tvars[i]); + // fprintf( stdout, "New value of var %d is now %f\n", i+1, new_val ); + *(ptr_params[i]) = new_val; + }else{ + error->all(FLERR, + "Variable for fix nve/manifold/rattle is invalid style"); + } + } + } +} + + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +int FixNVEManifoldRattle::dof(int igroup) +{ + int *mask = atom->mask; + int nlocal = atom->nlocal; + int natoms = 0; + for( int i = 0; i < nlocal; ++i ){ + if(mask[i] & groupbit) ++natoms; + } + + int dofs; + MPI_Allreduce( &natoms, &dofs, 1, MPI_INT, MPI_SUM, world ); + + // Make sure that, if there is just no or one atom, no dofs are subtracted, + // since for the first atom already 3 dofs are subtracted because of the + // centre of mass corrections: + if( dofs <= 1 ) dofs = 0; + stats.dofs_removed = dofs; + + return dofs; +} + + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +double FixNVEManifoldRattle::memory_usage() +{ + double bytes = 0.0; + + bytes += sizeof(statistics); + bytes += sizeof(*ptr_m) + sizeof(ptr_m); + bytes += nvars*sizeof(double) + sizeof(double*); + bytes += nvars*( sizeof(char*) + 3*sizeof(int) ); + return bytes; +} + + + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::initial_integrate(int vflag) +{ + update_var_params(); + nve_x_rattle(igroup, groupbit); +} + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::final_integrate() +{ + nve_v_rattle(igroup, groupbit); +} + + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::end_of_step() +{ + print_stats( "nve/manifold/rattle" ); +} + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::nve_x_rattle(int igroup, int groupbit) +{ + double dtfm; + // update v and x of atoms in group + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int natoms = 0; + + if (igroup == atom->firstgroup){ + nlocal = atom->nfirst; + } + + + if (rmass) { + for (int i = 0; i < nlocal; i++){ + if (mask[i] & groupbit){ + natoms++; + dtfm = dtf / rmass[i]; + rattle_manifold_x( x[i], v[i], f[i], dtv, dtfm, atom->tag[i] ); + } + } + } else { + for (int i = 0; i < nlocal; i++){ + if (mask[i] & groupbit) { + natoms++; + dtfm = dtf / mass[type[i]]; + rattle_manifold_x( x[i], v[i], f[i], dtv, dtfm, atom->tag[i] ); + } + } + } + + if( nevery > 0 ){ + // Count ALL atoms this fix works on: + MPI_Allreduce(&natoms,&stats.natoms,1,MPI_INT,MPI_SUM,world); + } +} + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::nve_v_rattle(int igroup, int groupbit) +{ + double dtfm; + + // update v of atoms in group + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + if (rmass) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + rattle_manifold_v( v[i], f[i], x[i], dtfm, atom->tag[i] ); + } + } + } else { + for (int i = 0; i < nlocal; i++){ + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + rattle_manifold_v( v[i], f[i], x[i], dtfm, atom->tag[i] ); + } + } + } +} + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::rattle_manifold_x(double *x, double *v, + double *f, double dtv, + double dtfm, tagint tagi ) +{ + /* + A RATTLE update for the position constraint. + Original update is x += dtv * v_1/2 + Now you do + v_1/2(lambda) = v_0 + dtfm * ( f + lambda*n_old ) + and solve + xold - xnew + dtv * v_1/2(lambda) = 0 + g(xnew) = 0 + for x and lambda. The lambda you find then gives v_1/2 as well. + */ + double xo[3]; // Previous position to update from. + double vo[3]; // Previous velocity to update from. + double l = 0; // Lagrangian multiplier for constraint forces. + double R[4]; // System that is 0. + double dx[4]; // Update that follows from Newton iteration. + double no[3]; // Normal at xo. + double nn[3]; // Normal at x, the new position. + double res; // Residual. + int iters = 0; // Iterations used + + double c = dtfm*dtv; // Used for iterating in the Newton loop: + double no_nn, nn_R; + + vo[0] = v[0]; + vo[1] = v[1]; + vo[2] = v[2]; + + xo[0] = x[0]; + xo[1] = x[1]; + xo[2] = x[2]; + + double gg = ptr_m->g_and_n(x,no); + nn[0] = no[0]; + nn[1] = no[1]; + nn[2] = no[2]; + + double vt[3]; + vt[0] = vo[0] + dtfm*f[0]; + vt[1] = vo[1] + dtfm*f[1]; + vt[2] = vo[2] + dtfm*f[2]; + double no_dt[3]; + no_dt[0] = dtfm*no[0]; + no_dt[1] = dtfm*no[1]; + no_dt[2] = dtfm*no[2]; + + // Assume that no_nn is roughly constant during iteration: + + const double c_inv = 1.0 / c; + + + while ( 1 ) { + v[0] = vt[0] - l*no_dt[0]; + v[1] = vt[1] - l*no_dt[1]; + v[2] = vt[2] - l*no_dt[2]; + + R[0] = xo[0] - x[0] + dtv * v[0]; + R[1] = xo[1] - x[1] + dtv * v[1]; + R[2] = xo[2] - x[2] + dtv * v[2]; + R[3] = gg; + + // Analytic solution to system J*(dx,dy,dz,dl)^T = R + // no_nn = no[0]*nn[0] + no[1]*nn[1] + no[2]*nn[2]; + nn_R = nn[0]*R[0] + nn[1]*R[1] + nn[2]*R[2]; + no_nn = no[0]*nn[0] + no[1]*nn[1] + no[2]*nn[2]; + double n_inv = 1.0 / no_nn; + + // fprintf( screen, "nn_R = %f, no_nn = %f\n", nn_R, no_nn ); + + dx[3] = -nn_R - R[3]; + dx[3] *= n_inv; + dx[0] = -R[0] - no[0]*dx[3]; + dx[1] = -R[1] - no[1]*dx[3]; + dx[2] = -R[2] - no[2]*dx[3]; + + dx[3] *= c_inv; + + + x[0] -= dx[0]; + x[1] -= dx[1]; + x[2] -= dx[2]; + l -= dx[3]; + + res = infnorm<4>(R); + ++iters; + + if( (res < tolerance) || (iters >= max_iter) ) break; + + // Update nn and g. + gg = ptr_m->g(x); + ptr_m->n(x,nn); + // gg = ptr_m->g(x); + } + + if( iters >= max_iter && res > tolerance ){ + char msg[2048]; + sprintf(msg,"Failed to constrain atom %d (x = (%f, %f, %f)! res = %e, iters = %d\n", + tagi, x[0], x[1], x[2], res, iters); + error->one(FLERR,msg); + } + + // "sync" x and v: + v[0] = vt[0] - l*no_dt[0]; + v[1] = vt[1] - l*no_dt[1]; + v[2] = vt[2] - l*no_dt[2]; + + stats.x_iters += iters; +} + + +/* ----------------------------------------------------------------------------- + ---------------------------------------------------------------------------*/ +void FixNVEManifoldRattle::rattle_manifold_v(double *v, double *f, + double *x, double dtfm, + tagint tagi ) +{ + /* + The original update was + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + Now you add the rattle-like update: + vold - vnew + dtfm * F + mu * n_new = 0 + dot( vnew, n_new ) = 0 + */ + double vo[3]; // V at t + 1/2 dt + double l = 0; // Lagrangian multiplier for constraint forces. + double R[4]; // System that is 0. + double dv[4]; // Update that follows from Newton iteration. + double n[3]; // Normal. + double res; // Residual. + int iters = 0; // Iterations used + + double c = dtfm; // Used for iterating in the Newton loop: + double nn2, nn_R; + + vo[0] = v[0]; + vo[1] = v[1]; + vo[2] = v[2]; + + // Initial guess is unconstrained update: + v[0] += dtfm*f[0]; + v[1] += dtfm*f[1]; + v[2] += dtfm*f[2]; + + ptr_m->n(x,n); + + double vt[3]; + vt[0] = vo[0] + dtfm*f[0]; + vt[1] = vo[1] + dtfm*f[1]; + vt[2] = vo[2] + dtfm*f[2]; + double no_dt[3]; + no_dt[0] = dtfm*n[0]; + no_dt[1] = dtfm*n[1]; + no_dt[2] = dtfm*n[2]; + + nn2 = n[0]*n[0] + n[1]*n[1] + n[2]*n[2]; + + const double n_inv = 1.0 / nn2; + const double c_inv = 1.0 / c; + + do{ + R[0] = vt[0] - v[0] - l * no_dt[0]; + R[1] = vt[1] - v[1] - l * no_dt[1]; + R[2] = vt[2] - v[2] - l * no_dt[2]; + R[3] = v[0]*n[0] + v[1]*n[1] + v[2]*n[2]; + + // Analytic solution to system J*(dx,dy,dz,dl)^T = R + nn_R = n[0]*R[0] + n[1]*R[1] + n[2]*R[2]; + + dv[3] = -nn_R - R[3]; + dv[3] *= n_inv; + dv[0] = -n[0]*dv[3] - R[0]; + dv[1] = -n[1]*dv[3] - R[1]; + dv[2] = -n[2]*dv[3] - R[2]; + dv[3] *= c_inv; + + v[0] -= dv[0]; + v[1] -= dv[1]; + v[2] -= dv[2]; + l -= dv[3]; + + res = infnorm<4>(R); + ++iters; + }while( (res > tolerance) && (iters < max_iter) ); + + if( iters >= max_iter && res >= tolerance ){ + char msg[2048]; + sprintf(msg,"Failed to constrain atom %d (x = (%f, %f, %f)! res = %e, iters = %d\n", + tagi, x[0], x[1], x[2], res, iters); + error->all(FLERR,msg); + } + + stats.v_iters += iters; +} diff --git a/src/USER-MANIFOLD/fix_nve_manifold_rattle.h b/src/USER-MANIFOLD/fix_nve_manifold_rattle.h new file mode 100644 index 0000000000..4bd17ab899 --- /dev/null +++ b/src/USER-MANIFOLD/fix_nve_manifold_rattle.h @@ -0,0 +1,172 @@ +/* ---------------------------------------------------------------------- + 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. + + Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology. + License: GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + This file is part of the user-manifold package written by + Stefan Paquay at the Eindhoven University of Technology. + This module makes it possible to do MD with particles constrained + to pretty arbitrary manifolds characterised by some constraint function + g(x,y,z) = 0 and its normal grad(g). The number of manifolds available + right now is limited but can be extended straightforwardly by making + a new class that inherits from manifold and implements all pure virtual + methods. + + Thanks to Remy Kusters for beta-testing! + +------------------------------------------------------------------------- */ + + +#ifdef FIX_CLASS + +FixStyle(nve/manifold/rattle,FixNVEManifoldRattle) + +#else + +#ifndef LMP_FIX_NVE_MANIFOLD_RATTLE_H +#define LMP_FIX_NVE_MANIFOLD_RATTLE_H + +#include "fix.h" +#include "manifold.h" + + +namespace LAMMPS_NS { + +// namespace user_manifold { + + class FixNVEManifoldRattle : public Fix { + public: + + struct statistics { + + statistics() : x_iters(0), v_iters(0), x_iters_per_atom(0), + v_iters_per_atom(0), natoms(0), dofs_removed(0), + last_out(0) {} + double x_iters, v_iters; + double x_iters_per_atom; + double v_iters_per_atom; + int natoms; + int dofs_removed; + bigint last_out; + }; + + FixNVEManifoldRattle(LAMMPS *, int &, char **, int = 1); + virtual ~FixNVEManifoldRattle(); + // All this stuff is interface, so you DO need to implement them. + // Just delegate them to the workhorse classes. + virtual int setmask(); + virtual void initial_integrate(int); + virtual void final_integrate(); + virtual void init(); + virtual void reset_dt(); + virtual void end_of_step(); + virtual int dof(int); + virtual void setup(int){} // Not needed for fixNVE but is for fixNVT + virtual double memory_usage(); + + protected: + + int nevery; + + double dtv, dtf; + double tolerance; + int max_iter; + + char **tstrs; + int *tvars; + int *tstyle; + int *is_var; + + statistics stats; + int update_style; + int nvars; + + user_manifold::manifold *ptr_m; + + + void print_stats( const char * ); + + int was_var( const char * ); + + virtual void update_var_params(); + virtual void rattle_manifold_x( double *, double *, double *, double, double, tagint ); + virtual void rattle_manifold_v( double *, double *, double *, double, tagint ); + + virtual void nve_x_rattle(int, int); + virtual void nve_v_rattle(int, int); + }; +} + + +#endif // LMP_FIX_NVE_MANIFOLD_RATTLE_H +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: There is no manifold named ... + +Self-explanatory. You requested a manifold whose name was not +registered at the factory. + +E: Manifold pointer was NULL for some reason! + +This indicates a bug. The factory was unable to properly create +the requested manifold even though it was registered. Send the +maintainer an e-mail. + +E: Manifold ... needs at least ... argument(s)! + +Self-explanatory. Provide enough arguments for the proper +creating of the requested manifold. + +E: Parameter pointer was NULL! + +This indicates a bug. The array that contains the parameters +could not be allocated. Send the maintainer an e-mail. + +E: Could not allocate space for arg! + +One of the arguments provided was too long (it contained +too many characters) + +E: Option ... needs ... argument(s): + +Self-explanatory. Read the documentation of this fix properly. + + +E: Illegal fix nve/manifold/rattle command! Option ... not recognized! + +Self-explanatory. The given option(s) do not exist. + +E: Variable name for fix nve/manifold/rattle does not exist + +Self-explanatory. + +E: Variable for fix nve/manifold/rattle is invalid style + +fix nve/manifold/rattle only supports equal style variables. + + + +*/ diff --git a/src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp b/src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp new file mode 100644 index 0000000000..3801508336 --- /dev/null +++ b/src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp @@ -0,0 +1,415 @@ +/* ---------------------------------------------------------------------- + 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. + + Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology. + License: GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + This file is part of the user-manifold package written by + Stefan Paquay at the Eindhoven University of Technology. + This module makes it possible to do MD with particles constrained + to pretty arbitrary manifolds characterised by some constraint function + g(x,y,z) = 0 and its normal grad(g). The number of manifolds available + right now is limited but can be extended straightforwardly by making + a new class that inherits from manifold and implements all pure virtual + methods. + + Thanks to Remy Kusters for beta-testing! + +------------------------------------------------------------------------- */ + + +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "atom.h" +#include "force.h" +#include "update.h" +#include "respa.h" +#include "error.h" +#include "group.h" +#include "math.h" +#include "input.h" +#include "variable.h" +#include "citeme.h" +#include "memory.h" +#include "comm.h" +#include "modify.h" +#include "compute.h" + +#include "fix_nvt_manifold_rattle.h" +#include "manifold.h" + + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace user_manifold; + +enum {CONSTANT,EQUAL}; +enum {NOBIAS,BIAS}; + + + + +static const char* cite_fix_nvt_manifold_rattle = + "fix nve/manifold/rattle command:\n\n" + "@article{paquay-2016,\n" + " author = {Paquay, Stefan and Kusters, Remy}, \n" + " eprint = {arXiv:1411.3019}, \n" + " title = {A method for molecular dynamics on curved surfaces}, \n" + " url = {http://arxiv.org/abs/1411.3019}, \n" + " year = 2016, \n" + " journal = {Arxiv preprint, \\url{http://arxiv.org/abs/1411.3019}}, \n" + " note = {To be published in Biophys. J.}, \n" + "}\n\n"; + + +/* ---------------------------------------------------------------------- */ + +FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg, + int error_on_unknown_keyword ) + : FixNVEManifoldRattle(lmp,narg,arg, 0) +{ + if (lmp->citeme) lmp->citeme->add(cite_fix_nvt_manifold_rattle); + + if( narg < 6 ) error->all(FLERR,"Illegal fix nvt/manifold/rattle command"); + + // Set all bits/settings: + dof_flag = 1; + dthalf = dt4 = dt8 = 0; + + t_start = t_stop = t_period = t_current = t_target = ke_target = 0.0; + t_freq = drag = tdrag_factor = 0; + + boltz = force->boltz, nktv2p = force->nktv2p; + tdof = 0; + mtchain = 3; + factor_eta = 0.0; + which = got_temp = 0; + + int argi = 6 + ptr_m->nparams(); + while( argi < narg ) + { + if( strcmp( arg[argi], "temp") == 0 ){ + if( argi+3 >= narg ) + error->all(FLERR,"Keyword 'temp' needs 3 arguments"); + + t_start = force->numeric(FLERR, arg[argi+1]); + t_stop = force->numeric(FLERR, arg[argi+2]); + t_period = force->numeric(FLERR, arg[argi+3]); + t_target = t_start; + got_temp = 1; + + argi += 4; + }else if( strcmp( arg[argi], "tchain" ) == 0 ){ + if( argi+1 >= narg ) + error->all(FLERR,"Keyword 'tchain' needs 1 argument"); + + mtchain = force->inumeric(FLERR, arg[argi+1]); + argi += 2; + }else if( error_on_unknown_keyword ){ + char msg[2048]; + sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]); + error->all(FLERR, msg); + }else{ + argi += 1; + } + } + + reset_dt(); + + if( !got_temp ) error->all(FLERR,"Fix nvt/manifold/rattle needs 'temp'!"); + + if( t_period < 0.0 ){ + error->all(FLERR,"Fix nvt/manifold/rattle damping parameter must be > 0.0"); + } + + // Create temperature compute: + const char *fix_id = arg[1]; + int n = strlen(fix_id)+6; + id_temp = new char[n]; + strcpy(id_temp,fix_id); + strcat(id_temp,"_temp"); + char **newarg = new char*[3]; + newarg[0] = id_temp; + newarg[1] = group->names[igroup]; + newarg[2] = (char*) "temp"; + + + modify->add_compute(3,newarg); + delete [] newarg; + int icompute = modify->find_compute(id_temp); + if( icompute < 0 ){ + error->all(FLERR,"Temperature ID for fix nvt/manifold/rattle " + "does not exist"); + } + temperature = modify->compute[icompute]; + if( temperature->tempbias ) which = BIAS; + else which = NOBIAS; + + // Set t_freq from t_period + t_freq = 1.0 / t_period; + + // Init Nosé-Hoover chain: + eta = new double[mtchain]; + eta_dot = new double[mtchain+1]; + eta_dotdot = new double[mtchain]; + eta_mass = new double[mtchain]; + eta_dot[mtchain] = 0.0; + + eta_dot[mtchain] = 0.0; + for( int ich = 0; ich < mtchain; ++ich ){ + eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0; + } + + +} + +/* ---------------------------------------------------------------------- */ + +FixNVTManifoldRattle::~FixNVTManifoldRattle() +{ + // Deallocate heap-allocated objects. + if( eta ) delete[] eta; + if( eta_dot ) delete[] eta_dot; + if( eta_dotdot ) delete[] eta_dotdot; + if( eta_mass ) delete[] eta_mass; + + modify->delete_compute(id_temp); + if( id_temp ) delete[] id_temp; +} + + + + +int FixNVTManifoldRattle::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + if( nevery > 0 ) mask |= END_OF_STEP; + + return mask; +} + + +/* -------------------------------------------------------------------------- + Check that force modification happens before position and velocity update. + Make sure respa is not used. +------------------------------------------------------------------------- */ +void FixNVTManifoldRattle::init() +{ + // Makes sure the manifold params are set initially. + update_var_params(); + + int icompute = modify->find_compute(id_temp); + if( icompute < 0 ){ + error->all(FLERR,"Temperature ID for fix nvt/manifold/rattle " + "does not exist"); + } + temperature = modify->compute[icompute]; + if( temperature->tempbias ) which = BIAS; + else which = NOBIAS; + +} + + + +void FixNVTManifoldRattle::setup(int vflag) +{ + compute_temp_target(); + + t_current = temperature->compute_scalar(); + tdof = temperature->dof; + + // Compute/set eta-masses: + double inv_t_freq2 = 1.0 / (t_freq*t_freq); + eta_mass[0] = tdof * boltz * t_target * inv_t_freq2; + for( int ich = 1; ich < mtchain; ++ich ){ + eta_mass[ich] = boltz * t_target * inv_t_freq2; + } + + for( int ich = 1; ich < mtchain; ++ich ){ + eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] - + boltz * t_target ) / eta_mass[ich]; + } +} + +void FixNVTManifoldRattle::compute_temp_target() +{ + + t_current = temperature->compute_scalar(); + tdof = temperature->dof; + + double delta = update->ntimestep - update->beginstep; + if (delta != 0.0){ + delta /= update->endstep - update->beginstep; + } + + tdof = temperature->dof; + t_target = t_start + delta * (t_stop-t_start); + ke_target = tdof * boltz * t_target; +} + +void FixNVTManifoldRattle::nhc_temp_integrate() +{ + int ich; + // t_current = temperature->compute_scalar(); + // tdof = temperature->dof; + compute_temp_target(); + + double expfac, kecurrent = tdof * boltz * t_current; + double inv_t_freq2 = 1.0 / (t_freq*t_freq); + eta_mass[0] = tdof * boltz * t_target * inv_t_freq2; + for( int ich = 1; ich < mtchain; ++ich ){ + eta_mass[ich] = boltz * t_target * inv_t_freq2; + } + + if( eta_mass[0] > 0.0 ){ + eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0]; + }else{ + eta_dotdot[0] = 0; + } + + for( ich = mtchain-1; ich > 0; --ich ){ + expfac = exp(-dt8*eta_dot[ich+1]); + eta_dot[ich] *= expfac; + eta_dot[ich] += eta_dotdot[ich] * dt4; + eta_dot[ich] *= tdrag_factor * expfac; + + } + + expfac = exp(-dt8*eta_dot[1]); + eta_dot[0] *= expfac; + eta_dot[0] += eta_dotdot[0] * dt4; + eta_dot[0] *= tdrag_factor * expfac; + + factor_eta = exp(-dthalf*eta_dot[0]); + + if( factor_eta == 0 ){ + char msg[2048]; + sprintf(msg, "WTF, factor_eta is 0! dthalf = %f, eta_dot[0] = %f", + dthalf, eta_dot[0]); + error->all(FLERR,msg); + } + + nh_v_temp(); + + t_current *= factor_eta*factor_eta; + kecurrent = tdof * boltz * t_current; + + if( eta_mass[0] > 0.0 ){ + eta_dotdot[0] = (kecurrent - ke_target) / eta_mass[0]; + }else{ + eta_dotdot[0] = 0.0; + } + + for( int ich = 1; ich < mtchain; ++ich ){ + eta[ich] += dthalf*eta_dot[ich]; + } + eta_dot[0] *= expfac; + eta_dot[0] += eta_dotdot[0]*dt4; + eta_dot[0] *= expfac; + + for( int ich = 1; ich < mtchain; ++ich ){ + expfac = exp(-dt8*eta_dot[ich+1]); + eta_dot[ich] *= expfac; + eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] + - boltz*t_target) / eta_mass[ich]; + eta_dot[ich] *= eta_dotdot[ich] * dt4; + eta_dot[ich] *= expfac; + } + +} + +void FixNVTManifoldRattle::nh_v_temp() +{ + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if( igroup == atom->firstgroup) nlocal = atom->nfirst; + + + + + if( which == NOBIAS ){ + for( int i = 0; i < nlocal; ++i ){ + if( mask[i] & groupbit ){ + v[i][0] *= factor_eta; + v[i][1] *= factor_eta; + v[i][2] *= factor_eta; + } + } + }else if( which == BIAS ){ + for( int i = 0; i < nlocal; ++i ){ + if( mask[i] & groupbit ){ + temperature->remove_bias(i,v[i]); + v[i][0] *= factor_eta; + v[i][1] *= factor_eta; + v[i][2] *= factor_eta; + temperature->restore_bias(i,v[i]); + } + } + } +} + + + + +// Most of this logic is based on fix_nh: +void FixNVTManifoldRattle::initial_integrate(int vflag) +{ + + update_var_params(); + + compute_temp_target(); + nhc_temp_integrate(); + + nve_x_rattle(igroup, groupbit); +} + +void FixNVTManifoldRattle::final_integrate() +{ + nve_v_rattle(igroup, groupbit); + + nhc_temp_integrate(); +} + + + +/* ---------------------------------------------------------------------- */ +void FixNVTManifoldRattle::reset_dt() +{ + FixNVEManifoldRattle::reset_dt(); + + dthalf = 0.5 * update->dt; + dt4 = 0.25 * update->dt; + dt8 = 0.125 * update->dt; + tdrag_factor = 1.0 - (update->dt * t_freq * drag); + +} + + + + + +double FixNVTManifoldRattle::memory_usage() +{ + double bytes = FixNVEManifoldRattle::memory_usage(); + bytes += (4*mtchain+1)*sizeof(double); + + return bytes; +} diff --git a/src/USER-MANIFOLD/fix_nvt_manifold_rattle.h b/src/USER-MANIFOLD/fix_nvt_manifold_rattle.h new file mode 100644 index 0000000000..a9d3e3122f --- /dev/null +++ b/src/USER-MANIFOLD/fix_nvt_manifold_rattle.h @@ -0,0 +1,155 @@ +/* ---------------------------------------------------------------------- + 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. + + Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology. + License: GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + This file is part of the user-manifold package written by + Stefan Paquay at the Eindhoven University of Technology. + This module makes it possible to do MD with particles constrained + to pretty arbitrary manifolds characterised by some constraint function + g(x,y,z) = 0 and its normal grad(g). The number of manifolds available + right now is limited but can be extended straightforwardly by making + a new class that inherits from manifold and implements all pure virtual + methods. + + Thanks to Remy Kusters for beta-testing! + +------------------------------------------------------------------------- */ + + +#ifdef FIX_CLASS + +FixStyle(nvt/manifold/rattle,FixNVTManifoldRattle) + +#else + +#ifndef LMP_FIX_NVT_MANIFOLD_RATTLE_H +#define LMP_FIX_NVT_MANIFOLD_RATTLE_H + +#include "fix_nve_manifold_rattle.h" + + + +/* + FixNVTManifoldRattle works by wrapping some Nose-Hoover thermostat routines + around the time integration functions of FixNVEManifoldRattle. +*/ +namespace LAMMPS_NS { + +// namespace user_manifold { + + class FixNVTManifoldRattle : public FixNVEManifoldRattle { + public: + FixNVTManifoldRattle(LAMMPS *, int, char **, int = 1); + virtual ~FixNVTManifoldRattle(); + + virtual void initial_integrate(int); + virtual void final_integrate(); + virtual void init(); + virtual void reset_dt(); + virtual int setmask(); + virtual void setup(int); // Not needed for fixNVE but is for fixNVT + virtual double memory_usage(); + + + + protected: + + void compute_temp_target(); + void nhc_temp_integrate(); + void nh_v_temp(); + + double dthalf, dt4, dt8; + + char *id_temp; + class Compute* temperature; + double t_start,t_stop, t_period; + double t_current,t_target,ke_target; + double t_freq, drag, tdrag_factor; + double boltz,nktv2p,tdof; + double *eta,*eta_dot; + double *eta_dotdot; + double *eta_mass; + int mtchain; + double factor_eta; + int which, got_temp; + + const char *fix_id; + }; +} + + + + +#endif // LMP_FIX_NVE_MANIFOLD_RATTLE_H +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: There is no manifold named ... + +Self-explanatory. You requested a manifold whose name was not +registered at the factory. + +E: Manifold pointer was NULL for some reason! + +This indicates a bug. The factory was unable to properly create +the requested manifold even though it was registered. Send the +maintainer an e-mail. + +E: Manifold ... needs at least ... argument(s)! + +Self-explanatory. Provide enough arguments for the proper +creating of the requested manifold. + +E: Parameter pointer was NULL! + +This indicates a bug. The array that contains the parameters +could not be allocated. Send the maintainer an e-mail. + +E: Could not allocate space for arg! + +One of the arguments provided was too long (it contained +too many characters) + +E: Option ... needs ... argument(s): + +Self-explanatory. Read the documentation of this fix properly. + + +E: Illegal fix nve/manifold/rattle command! Option ... not recognized! + +Self-explanatory. The given option(s) do not exist. + +E: Variable name for fix nve/manifold/rattle does not exist + +Self-explanatory. + +E: Variable for fix nve/manifold/rattle is invalid style + +fix nve/manifold/rattle only supports equal style variables. + + + +*/ diff --git a/src/USER-MANIFOLD/manifold.h b/src/USER-MANIFOLD/manifold.h new file mode 100644 index 0000000000..82a565af38 --- /dev/null +++ b/src/USER-MANIFOLD/manifold.h @@ -0,0 +1,100 @@ +/* ---------------------------------------------------------------------- + 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, s.paquay@tue.nl + 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_H +#define LMP_MANIFOLD_H + +#include "pointers.h" +#include + +namespace LAMMPS_NS { +namespace user_manifold { + + // Abstract base class. + class manifold : protected Pointers { + public: + manifold(class LAMMPS* lmp) : Pointers(lmp), params(NULL){ } + virtual ~manifold(){ delete[] params; } + virtual double g( const double * ) = 0; + virtual void n( const double *, double * ) = 0; + + + // Variant of g that computes n at the same time. + virtual double g_and_n( const double *x, double *nn ) + { + this->n(x,nn); + return g(x); + } + + virtual const char *id() = 0; + + virtual void set_atom_id( tagint a_id ){} + virtual int nparams() = 0; + double **get_params(){ return ¶ms; }; + + // Overload if any initialisation depends on params: + virtual void post_param_init(){} + virtual void checkup(){} // Some diagnostics... + protected: + double *params; + }; + + + + // Some utility functions that are templated, so I implement them + // here in the header. + template< unsigned int size > inline + double infnorm( double *vect ) + { + double largest = fabs( vect[0] ); + for( unsigned int i = 1; i < size; ++i ){ + double c = fabs( vect[i] ); + largest = ( c > largest ) ? c : largest; + } + return largest; + } + + inline double dot( double *a, double *b ){ + return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; + } + + +} // namespace LAMMPS_NS + +} + + +#endif // LMP_MANIFOLD_H diff --git a/src/USER-MANIFOLD/manifold_cylinder.cpp b/src/USER-MANIFOLD/manifold_cylinder.cpp new file mode 100644 index 0000000000..c13581c640 --- /dev/null +++ b/src/USER-MANIFOLD/manifold_cylinder.cpp @@ -0,0 +1,27 @@ +#include "manifold_cylinder.h" + +using namespace LAMMPS_NS; + +using namespace user_manifold; + + +manifold_cylinder::manifold_cylinder( LAMMPS *lmp, int argc, + char **argv ) : manifold(lmp) +{} + + +double manifold_cylinder::g( const double *x ) +{ + double R = params[0]; + double r2 = x[0]*x[0] + x[1]*x[1]; + return R*R - r2; +} + +void manifold_cylinder::n( const double *x, double *n ) +{ + n[0] = -2*x[0]; + n[1] = -2*x[1]; + n[2] = 0.0; +} + + diff --git a/src/USER-MANIFOLD/manifold_cylinder.h b/src/USER-MANIFOLD/manifold_cylinder.h new file mode 100644 index 0000000000..f8cddadf66 --- /dev/null +++ b/src/USER-MANIFOLD/manifold_cylinder.h @@ -0,0 +1,29 @@ +#ifndef LMP_MANIFOLD_CYLINDER_H +#define LMP_MANIFOLD_CYLINDER_H + +#include "manifold.h" + +// A normal cylinder + + +namespace LAMMPS_NS { + +namespace user_manifold { + + class manifold_cylinder : public manifold { + public: + enum { NPARAMS = 1 }; // Number of parameters. + manifold_cylinder( LAMMPS *lmp, int, char ** ); + virtual ~manifold_cylinder(){} + virtual double g( const double *x ); + virtual void n( const double *x, double *n ); + static const char *type(){ return "cylinder"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + +} + +#endif // LMP_MANIFOLD_CYLINDER_H diff --git a/src/USER-MANIFOLD/manifold_cylinder_dent.cpp b/src/USER-MANIFOLD/manifold_cylinder_dent.cpp new file mode 100644 index 0000000000..687e311c9c --- /dev/null +++ b/src/USER-MANIFOLD/manifold_cylinder_dent.cpp @@ -0,0 +1,44 @@ +#include "manifold_cylinder_dent.h" +#include "math_const.h" + +#include + + +using namespace LAMMPS_NS; + +using namespace user_manifold; + +manifold_cylinder_dent::manifold_cylinder_dent( LAMMPS *lmp, int argc, + char **argv ) : manifold(lmp) +{} + + +double manifold_cylinder_dent::g( const double *x ) +{ + double l = params[1], R = params[0], a = params[2]; + double r2 = x[1]*x[1] + x[0]*x[0]; + if( fabs(x[2]) < 0.5*l ){ + double k = MathConst::MY_2PI / l; + double c = R - 0.5*a*( 1.0 + cos(k*x[2]) ); + return c*c - r2; + }else{ + return R*R - r2; + } +} + + +void manifold_cylinder_dent::n( const double *x, double *n ) +{ + double l = params[1], R = params[0], a = params[2]; + if( fabs(x[2]) < 0.5*l ){ + double k = MathConst::MY_2PI / l; + double c = R - 0.5*a*(1.0 + cos(k*x[2])); + n[0] = -2*x[0]; + n[1] = -2*x[1]; + n[2] = c*a*k*sin(k*x[2]); + }else{ + n[0] = -2*x[0]; + n[1] = -2*x[1]; + n[2] = 0.0; + } +} diff --git a/src/USER-MANIFOLD/manifold_cylinder_dent.h b/src/USER-MANIFOLD/manifold_cylinder_dent.h new file mode 100644 index 0000000000..33167d670b --- /dev/null +++ b/src/USER-MANIFOLD/manifold_cylinder_dent.h @@ -0,0 +1,27 @@ +#ifndef LMP_MANIFOLD_CYLINDER_DENT_H +#define LMP_MANIFOLD_CYLINDER_DENT_H + +#include "manifold.h" + + +namespace LAMMPS_NS { + +namespace user_manifold { + + class manifold_cylinder_dent : public manifold { + public: + manifold_cylinder_dent( LAMMPS *lmp, int, char ** ); + enum { NPARAMS = 3 }; // Number of parameters. + virtual ~manifold_cylinder_dent(){} + virtual double g( const double *x ); + virtual void n( const double *x, double *n ); + static const char *type(){ return "cylinder/dent"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + +} + +#endif // LMP_MANIFOLD_CYLINDER_DENT_H diff --git a/src/USER-MANIFOLD/manifold_dumbbell.cpp b/src/USER-MANIFOLD/manifold_dumbbell.cpp new file mode 100644 index 0000000000..ed1583ac9a --- /dev/null +++ b/src/USER-MANIFOLD/manifold_dumbbell.cpp @@ -0,0 +1,43 @@ +#include "manifold_dumbbell.h" + +#include + +using namespace LAMMPS_NS; + +using namespace user_manifold; + +manifold_dumbbell::manifold_dumbbell( LAMMPS *lmp, int argc, char **argv ) : manifold(lmp) +{} + + +/* + * Equation for dumbbell is: + * + * -x^2 - y^2 + (a^2 - (z/c)^2)*( 1 + A*sin(const(z) *z^2) )^4 + * const(z) = (z < 0) ? B2 : B + * params[4] = { a, A, B, c } + */ + + +double manifold_dumbbell::g( const double *x ) +{ + double a = params[0]; + double A = params[1]; + double B = params[2]; + double c = params[3]; + return -1.0*(x[0]*x[0]+x[1]*x[1])+(a*a-x[2]*x[2]/(c*c))*(1.0+pow(A*sin(B *x[2]*x[2]),4)); +} + +void manifold_dumbbell::n( const double *x, double *nn ) +{ + double a = params[0]; + double A = params[1]; + double B = params[2]; + double c = params[3]; + nn[0] = -2.0*x[0]; + nn[1] = -2.0*x[1]; + nn[2] = 8.0*A*A*A*A*B*x[2]*(a*a-(x[2]*x[2]/(c*c)))*cos(B*x[2]*x[2]) * + pow(sin(B*x[2]*x[2]),3)-2*x[2]*(1+pow(A*sin(B*x[2]*x[2]),4))/(c*c); +} + + diff --git a/src/USER-MANIFOLD/manifold_dumbbell.h b/src/USER-MANIFOLD/manifold_dumbbell.h new file mode 100644 index 0000000000..fea718b487 --- /dev/null +++ b/src/USER-MANIFOLD/manifold_dumbbell.h @@ -0,0 +1,31 @@ +#ifndef LMP_MANIFOLD_DUMBBELL_H +#define LMP_MANIFOLD_DUMBBELL_H + +#include "manifold.h" + + + +namespace LAMMPS_NS { + +namespace user_manifold { + + // A dendritic dumbbell approximation: + class manifold_dumbbell : public manifold { + public: + enum { NPARAMS = 4 }; // Number of parameters. + manifold_dumbbell( LAMMPS *lmp, int, char ** ); + virtual ~manifold_dumbbell(){} + virtual double g ( const double *x ); + virtual void n ( const double *x, double *nn ); + + static const char* type(){ return "dumbbell"; } + virtual const char *id(){ return type(); } + + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + +} + +#endif // LMP_MANIFOLD_DUMBBELL_H diff --git a/src/USER-MANIFOLD/manifold_ellipsoid.cpp b/src/USER-MANIFOLD/manifold_ellipsoid.cpp new file mode 100644 index 0000000000..f590123dc3 --- /dev/null +++ b/src/USER-MANIFOLD/manifold_ellipsoid.cpp @@ -0,0 +1,28 @@ +#include "manifold_ellipsoid.h" + +using namespace LAMMPS_NS; + +using namespace user_manifold; + +manifold_ellipsoid::manifold_ellipsoid( LAMMPS *lmp, int narg, char **argv ) : manifold(lmp) +{} + + +double manifold_ellipsoid::g( const double *x ) +{ + const double ai = 1.0 / params[0]; + const double bi = 1.0 / params[1]; + const double ci = 1.0 / params[2]; + return x[0]*x[0]*ai*ai + x[1]*x[1]*bi*bi + x[2]*x[2]*ci*ci - 1.0; +} + +void manifold_ellipsoid::n( const double *x, double * n ) +{ + const double ai = 1.0 / params[0]; + const double bi = 1.0 / params[1]; + const double ci = 1.0 / params[2]; + n[0] = 2*x[0]*ai*ai; + n[1] = 2*x[1]*bi*bi; + n[2] = 2*x[2]*ci*ci; +} + diff --git a/src/USER-MANIFOLD/manifold_ellipsoid.h b/src/USER-MANIFOLD/manifold_ellipsoid.h new file mode 100644 index 0000000000..e23dbec72a --- /dev/null +++ b/src/USER-MANIFOLD/manifold_ellipsoid.h @@ -0,0 +1,28 @@ +#ifndef LMP_MANIFOLD_ELLIPSOID_H +#define LMP_MANIFOLD_ELLIPSOID_H + +#include "manifold.h" + + +namespace LAMMPS_NS { + +namespace user_manifold { + // An ellipsoid: + class manifold_ellipsoid : public manifold { + public: + enum { NPARAMS = 3 }; + manifold_ellipsoid( LAMMPS *lmp, int, char ** ); + virtual ~manifold_ellipsoid(){} + virtual double g( const double *x ); + virtual void n( const double *x, double *n ); + + static const char* type(){ return "ellipsoid"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + +} + +#endif // LMP_MANIFOLD_ELLIPSOID_H diff --git a/src/USER-MANIFOLD/manifold_factory.cpp b/src/USER-MANIFOLD/manifold_factory.cpp new file mode 100644 index 0000000000..57c6e8305f --- /dev/null +++ b/src/USER-MANIFOLD/manifold_factory.cpp @@ -0,0 +1,73 @@ +/* ---------------------------------------------------------------------- + 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. + + Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology. + License: GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + This file is part of the user-manifold package written by + Stefan Paquay at the Eindhoven University of Technology. + This module makes it possible to do MD with particles constrained + to pretty arbitrary manifolds characterised by some constraint function + g(x,y,z) = 0 and its normal grad(g). The number of manifolds available + right now is limited but can be extended straightforwardly by making + a new class that inherits from manifold and implements all pure virtual + methods. + + Thanks to Remy Kusters for beta-testing! + +------------------------------------------------------------------------- */ + +#include "manifold_factory.h" + +#include "manifold_cylinder.h" +#include "manifold_cylinder_dent.h" +#include "manifold_dumbbell.h" +#include "manifold_ellipsoid.h" +#include "manifold_plane.h" +#include "manifold_plane_wiggle.h" +#include "manifold_sphere.h" +#include "manifold_supersphere.h" +#include "manifold_spine.h" +#include "manifold_thylakoid.h" +#include "manifold_torus.h" + +using namespace LAMMPS_NS; +using namespace user_manifold; + + +manifold* LAMMPS_NS::user_manifold::create_manifold(const char *mname, + LAMMPS *lmp, + int narg, char **arg ) +{ + manifold *man = NULL; + 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 ); + + + return man; +} + diff --git a/src/USER-MANIFOLD/manifold_factory.h b/src/USER-MANIFOLD/manifold_factory.h new file mode 100644 index 0000000000..a2f21861eb --- /dev/null +++ b/src/USER-MANIFOLD/manifold_factory.h @@ -0,0 +1,110 @@ +/* ---------------------------------------------------------------------- + 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. + + Copyright (2013-2014) Stefan Paquay, Eindhoven University of Technology. + License: GNU General Public License. + + See the README file in the top-level LAMMPS directory. + + This file is part of the user-manifold package written by + Stefan Paquay at the Eindhoven University of Technology. + This module makes it possible to do MD with particles constrained + to pretty arbitrary manifolds characterised by some constraint function + g(x,y,z) = 0 and its normal grad(g). The number of manifolds available + right now is limited but can be extended straightforwardly by making + a new class that inherits from manifold and implements all pure virtual + methods. + + Thanks to Remy Kusters for beta-testing! + +------------------------------------------------------------------------- */ + + +#ifndef LMP_MANIFOLD_FACTORY_H +#define LMP_MANIFOLD_FACTORY_H + + +#include "manifold.h" +#include + +/* + * Defining USE_PHONY_LAMMPS makes sure that none of the LAMMPS classes are + * included/compiled. This is done in order to allow other programs to use + * the manifold_factory without compiling all of LAMMPS itself. The relevant + * classes/functions are replaced with dummy ones defined in this #ifdef-block: + */ +#ifdef USE_PHONY_LAMMPS +# ifdef __GNUG__ +# warning Not compiling actual LAMMPS classes! +# endif + + +struct Error { + void all(const char *fname, int line, const char* msg) + { + fprintf(stderr,"ERROR: %s (%s:%d)",msg,fname,line); + std::terminate(); + } + + void one(const char *fname, int line, const char* msg) + { + fprintf(stderr,"ERROR: %s (%s:%d)",msg,fname,line); + std::terminate(); + } +}; + +struct LAMMPS { }; + +struct Pointers +{ + Pointers(LAMMPS *) : error( &e ){} + Error e; + Error *error; +}; + +static FILE *screen = fopen("/dev/stdout","w"); + +#define FLERR __FILE__,__LINE__ // Equivalent to definition in pointers.h +#endif // USE_PHONY_LAMMPS + + + +/* Here the actual implementation of LAMMPS-related functions begins. */ + +namespace LAMMPS_NS { + +namespace user_manifold { + + // Templated, so needs to be in header. + template + void make_manifold_if( manifold **man_ptr, const char *name, + LAMMPS *lmp, int narg, char **arg ) + { + if( strcmp( m_type::type(), name ) == 0 ){ + if( *man_ptr == NULL ){ + *man_ptr = new m_type(lmp, narg, arg); + } + } + } + + manifold* create_manifold(const char *, LAMMPS *, + int , char ** ); + +} // namespace user_manifold + +} // namespace LAMMPS_NS + + +#endif // LMP_MANIFOLD_FACTORY_H diff --git a/src/USER-MANIFOLD/manifold_plane.cpp b/src/USER-MANIFOLD/manifold_plane.cpp new file mode 100644 index 0000000000..d33617fbfb --- /dev/null +++ b/src/USER-MANIFOLD/manifold_plane.cpp @@ -0,0 +1,24 @@ +#include "manifold_plane.h" + +using namespace LAMMPS_NS; + +using namespace user_manifold; + +manifold_plane::manifold_plane( LAMMPS *lmp, int argc, char **argv ) : + manifold(lmp) +{} + +double manifold_plane::g( const double *x ) +{ + double a = params[0], b = params[1], c = params[2]; + double x0 = params[3], y0 = params[4], z0 = params[5]; + return a*(x[0] - x0) + b*(x[1] - y0) + c*(x[2] - z0); +} + + +void manifold_plane::n( const double *x, double *n ) +{ + n[0] = params[0]; + n[1] = params[1]; + n[2] = params[2]; +} diff --git a/src/USER-MANIFOLD/manifold_plane.h b/src/USER-MANIFOLD/manifold_plane.h new file mode 100644 index 0000000000..738c632b17 --- /dev/null +++ b/src/USER-MANIFOLD/manifold_plane.h @@ -0,0 +1,30 @@ +#ifndef LMP_MANIFOLD_PLANE_H +#define LMP_MANIFOLD_PLANE_H + +#include "manifold.h" + + +namespace LAMMPS_NS { + +namespace user_manifold { + + + // A 2D plane + class manifold_plane : public manifold { + public: + enum { NPARAMS = 6 }; // Number of parameters. + manifold_plane( LAMMPS *lmp, int, char ** ); + virtual ~manifold_plane(){} + virtual double g( const double *x ); + virtual void n( const double *x, double *n ); + static const char *type(){ return "plane"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + + +} + +#endif // LMP_MANIFOLD_PLANE_H diff --git a/src/USER-MANIFOLD/manifold_plane_wiggle.cpp b/src/USER-MANIFOLD/manifold_plane_wiggle.cpp new file mode 100644 index 0000000000..fd50e774da --- /dev/null +++ b/src/USER-MANIFOLD/manifold_plane_wiggle.cpp @@ -0,0 +1,28 @@ +#include "manifold_plane_wiggle.h" + +#include + +using namespace LAMMPS_NS; +using namespace user_manifold; + +manifold_plane_wiggle::manifold_plane_wiggle( LAMMPS *lmp, int argc, char **argv ) : + manifold(lmp) +{} + + +double manifold_plane_wiggle::g( const double *x ) +{ + double a = params[0]; + double w = params[1]; + return x[2] - a*sin(w*x[0]); +} + + +void manifold_plane_wiggle::n( const double *x, double *n ) +{ + double a = params[0]; + double w = params[1]; + n[2] = 1; + n[1] = 0.0; + n[0] = -a*w*cos(x[0]); +} diff --git a/src/USER-MANIFOLD/manifold_plane_wiggle.h b/src/USER-MANIFOLD/manifold_plane_wiggle.h new file mode 100644 index 0000000000..fe84713d7d --- /dev/null +++ b/src/USER-MANIFOLD/manifold_plane_wiggle.h @@ -0,0 +1,28 @@ +#ifndef LMP_MANIFOLD_PLANE_WIGGLE_H +#define LMP_MANIFOLD_PLANE_WIGGLE_H + +#include "manifold.h" + + +namespace LAMMPS_NS { + +namespace user_manifold { + + // A 2D wiggly/wave-y plane (Like z = A cos(kx)) + class manifold_plane_wiggle : public manifold { + public: + enum { NPARAMS = 2 }; // Number of parameters. + manifold_plane_wiggle( LAMMPS *lmp, int, char ** ); + virtual ~manifold_plane_wiggle(){} + virtual double g( const double *x ); + virtual void n( const double *x, double *n ); + static const char *type(){ return "plane/wiggle"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + +} + +#endif // LMP_MANIFOLD_PLANE_WIGGLE_H diff --git a/src/USER-MANIFOLD/manifold_sphere.h b/src/USER-MANIFOLD/manifold_sphere.h new file mode 100644 index 0000000000..5d4ae55fda --- /dev/null +++ b/src/USER-MANIFOLD/manifold_sphere.h @@ -0,0 +1,59 @@ +#ifndef LMP_MANIFOLD_SPHERE_H +#define LMP_MANIFOLD_SPHERE_H + +#include "manifold.h" + +namespace LAMMPS_NS { + +namespace user_manifold { + + + // A sphere: + class manifold_sphere : public manifold { + public: + enum { NPARAMS = 1 }; + manifold_sphere( LAMMPS *lmp, int, char ** ) : manifold(lmp){} + + virtual ~manifold_sphere(){} + virtual double g( const double *x ) + { + double R = params[0]; + double r2 = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]; + return r2 - R*R; + } + + virtual double g_and_n( const double *x, double *nn ) + { + double R = params[0]; + double r2 = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]; + nn[0] = 2*x[0]; + nn[1] = 2*x[1]; + nn[2] = 2*x[2]; + + return r2 - R*R; + } + + virtual void n( const double *x, double *nn ) + { + nn[0] = 2*x[0]; + nn[1] = 2*x[1]; + nn[2] = 2*x[2]; + } + + virtual void H( double *x, double h[3][3] ) + { + h[0][1] = h[0][2] = h[1][0] = h[1][2] = h[2][0] = h[2][1] = 0.0; + h[0][0] = h[1][1] = h[2][2] = 2.0; + } + + static const char* type(){ return "sphere"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + +} + + +#endif // LMP_MANIFOLD_SPHERE_H diff --git a/src/USER-MANIFOLD/manifold_spine.cpp b/src/USER-MANIFOLD/manifold_spine.cpp new file mode 100644 index 0000000000..1081ab673e --- /dev/null +++ b/src/USER-MANIFOLD/manifold_spine.cpp @@ -0,0 +1,163 @@ +#include "manifold_spine.h" + +#include +#include "math_special.h" + +using namespace LAMMPS_NS; +using namespace user_manifold; + + + +manifold_spine::manifold_spine( LAMMPS *lmp, int argc, char **argv ) + : manifold(lmp) +{ + power = 4; +} + + + +manifold_spine_two::manifold_spine_two( LAMMPS *lmp, int narg, char **argv ) + : manifold_spine( lmp, narg, argv ) +{ + power = 2; +} + + + +/* + * Equation for spine is: + * + * -x^2 - y^2 + (a^2 - (z/c)^2)*( 1 + (A*sin(const(z) *z^2))^4 ) + * const(z) = (z < 0) ? B2 : B + * params[5] = { a, A, B, B2, c } + */ + +double manifold_spine::g_and_n( const double *x, double *nn ) +{ + double a = params[0]; + double A = params[1]; + double B = params[2]; + double B2 = params[3]; + double c = params[4]; + + double cc, BB, z2, xy2; + double c2, As, Ac, azc, Apart; + double AMs, AMc; + double dazc, dAMs; + + if( x[2] > 0 ){ + BB = B; + cc = c; + }else{ + BB = B2; + cc = 1.0; + } + + xy2 = x[0]*x[0] + x[1]*x[1]; + z2 = x[2]*x[2]; + c2 = 1.0/( cc*cc ); + + azc = a*a - z2*c2; + As = sin( BB*z2 ); + Ac = cos( BB*z2 ); + As *= A; + Ac *= A; + + Apart = MathSpecial::powint( As, power-1 ); // Much more efficient! =D + //Apart = pow( As, power-1 ); + AMs = Apart*As; + AMc = Apart*Ac; + + dAMs = power * AMc * 2.0*BB*x[2]; + dazc = -2.0*x[2]*c2; + + nn[0] = -2*x[0]; + nn[1] = -2*x[1]; + nn[2] = azc * dAMs + ( 1.0 + AMs ) * dazc; + + return -xy2 + azc * ( 1.0 + AMs ); +} + + + + +void manifold_spine::n( const double *x, double *nn ) +{ + double a = params[0]; + double A = params[1]; + double B = params[2]; + double B2 = params[3]; + double c = params[4]; + + double cc, BB, z2; + double c2, As, Ac, azc, Apart; + double AMs, AMc; + double dazc, dAMs; + + if( x[2] > 0 ){ + BB = B; + cc = c; + }else{ + BB = B2; + cc = 1.0; + } + + z2 = x[2]*x[2]; + c2 = 1.0/( cc*cc ); + + azc = a*a - z2*c2; + As = sin( BB*z2 ); + Ac = cos( BB*z2 ); + As *= A; + Ac *= A; + + Apart = MathSpecial::powint( As, power-1 ); // Much more efficient! =D + //Apart = pow( As, power-1 ); + AMs = Apart*As; + AMc = Apart*Ac; + + dAMs = power * AMc * 2.0*BB*x[2]; + dazc = -2.0*x[2]*c2; + + nn[0] = -2*x[0]; + nn[1] = -2*x[1]; + nn[2] = azc * dAMs + ( 1.0 + AMs ) * dazc; + +} + + +double manifold_spine::g( const double *x ) +{ + double a = params[0]; + double A = params[1]; + double B = params[2]; + double B2 = params[3]; + double c = params[4]; + + double cc, BB, z2, xy2; + double c2, As, azc, Apart; + double AMs; + + if( x[2] > 0 ){ + BB = B; + cc = c; + }else{ + BB = B2; + cc = 1.0; + } + + xy2 = x[0]*x[0] + x[1]*x[1]; + z2 = x[2]*x[2]; + c2 = 1.0/( cc*cc ); + + azc = a*a - z2*c2; + As = sin( BB*z2 ); + As *= A; + + Apart = MathSpecial::powint( As, power-1 ); + AMs = Apart*As; + return -xy2 + azc * ( 1.0 + AMs ); + +} + + diff --git a/src/USER-MANIFOLD/manifold_spine.h b/src/USER-MANIFOLD/manifold_spine.h new file mode 100644 index 0000000000..ba5e82a6bb --- /dev/null +++ b/src/USER-MANIFOLD/manifold_spine.h @@ -0,0 +1,43 @@ +#ifndef LMP_MANIFOLD_SPINE_H +#define LMP_MANIFOLD_SPINE_H + +#include "manifold.h" + + +namespace LAMMPS_NS { + +namespace user_manifold { + + // A dendritic spine approximation: + class manifold_spine : public manifold { + public: + enum { NPARAMS = 5 }; // Number of parameters. + manifold_spine( LAMMPS *lmp, int, char ** ); + virtual ~manifold_spine(){} + virtual double g ( const double *x ); + virtual void n ( const double *x, double *nn ); + virtual double g_and_n( const double *x, double *nn ); + + static const char* type(){ return "spine"; } + virtual const char *id(){ return type(); } + + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + protected: + int power; + }; + + class manifold_spine_two : public manifold_spine { + public: + manifold_spine_two( LAMMPS *lmp, int, char **); + + static const char* type(){ return "spine/two"; } + virtual const char *id(){ return type(); } + + }; +} + +} + + +#endif // LMP_MANIFOLD_SPINE_H diff --git a/src/USER-MANIFOLD/manifold_supersphere.h b/src/USER-MANIFOLD/manifold_supersphere.h new file mode 100644 index 0000000000..1d9c62d7a2 --- /dev/null +++ b/src/USER-MANIFOLD/manifold_supersphere.h @@ -0,0 +1,58 @@ +#ifndef LMP_MANIFOLD_SUPERSPHERE_H +#define LMP_MANIFOLD_SUPERSPHERE_H + +#include "manifold.h" + +namespace LAMMPS_NS { + + +namespace user_manifold { + + // A sphere: + class manifold_supersphere : public manifold { + public: + enum { NPARAMS = 2 }; + manifold_supersphere( LAMMPS *lmp, int, char ** ) : manifold(lmp){} + + virtual ~manifold_supersphere(){} + + double my_sign( double a ) + { + return (a > 0) - (a < 0); + } + + virtual double g( const double *x ) + { + double R = params[0]; + double q = params[1]; + double xx = fabs(x[0]); + double yy = fabs(x[1]); + double zz = fabs(x[2]); + + double rr = pow(xx,q) + pow(yy,q) + pow(zz,q); + + return rr - pow(R,q); + } + + virtual void n( const double *x, double *nn ) + { + double q = params[1]; + double xx = fabs(x[0]); + double yy = fabs(x[1]); + double zz = fabs(x[2]); + + nn[0] = q * my_sign(x[0])*pow(xx,q-1); + nn[1] = q * my_sign(x[1])*pow(yy,q-1); + nn[2] = q * my_sign(x[2])*pow(zz,q-1); + } + + static const char* type(){ return "supersphere"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; +} + +} + +#endif diff --git a/src/USER-MANIFOLD/manifold_thylakoid.cpp b/src/USER-MANIFOLD/manifold_thylakoid.cpp new file mode 100644 index 0000000000..719a2d3281 --- /dev/null +++ b/src/USER-MANIFOLD/manifold_thylakoid.cpp @@ -0,0 +1,618 @@ +#include "manifold_thylakoid.h" +#include + +#include "comm.h" +#include "domain.h" // For some checks regarding the simulation box. +#include "error.h" + + +#define MANIFOLD_THYLAKOID_DEBUG + +using namespace LAMMPS_NS; +using namespace user_manifold; + + +manifold_thylakoid::manifold_thylakoid( LAMMPS *lmp, int narg, char ** arg) + : manifold(lmp) +{ + // You can NOT depend on proper construction of the domains in + // the constructor, because the params are set by the function that + // calls the factory-construction method. Instead, the constructing + // fix should call post_param_init(); +} + + + +manifold_thylakoid::~manifold_thylakoid() +{ + for( std::size_t i = 0; i < parts.size(); ++i ){ + delete parts[i]; + } +} + + + +void manifold_thylakoid::post_param_init() +{ + // Set coefficients: + lT = 3; // Radius of cylinder edges of stacks + LT = 15; // Size of faces of cylinder stacks + pad = 3; // Padding to prevent improper look-ups + wB = 3.0; + LB = 10.0; + lB = 3.0; + + wB = params[0]; + LB = params[1]; + lB = params[2]; + + if( comm->me == 0 ){ + fprintf(screen,"My params are now: lT = %f, LT = %f, pad = %f, " + "wB = %f, LB = %f, lB = %f\n", lT, LT, pad, wB, LB, lB ); + fprintf(screen,"Calling init_domains() from post_param_init().\n"); + } + init_domains(); + checkup(); +} + +void manifold_thylakoid::checkup() +{ + if( comm->me == 0 ){ + fprintf(screen,"This is checkup of thylakoid %p\n", this); + fprintf(screen,"I have %ld parts. They are:\n", parts.size()); + for( int i = 0; i < parts.size(); ++i ){ + fprintf(screen, "[%f, %f] x [%f, %f] x [%f, %f]\n", + parts[i]->xlo, parts[i]->xhi, + parts[i]->ylo, parts[i]->yhi, + parts[i]->zlo, parts[i]->zhi ); + } + fprintf(screen,"My params are:\n"); + for( int i = 0; i < NPARAMS; ++i ){ + fprintf(screen,"%f\n", params[i]); + } + } +} + + +double manifold_thylakoid::g( const double *x ) +{ + int err = 0; + std::size_t idx; + thyla_part *p = get_thyla_part(x,&err,&idx); + if(err){ + char msg[2048]; + sprintf(msg,"Error getting thyla_part for x = (%f, %f, %f)",x[0],x[1],x[2]); + error->one(FLERR,msg); + } + double con_val = p->g(x); + if( isfinite(con_val) ){ + return con_val; + }else{ + char msg[2048]; + sprintf(msg,"Error, thyla_part of type %d returned %f as constraint val!", + p->type, con_val); + error->one(FLERR,msg); + return 0; + } +} + +void manifold_thylakoid::n( const double *x, double *n ) +{ + int err = 0; + std::size_t idx; + thyla_part *p = get_thyla_part(x,&err,&idx); + if(err){ + char msg[2048]; + sprintf(msg,"Error getting thyla_part for x = (%f, %f, %f)",x[0],x[1],x[2]); + error->one(FLERR,msg); + } + p->n(x,n); + if( isfinite(n[0]) && isfinite(n[1]) && isfinite(n[2]) ){ + return; + }else{ + char msg[2048]; + sprintf(msg,"Error, thyla_part of type %d returned (%f,%f,%f) as gradient!", + p->type, n[0], n[1], n[2]); + error->one(FLERR,msg); + } +} + +thyla_part *manifold_thylakoid::get_thyla_part( const double *x, int *err_flag, std::size_t *idx ) +{ + + for( std::size_t i = 0; i < parts.size(); ++i ){ + thyla_part *p = parts[i]; + if( is_in_domain(p,x) ){ + if( idx != NULL ) *idx = i; + return p; + } + } + char msg[2048]; + sprintf(msg,"Could not find thyla_part for x = (%f,%f,%f)", x[0],x[1],x[2]); + error->one(FLERR,msg); + return NULL; +} + + + + + + +void manifold_thylakoid::init_domains() +{ + if( wB + 2*lB > LT ){ + char msg[2048]; + sprintf(msg,"LT = %f not large enough to accomodate bridge with " + "wB = %f and lB = %f! %f > %f\n", LT, wB, lB, wB + 2*lB, LT); + error->one(FLERR,msg); + } + + // Determine some constant coordinates: + x0 = -( 0.5*LB + lB + lT + LT + lT + pad); + y0 = -( 0.5*LT + lT + pad ); + z0 = -15; +#ifdef MANIFOLD_THYLAKOID_DEBUG + if( comm->me == 0 ){ + fprintf(screen,"x0, y0, z0 = %f, %f, %f\n",x0,y0,z0); + } +#endif // MANIFOLD_THYLAKOID_DEBUG + +#ifndef USE_PHONY_LAMMPS + if( x0 < domain->boxlo[0] ){ + char msg[2048]; + sprintf(msg,"Thylakoid expects xlo of at most %f, but found %f", + x0, domain->boxlo[0]); + error->one(FLERR,msg); + } + if( y0 < domain->boxlo[1] ){ + char msg[2048]; + sprintf(msg,"Thylakoid expects ylo of at most %f, but found %f", + y0, domain->boxlo[1]); + error->one(FLERR,msg); + } +#endif + + // Add some padding to prevent improper lookups. + z0 -= pad; + + x1 = -x0; + y1 = -y0; + z1 = -z0; + + Lx = x1 - x0; + Ly = y1 - y0; + Lz = z1 - z0; + +#ifndef USE_PHONY_LAMMPS + char msg[2048]; + if(x1 > domain->boxhi[0]){ + sprintf(msg,"Expected xhi larger than current box has: %f > %f", + x1, domain->boxhi[0]); + error->one(FLERR,msg); + } + if(y1 > domain->boxhi[1]){ + sprintf(msg,"Expected yhi larger than current box has: %f > %f", + y1, domain->boxhi[1]); + error->one(FLERR,msg); + } + // if(z1 > domain->boxhi[2]){ + // sprintf(msg,"Expected zhi larger than current box has: %f > %f", + // z1, domain->boxhi[2]); + // error->one(FLERR,msg); + // } +#endif + + // Create and add the manifold parts to the array. + thyla_part *p; + + + // Determine coordinates of domain boundaries and centres of "mass": + thyla_part_geom cllb, cllt, clrb, clrt; // Left thylakoid cylinder parts + thyla_part_geom pll, plb, plt, plr; // Left thylakoid plane parts + + thyla_part_geom crlb, crlt, crrb, crrt; // Right thylakoid cylinder parts + thyla_part_geom prl, prb, prt, prr; // Right thylakoid plane parts + + thyla_part_geom bl, br, bc; // Bridge left, right connectors and cylinder. + + + + + + // The bridge is three parts. + // 1. A connector between the right face of the left grana and a cylinder + // 2. A connector between the left face of the right grana and a cylinder + // 3. The aforementioned cylinder. + + // 1.: + // Args: pt, X0, R0, R, x0, y0, z0, sign + double rB = 0.5*wB; + double Reff = rB + lB; + double safety_fac = 1.0; + bl.pt[0] = -0.5*LB; + bl.pt[1] = 0; + bl.pt[2] = 0; + + bl.lo[0] = bl.pt[0] - safety_fac*lB; + bl.lo[1] = -(1.0 + safety_fac) * Reff; + bl.lo[2] = -(1.0 + safety_fac) * Reff; + + bl.hi[0] = bl.pt[0]; + bl.hi[1] = (1.0 + safety_fac) * Reff; + bl.hi[2] = (1.0 + safety_fac) * Reff; + + // double X0, double R0, double R, double s, +#ifdef MANIFOLD_THYLAKOID_DEBUG + if( comm->me == 0 ){ + fprintf(screen,"x0, r0, R = %f, %f, %f\n", bl.pt[0], rB, lB); + } +#endif // MANIFOLD_THYLAKOID_DEBUG + p = make_cyl_to_plane_part(bl.pt[0], rB, lB, -1, bl.pt); + set_domain(p, bl.lo, bl.hi); + parts.push_back(p); + + // 2.: + br.pt[0] = 0.5*LB; + br.pt[1] = 0; + br.pt[2] = 0; + + br.lo[0] = br.pt[0]; + br.lo[1] = -(1.0 + safety_fac) * Reff; + br.lo[2] = -(1.0 + safety_fac) * Reff; + + br.hi[0] = br.pt[0] + safety_fac*lB; + br.hi[1] = (1.0 + safety_fac) * Reff; + br.hi[2] = (1.0 + safety_fac) * Reff; + + // double X0, double R0, double R, double s, +#ifdef MANIFOLD_THYLAKOID_DEBUG + if( comm->me == 0 ){ + fprintf(screen,"x0, r0, R = %f, %f, %f\n", br.pt[0], rB, lB); + } +#endif // MANIFOLD_THYLAKOID_DEBUG + p = make_cyl_to_plane_part(br.pt[0], rB, lB, 1, br.pt); + set_domain(p, br.lo, br.hi); + parts.push_back(p); + + + + // 3.: + // Cylinder in between: + bc.pt[0] = 0; + bc.pt[1] = 0; + bc.pt[2] = 0; + + bc.lo[0] = bl.pt[0]; + bc.lo[1] = -Reff; + bc.lo[2] = -Reff; + + bc.hi[0] = br.pt[0]; + bc.hi[1] = Reff; + bc.hi[2] = Reff; + + p = make_cyl_part( 0, 1, 1, bc.pt, rB ); + set_domain( p, bc.lo, bc.hi ); +#ifdef MANIFOLD_THYLAKOID_DEBUG + if( comm->me == 0 ){ + fprintf(screen,"Cylinder lives on [ %f x %f ] x [ %f x %f ] x [ %f x %f]\n", + bc.lo[0], bc.hi[0], bc.lo[1], bc.hi[1], bc.lo[2], bc.hi[2]); + } +#endif // MANIFOLD_THYLAKOID_DEBUG + + parts.push_back(p); + + + // The stack on the left: + cllb.lo[0] = x0; + cllb.lo[1] = y0; + cllb.lo[2] = z0; + cllb.pt[0] = x0 + pad + lT; + cllb.pt[1] = y0 + pad + lT; + cllb.pt[2] = 0; + cllb.hi[0] = cllb.pt[0]; + cllb.hi[1] = cllb.pt[1]; + cllb.hi[2] = z1; + + p = make_cyl_part(1,1,0,cllb.pt,lT); + set_domain(p, cllb.lo, cllb.hi); + parts.push_back(p); + + // left left top cylinder + cllt = cllb; + cllt.lo[1] = y1 - pad - lT; + cllt.hi[1] = y1; + cllt.pt[1] = cllb.pt[1] + LT; + + p = make_cyl_part(1,1,0,cllt.pt,lT); + set_domain(p, cllt.lo, cllt.hi); + parts.push_back(p); + + // left right bottom cylinder + clrb = cllb; + clrb.pt[0] += LT; + clrb.lo[0] = clrb.pt[0]; + clrb.hi[0] = clrb.lo[0] + lT + lB; + + p = make_cyl_part(1,1,0,clrb.pt,lT); + set_domain(p, clrb.lo, clrb.hi); + parts.push_back(p); + + // left right top cylinder + clrt = clrb; + clrt.pt[1] += LT; + clrt.lo[1] = y1 - pad - lT; + clrt.hi[1] = y1; + + p = make_cyl_part(1,1,0,clrt.pt,lT); + set_domain(p, clrt.lo, clrt.hi); + parts.push_back(p); + + + // left left plane + pll.pt[0] = x0 + pad; + pll.pt[1] = 0; + pll.pt[2] = 0; + pll.lo[0] = x0; + pll.lo[1] = cllb.pt[1]; + pll.lo[2] = z0; + pll.hi[0] = pll.lo[0] + pad + lT; + pll.hi[1] = pll.lo[1] + LT; + pll.hi[2] = z1; + + p = make_plane_part(1,0,0,pll.pt); + set_domain(p, pll.lo, pll.hi); + parts.push_back(p); + + // left bottom plane + plb.pt[0] = x0 + pad + lT + 0.5*LT; + plb.pt[1] = y0 + pad; + plb.pt[2] = 0; + plb.lo[0] = x0 + pad + lT; + plb.lo[1] = y0; + plb.lo[2] = z0; + plb.hi[0] = plb.lo[0] + LT; + plb.hi[1] = plb.lo[1] + pad + lT; + plb.hi[2] = z1; + + p = make_plane_part(0,1,0,plb.pt); + set_domain(p, plb.lo, plb.hi); + parts.push_back(p); + + // left top plane + plt = plb; + plt.lo[1] = cllb.pt[1] + LT; + plt.hi[1] = y1; + plt.pt[1] = y1 - pad; + + p = make_plane_part(0,1,0,plt.pt); + set_domain(p, plt.lo, plt.hi); + parts.push_back(p); + + // left right plane + plr = pll; + plr.lo[0] = bl.lo[0] - 0.5; + plr.lo[1] = y0 - pad; + plr.hi[0] = bl.lo[0] + 0.5; + plr.hi[1] = y1 + pad; + plr.pt[0] = bl.pt[0] - lB; + plr.pt[1] = 0.0; + plr.pt[2] = 0.0; + plr.hi[2] = z1 + pad; + plr.lo[2] = z0 - pad; + + p = make_plane_part(1,0,0,plr.pt); + set_domain(p, plr.lo, plr.hi); + parts.push_back(p); + + // Check if this plane lines up with bl: + if( fabs(plr.pt[0] - bl.pt[0] + lB) > 1e-8 ){ + char msg[2048]; + sprintf(msg,"Origins of plane left right and bridge left misaligned! %f != %f!\n", + plr.pt[0], bl.pt[0] - lB ); + error->one(FLERR,msg); + } + + // Now, for the right stack, you can mirror the other... + // To mirror them you need to invert lo[0] and hi[0] and flip their sign. + + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crlb, &cllb); + p = make_cyl_part(1,1,0,crlb.pt,lT); + set_domain(p, crlb.lo, crlb.hi); + parts.push_back(p); + + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crlt, &cllt); + p = make_cyl_part(1,1,0,crlt.pt,lT); + set_domain(p, crlt.lo, crlt.hi); + parts.push_back(p); + + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crrb, &clrb); + p = make_cyl_part(1,1,0,crrb.pt,lT); + set_domain(p, crrb.lo, crrb.hi); + parts.push_back(p); + + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crrt, &clrt); + p = make_cyl_part(1,1,0,crrt.pt,lT); + set_domain(p, crrt.lo, crrt.hi); + parts.push_back(p); + + + + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prl , &pll ); + p = make_plane_part(1,0,0,prl.pt); + set_domain(p, prl.lo, prl.hi); + parts.push_back(p); + + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prb , &plb ); + p = make_plane_part(0,1,0,prb.pt); + set_domain(p, prb.lo, prb.hi); + parts.push_back(p); + + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prt , &plt ); + p = make_plane_part(0,1,0,prt.pt); + set_domain(p, prt.lo, prt.hi); + parts.push_back(p); + + // Careful, this one is wrongly named. + thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prr, &plr); + p = make_plane_part(1,0,0,prr.pt); + set_domain(p, prr.lo, prr.hi); + parts.push_back(p); + + if( fabs(prr.pt[0] - br.pt[0] - lB) > 1e-8 ){ + char msg[2048]; + sprintf(msg,"Origins of plane left right and bridge left misaligned! %f != %f!\n", + prr.pt[0], br.pt[0] + lB); + error->one(FLERR,msg); + } + + // For debugging, print the domains and coms: +#ifdef MANIFOLD_THYLAKOID_DEBUG + if( comm->me == 0 ){ + FILE *fp_doms = fopen("test_doms.dat","w"); + FILE *fp_coms = fopen("test_coms.dat","w"); + print_part_data(fp_doms, fp_coms); + fclose(fp_doms); + fclose(fp_coms); + } +#endif // MANIFOLD_THYLAKOID_DEBUG +} + + +void manifold_thylakoid::set_domain( thyla_part *p, const std::vector &lo, + const std::vector &hi ) +{ +#ifdef MANIFOLD_THYLAKOID_DEBUG + if( comm->me == 0 ){ + fprintf(screen,"Adding part with domain [%f, %f] x [%f, %f] x [%f, %f]\n", + lo[0],hi[0],lo[1],hi[1],lo[2],hi[2] ); + } +#endif // MANIFOLD_THYLAKOID_DEBUG + if( lo[0] >= hi[0] ){ + char msg[2048]; + sprintf(msg,"xlo >= xhi (%f >= %f)",lo[0],hi[0]); + error->one(FLERR,msg); + }else if( lo[1] >= hi[1] ){ + char msg[2048]; + sprintf(msg,"ylo >= yhi (%f >= %f)",lo[1],hi[1]); + error->one(FLERR,msg); + }else if( lo[2] >= hi[2] ){ + char msg[2048]; + sprintf(msg,"zlo >= zhi (%f >= %f)",lo[2],hi[2]); + error->one(FLERR,msg); + } + p->xlo = lo[0]; + p->ylo = lo[1]; + p->zlo = lo[2]; + + p->xhi = hi[0]; + p->yhi = hi[1]; + p->zhi = hi[2]; +} + +int manifold_thylakoid::is_in_domain( thyla_part *part, const double *x ) +{ + bool domain_ok = (x[0] >= part->xlo) && (x[0] <= part->xhi) && + (x[1] >= part->ylo) && (x[1] <= part->yhi) && + (x[2] >= part->zlo) && (x[2] <= part->zhi); + + if( !domain_ok ) return false; + + // From here on out, domain is ok. + + if( part->type == thyla_part::THYLA_TYPE_CYL_TO_PLANE ){ + + double R0 = part->params[1]; + double R = part->params[2]; + double y = x[1]; + double z = x[2]; + double dist2 = y*y + z*z; + double RR = R0+R; + double RR2 = RR*RR; + + + if( dist2 < RR2 ){ + return true; + }else{ + // Domain was ok, but radius not. + return false; + } + }else{ + return true; + } +} + + +thyla_part *manifold_thylakoid::make_plane_part (double a, double b, double c, + const std::vector &pt ) +{ + double args[7]; + args[0] = a; + args[1] = b; + args[2] = c; + args[3] = pt[0]; + args[4] = pt[1]; + args[5] = pt[2]; + thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_PLANE,args,0,0,0,0,0,0); + return p; +} + +thyla_part *manifold_thylakoid::make_cyl_part (double a, double b, double c, + const std::vector &pt, double R) +{ + double args[7]; + args[0] = a; + args[1] = b; + args[2] = c; + args[3] = pt[0]; + args[4] = pt[1]; + args[5] = pt[2]; + args[6] = R; + thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_CYL,args,0,0,0,0,0,0); + return p; +} + + +thyla_part *manifold_thylakoid::make_sphere_part(const std::vector &pt, double R) +{ + double args[7]; + args[0] = R; + args[1] = pt[0]; + args[2] = pt[1]; + args[3] = pt[2]; + thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_SPHERE,args,0,0,0,0,0,0); + return p; +} + + +thyla_part *manifold_thylakoid::make_cyl_to_plane_part(double X0, double R0, double R, + double s, const std::vector &pt ) +{ + double args[7]; + args[0] = X0; + args[1] = R0; + args[2] = R; + args[3] = pt[0]; + args[4] = pt[1]; + args[5] = pt[2]; + args[6] = s; + thyla_part *p = new thyla_part(thyla_part::THYLA_TYPE_CYL_TO_PLANE,args,0,0,0,0,0,0); + return p; +} + + + + +void manifold_thylakoid::print_part_data( FILE *fp_doms, FILE *fp_coms ) +{ + for( std::size_t i = 0; i < parts.size(); ++i ){ + thyla_part *p = parts[i]; + fprintf(fp_doms, "%f %f\n", p->xlo, p->ylo); + fprintf(fp_doms, "%f %f\n", p->xlo, p->yhi); + fprintf(fp_doms, "%f %f\n", p->xhi, p->yhi); + fprintf(fp_doms, "%f %f\n", p->xhi, p->ylo); + fprintf(fp_doms, "%f %f\n\n",p->xlo, p->ylo); + fprintf(fp_coms, "%f %f\n", p->x0, p->y0 ); + } +} + + diff --git a/src/USER-MANIFOLD/manifold_thylakoid.h b/src/USER-MANIFOLD/manifold_thylakoid.h new file mode 100644 index 0000000000..231525950d --- /dev/null +++ b/src/USER-MANIFOLD/manifold_thylakoid.h @@ -0,0 +1,71 @@ +#ifndef LMP_MANIFOLD_THYLAKOID_H +#define LMP_MANIFOLD_THYLAKOID_H + +#include "manifold.h" +#include +#include + +#include "manifold_thylakoid_shared.h" + +namespace LAMMPS_NS { + +namespace user_manifold { + + + class manifold_thylakoid : public manifold { + public: + enum { NPARAMS = 3 }; + manifold_thylakoid( LAMMPS *lmp, int, char ** ); + virtual ~manifold_thylakoid(); + + virtual double g( const double *x ); + virtual void n( const double *x, double *n ); + + static const char* type(){ return "thylakoid"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + + + virtual void post_param_init(); + virtual void checkup(); // Some diagnostics... + private: + void init_domains(); + + thyla_part *get_thyla_part( const double *x, int *err_flag, std::size_t *idx = NULL ); + int is_in_domain( thyla_part *p, const double *x ); + void check_overlap(); + std::vector parts; + + thyla_part *make_plane_part (double a, double b, double c, + const std::vector &pt); + thyla_part *make_cyl_part (double a, double b, double c, + const std::vector &pt, double R); + thyla_part *make_sphere_part(const std::vector &pt, double R); + thyla_part *make_cyl_to_plane_part(double X0, double R0, double R, double s, + const std::vector &pt ); + + + void set_domain( thyla_part *p, const std::vector &lo, + const std::vector &hi ); + + void print_part_data( FILE *fp_doms, FILE *fp_coms ); + + // Coefficients for the thylakoid model. At the moment it is just + // a cylinder, we slowly expand it. + double pad; // Padding to make sure periodic images are mapped back properly. + double LB, lT, lB, wB, LT; + + // Domain size: + double x0, y0, z0; + double x1, y1, z1; + double Lx, Ly, Lz; + }; + + + +} // namespace LAMMPS_NS + +} + +#endif // LMP_MANIFOLD_THYLAKOID_H diff --git a/src/USER-MANIFOLD/manifold_thylakoid_shared.cpp b/src/USER-MANIFOLD/manifold_thylakoid_shared.cpp new file mode 100644 index 0000000000..5b87ba7c2b --- /dev/null +++ b/src/USER-MANIFOLD/manifold_thylakoid_shared.cpp @@ -0,0 +1,228 @@ +#include "manifold_thylakoid_shared.h" +#include + +using namespace LAMMPS_NS; +using namespace user_manifold; + + +thyla_part::thyla_part( int type, double *args, double xlo, double ylo, double zlo, + double xhi, double yhi, double zhi ) + : type(type), xlo(xlo), xhi(xhi), + ylo(ylo), yhi(yhi), zlo(zlo), zhi(zhi) +{ + switch(type){ + case THYLA_TYPE_PLANE: // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0 + params[0] = args[0]; // a + params[1] = args[1]; // b + params[2] = args[2]; // c + params[3] = args[3]; // x0 + params[4] = args[4]; // y0 + params[5] = args[5]; // z0 + break; + case THYLA_TYPE_SPHERE: // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0 + params[0] = args[0]; // R + params[1] = args[1]; // x0 + params[2] = args[2]; // y0 + params[3] = args[3]; // z0 + break; + case THYLA_TYPE_CYL: // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0 + params[0] = args[0]; // a + params[1] = args[1]; // b + params[2] = args[2]; // c + params[3] = args[3]; // x0 + params[4] = args[4]; // y0 + params[5] = args[5]; // z0 + params[6] = args[6]; // R + if( (args[0] != 0.0) && (args[1] != 0.0) && (args[2] != 0.0) ){ + err_flag = -1; + return; + } + // The others should be 1. + if( (args[0] != 1.0) && (args[0] != 0.0) && + (args[1] != 1.0) && (args[1] != 0.0) && + (args[2] != 1.0) && (args[2] != 0.0) ){ + err_flag = -1; + } + break; + case THYLA_TYPE_CYL_TO_PLANE: + /* + * Funky bit that connects a cylinder to a plane. + * It is what you get by rotating the equation + * r(x) = R0 + R*( 1 - sqrt( 1 - ( ( X0 - x ) /R )^2 ) ) around the x-axis. + * I kid you not. + * + * The shape is symmetric so you have to set whether it is the "left" or + * "right" end by truncating it properly. It is designed to run from + * X0 to X0 + R if "right" or X0 - R to X0 if "left". + * + * As params it takes X0, R0, R, and a sign that determines whether it is + * "left" (args[3] < 0) or "right" (args[3] > 0). + * + * The args are: X0, R0, R, x0, y0, z0, sign. + */ + params[0] = args[0]; + params[1] = args[1]; + params[2] = args[2]; + params[3] = args[3]; + params[4] = args[4]; + params[5] = args[5]; + params[6] = args[6]; + break; + default: + err_flag = -1; + } + x0 = (type == THYLA_TYPE_SPHERE) ? params[1] : params[3]; + y0 = (type == THYLA_TYPE_SPHERE) ? params[2] : params[4]; + z0 = (type == THYLA_TYPE_SPHERE) ? params[3] : params[5]; +} + + + +thyla_part::~thyla_part() +{} + +double thyla_part::g(const double *x) +{ + switch(type){ + case THYLA_TYPE_PLANE:{ // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0 + double a = params[0]; + double b = params[1]; + double c = params[2]; + double dx = x[0] - params[3]; + double dy = x[1] - params[4]; + double dz = x[2] - params[5]; + return a*dx + b*dy + c*dz; + break; + } + case THYLA_TYPE_SPHERE:{ // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0 + double R2 = params[0]*params[0]; + double dx = x[0] - params[1]; + double dy = x[1] - params[2]; + double dz = x[2] - params[3]; + return dx*dx + dy*dy + dz*dz - R2; + + break; + } + case THYLA_TYPE_CYL:{ // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0 + double a = params[0]; + double b = params[1]; + double c = params[2]; + double X0 = params[3]; + double Y0 = params[4]; + double Z0 = params[5]; + double R = params[6]; + double dx = x[0] - X0; + double dy = x[1] - Y0; + double dz = x[2] - Z0; + + return a*dx*dx + b*dy*dy + c*dz*dz - R*R; + break; + } + case THYLA_TYPE_CYL_TO_PLANE:{ + double X0 = params[0]; + double R0 = params[1]; + double R = params[2]; + + // Determine the centre of the sphere. + double dx = (x[0] - X0); + double dyz = sqrt( x[1]*x[1] + x[2]*x[2] ); + double rdyz = dyz - (R0 + R); + double norm = 1.0 / (2.0*R); + // Maybe sign is important here... + double g = norm*(dx*dx + rdyz*rdyz - R*R); + return g; + + } + default: + err_flag = -1; + return 0; // Mostly to get rid of compiler werrors. + break; + } +} + + + +void thyla_part::n( const double *x, double *n ) +{ + switch(type){ + case THYLA_TYPE_PLANE:{ // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0 + double a = params[0]; + double b = params[1]; + double c = params[2]; + n[0] = a; + n[1] = b; + n[2] = c; + break; + } + case THYLA_TYPE_SPHERE:{ // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0 + double dx = x[0] - params[1]; + double dy = x[1] - params[2]; + double dz = x[2] - params[3]; + n[0] = 2*dx; + n[1] = 2*dy; + n[2] = 2*dz; + break; + } + case THYLA_TYPE_CYL:{ // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0 + double a = params[0]; + double b = params[1]; + double c = params[2]; + double X0 = params[3]; + double Y0 = params[4]; + double Z0 = params[5]; + double dx = x[0] - X0; + double dy = x[1] - Y0; + double dz = x[2] - Z0; + + n[0] = 2*a*dx; + n[1] = 2*b*dy; + n[2] = 2*c*dz; + break; + } + case THYLA_TYPE_CYL_TO_PLANE:{ + double X0 = params[0]; + double R0 = params[1]; + double R = params[2]; + double s = (params[6] > 0.0) ? 1.0 : -1.0; + + // Determine the centre of the sphere. + double dx = s*(x[0] - X0); + double ryz = sqrt( x[1]*x[1] + x[2]*x[2] ); + // Maybe sign is important here... + // Normalize g and n so that the normal is continuous: + double norm = 1.0 / (2.0 * R); + + n[0] = s*2*dx*norm; + + double const_part = 1.0 - (R0 + R) / ryz; + n[1] = 2*x[1]*const_part*norm; + n[2] = 2*x[2]*const_part*norm; + break; + } + default: + err_flag = -1; + break; + } +} + + +void thyla_part_geom::mirror( unsigned int axis, thyla_part_geom *m, + const thyla_part_geom *o ) +{ + // Since dir is already the index of the array this is really simple: + m->lo[0] = o->lo[0]; + m->lo[1] = o->lo[1]; + m->lo[2] = o->lo[2]; + m->pt[0] = o->pt[0]; + m->pt[1] = o->pt[1]; + m->pt[2] = o->pt[2]; + m->hi[0] = o->hi[0]; + m->hi[1] = o->hi[1]; + m->hi[2] = o->hi[2]; + + m->lo[axis] = -o->hi[axis]; + m->hi[axis] = -o->lo[axis]; + m->pt[axis] = -o->pt[axis]; +} + + diff --git a/src/USER-MANIFOLD/manifold_thylakoid_shared.h b/src/USER-MANIFOLD/manifold_thylakoid_shared.h new file mode 100644 index 0000000000..756f392bbb --- /dev/null +++ b/src/USER-MANIFOLD/manifold_thylakoid_shared.h @@ -0,0 +1,63 @@ +#ifndef MANIFOLD_THYLAKOID_SHARED_H +#define MANIFOLD_THYLAKOID_SHARED_H + + +#include "lmptype.h" +#include + +namespace LAMMPS_NS { + +namespace user_manifold { + + + // The thylakoid is composed of many parts + struct thyla_part { + enum thyla_types { + THYLA_TYPE_PLANE, + THYLA_TYPE_SPHERE, + THYLA_TYPE_CYL, + THYLA_TYPE_CYL_TO_PLANE + }; + + thyla_part( int type, double *args, double xlo, double ylo, double zlo, + double xhi, double yhi, double zhi ); + thyla_part() : type(-1), x0(-1337), y0(-1337), z0(-1337){} + ~thyla_part(); + + double g( const double *x ); + void n( const double *x, double *n ); + + int type; + double params[7]; + double tol; + int maxit; + + int err_flag; + tagint a_id; + + double xlo, xhi, ylo, yhi, zlo, zhi; + double x0, y0, z0; + + }; // struct thyla_part + + + + struct thyla_part_geom { + thyla_part_geom() : pt(3), lo(3), hi(3){} + std::vector pt, lo, hi; + + // Function for mirroring thyla_geoms: + enum DIRS { DIR_X, DIR_Y, DIR_Z }; + static void mirror( unsigned int axis, thyla_part_geom *m, + const thyla_part_geom *o ); + + }; // struct thyla_part_geom + + +} // namespace user_manifold + + +} // namespace LAMMPS_NS + + +#endif // MANIFOLD_THYLAKOID_SHARED_H diff --git a/src/USER-MANIFOLD/manifold_torus.cpp b/src/USER-MANIFOLD/manifold_torus.cpp new file mode 100644 index 0000000000..b46998b28d --- /dev/null +++ b/src/USER-MANIFOLD/manifold_torus.cpp @@ -0,0 +1,40 @@ +#include +#include "manifold_torus.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace user_manifold; + + +manifold_torus::manifold_torus( LAMMPS *lmp, int argc, char **argv ) : manifold(lmp) +{} + + +double manifold_torus::g( const double *x ) +{ + double R = params[0]; + double r = params[1]; + if( R < r ){ + error->all(FLERR,"Large radius < small radius!"); + } + + double rad = sqrt(x[0]*x[0] + x[1]*x[1]); + double c = R - rad; + return c*c + x[2]*x[2] - r*r; +} + +void manifold_torus::n( const double *x, double *n ) +{ + double R = params[0]; + double r = params[1]; + if( R < r ){ + error->all(FLERR,"Large radius < small radius!"); + } + + double rad = sqrt(x[0]*x[0] + x[1]*x[1]); + double c = R - rad; + double fac = c / rad; + n[0] = -2.0 * fac * x[0]; + n[1] = -2.0 * fac * x[1]; + n[2] = 2.0 * x[2]; +} diff --git a/src/USER-MANIFOLD/manifold_torus.h b/src/USER-MANIFOLD/manifold_torus.h new file mode 100644 index 0000000000..0c3bd6e2eb --- /dev/null +++ b/src/USER-MANIFOLD/manifold_torus.h @@ -0,0 +1,32 @@ +#ifndef LMP_MANIFOLD_TORUS_H +#define LMP_MANIFOLD_TORUS_H + +#include "manifold.h" + + +namespace LAMMPS_NS { + +namespace user_manifold { + + + class manifold_torus : public manifold { + public: + enum {NPARAMS=2}; + manifold_torus( LAMMPS *, int, char ** ); + ~manifold_torus(){} + virtual double g( const double *x ); + virtual void n( const double *x, double *n ); + + static const char *type(){ return "torus"; } + virtual const char *id(){ return type(); } + static int expected_argc(){ return NPARAMS; } + virtual int nparams(){ return NPARAMS; } + }; + +} + +} + + + +#endif // LMP_MANIFOLD_TORUS_H