mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14994 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
244889aed4
commit
2b75c78f69
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 <your_architecture>
|
||||
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
|
||||
+==============================================================================+
|
||||
|
|
@ -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 <name> <group> 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);
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
@ -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<double>(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;
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
|
||||
|
||||
*/
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
|
||||
|
||||
*/
|
||||
|
|
@ -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 <math.h>
|
||||
|
||||
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
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
#include "manifold_cylinder_dent.h"
|
||||
#include "math_const.h"
|
||||
|
||||
#include <math.h>
|
||||
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,43 @@
|
|||
#include "manifold_dumbbell.h"
|
||||
|
||||
#include <math.h>
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -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
|
||||
|
|
@ -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;
|
||||
}
|
||||
|
||||
|
|
@ -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
|
||||
|
|
@ -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<manifold_cylinder> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_cylinder_dent> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_dumbbell> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_ellipsoid> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_plane> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_plane_wiggle> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_sphere> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_supersphere> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_spine> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_spine_two> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_thylakoid> ( &man, mname, lmp, narg, arg );
|
||||
make_manifold_if<manifold_torus> ( &man, mname, lmp, narg, arg );
|
||||
|
||||
|
||||
return man;
|
||||
}
|
||||
|
||||
|
|
@ -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 <string.h>
|
||||
|
||||
/*
|
||||
* 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 <typename m_type>
|
||||
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
|
||||
|
|
@ -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];
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,28 @@
|
|||
#include "manifold_plane_wiggle.h"
|
||||
|
||||
#include <math.h>
|
||||
|
||||
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]);
|
||||
}
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,163 @@
|
|||
#include "manifold_spine.h"
|
||||
|
||||
#include <math.h>
|
||||
#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 );
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
@ -0,0 +1,618 @@
|
|||
#include "manifold_thylakoid.h"
|
||||
#include <math.h>
|
||||
|
||||
#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<double> &lo,
|
||||
const std::vector<double> &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<double> &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<double> &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<double> &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<double> &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 );
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,71 @@
|
|||
#ifndef LMP_MANIFOLD_THYLAKOID_H
|
||||
#define LMP_MANIFOLD_THYLAKOID_H
|
||||
|
||||
#include "manifold.h"
|
||||
#include <vector>
|
||||
#include <stdio.h>
|
||||
|
||||
#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<thyla_part*> parts;
|
||||
|
||||
thyla_part *make_plane_part (double a, double b, double c,
|
||||
const std::vector<double> &pt);
|
||||
thyla_part *make_cyl_part (double a, double b, double c,
|
||||
const std::vector<double> &pt, double R);
|
||||
thyla_part *make_sphere_part(const std::vector<double> &pt, double R);
|
||||
thyla_part *make_cyl_to_plane_part(double X0, double R0, double R, double s,
|
||||
const std::vector<double> &pt );
|
||||
|
||||
|
||||
void set_domain( thyla_part *p, const std::vector<double> &lo,
|
||||
const std::vector<double> &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
|
||||
|
|
@ -0,0 +1,228 @@
|
|||
#include "manifold_thylakoid_shared.h"
|
||||
#include <math.h>
|
||||
|
||||
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];
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,63 @@
|
|||
#ifndef MANIFOLD_THYLAKOID_SHARED_H
|
||||
#define MANIFOLD_THYLAKOID_SHARED_H
|
||||
|
||||
|
||||
#include "lmptype.h"
|
||||
#include <vector>
|
||||
|
||||
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<double> 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
|
||||
|
|
@ -0,0 +1,40 @@
|
|||
#include <math.h>
|
||||
#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];
|
||||
}
|
||||
|
|
@ -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
|
||||
Loading…
Reference in New Issue