// This file is part of the ESPResSo distribution (http://www.espresso.mpg.de).
// It is therefore subject to the ESPResSo license agreement which you accepted upon receiving the distribution
// and by which you are legally bound while utilizing this file in any form or way.
// There is NO WARRANTY, not even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// You should have received a copy of that license along with this program;
// if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
// write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
// Copyright (c) 2002-2006; all rights reserved unless otherwise stated.
#ifndef IA_DATA_H
#define IA_DATA_H
/** \file interaction_data.h
Various procedures concerning interactions between particles.
Responsible:
Axel
interaction_data.h contains the parser \ref #inter for the
\ref tcl_inter Tcl command. Therefore the parsing of bonded and nonbonded
interaction definition commands both is done here. It also contains
procedures for low-level interactions setup and some helper functions.
Moreover it contains code for the treatments of constraints.
*/
#include
#include "utils.h"
#include "particle_data.h" /* needed for constraints */
/** \name Type codes of bonded interactions
Enumeration of implemented bonded interactions.
*/
/************************************************************/
/address@hidden/
/** This bonded interaction was not set. */
#define BONDED_IA_NONE -1
/** Type of bonded interaction is a FENE potential
(to be combined with Lennard Jones). */
#define BONDED_IA_FENE 0
/** Type of bonded interaction is a HARMONIC potential. */
#define BONDED_IA_HARMONIC 1
/** Type of bonded interaction is a bond angle potential. */
#define BONDED_IA_ANGLE 2
/** Type of bonded interaction is a dihedral potential. */
#define BONDED_IA_DIHEDRAL 3
/** Type of tabulated bonded interaction potential,
may be of bond length, of bond angle or of dihedral type. */
#define BONDED_IA_TABULATED 4
/** Type of bonded interaction is a (-LJ) potential. */
#define BONDED_IA_SUBT_LJ 5
/** Type of a Rigid/Constrained bond*/
#define BONDED_IA_RIGID_BOND 6
/** Type of a virtual bond*/
#define BONDED_IA_VIRTUAL_BOND 7
/** Type of bonded interaction is a dihedral potential.
Identical to BONDED_IA_DIHEDRAL except quaternions of two particles give the dihedral angle
In this case we need five sets of interetions, three for crankshift and two for the orthogonal
direction. */
#define BONDED_IA_GB_DIHEDRAL 8
/* Identical to BONDED_IA_SUBT_LJ except quaternions of two particles give the dihedral angle */
#define BONDED_IA_SUBT_GB 9
/* Specify tabulated bonded interactions */
#define TAB_UNKNOWN 0
#define TAB_BOND_LENGTH 1
#define TAB_BOND_ANGLE 2
#define TAB_BOND_DIHEDRAL 3
/address@hidden/
/** \name Type codes for the type of Coulomb interaction
Enumeration of implemented methods for the electrostatic
interaction.
*/
/************************************************************/
/address@hidden/
#ifdef ELECTROSTATICS
/** Coulomb interation switched off (NONE). */
#define COULOMB_NONE 0
/** Coulomb method is Debye-Hueckel. */
#define COULOMB_DH 1
/** Coulomb method is Debye-Hueckel with parallel separate calculation. */
#define COULOMB_DH_PW 2
/** Coulomb method is P3M. */
#define COULOMB_P3M 3
/** Coulomb method is one-dimensional MMM */
#define COULOMB_MMM1D 4
/** Coulomb method is two-dimensional MMM */
#define COULOMB_MMM2D 5
/** Coulomb method is "Maggs" */
#define COULOMB_MAGGS 6
/** Coulomb method is standard Ewald */
#define COULOMB_EWALD 7
/** Coulomb method is P3M plus ELC. */
#define COULOMB_ELC_P3M 8
/** Coulomb method is Reaction-Field. */
#define COULOMB_RF 9
/** Coulomb method is Reaction-Field BUT as interactions */
#define COULOMB_INTER_RF 10
#endif
/address@hidden/
#ifdef MAGNETOSTATICS
/** \name Type codes for the type of dipolar interaction
Enumeration of implemented methods for the magnetostatic
interaction.
*/
/************************************************************/
/address@hidden/
/** dipolar interation switched off (NONE). */
#define DIPOLAR_NONE 0
/** dipolar method is P3M. */
#define DIPOLAR_P3M 1
/** Dipolar method is P3M plus DLC. */
#define DIPOLAR_MDLC_P3M 2
/** Dipolar method is all with all and no replicas */
#define DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA 3
/** Dipolar method is magnetic dipolar direct sum */
#define DIPOLAR_DS 4
/** Dipolar method is direct sum plus DLC. */
#define DIPOLAR_MDLC_DS 5
/address@hidden/
#endif
/** \name Type codes for constraints
Enumeration of implemented constraint types.
*/
/************************************************************/
/address@hidden/
/** No constraint applied */
#define CONSTRAINT_NONE 0
/** wall constraint applied */
#define CONSTRAINT_WAL 1
/** spherical constraint applied */
#define CONSTRAINT_SPH 2
/** (finite) cylinder shaped constraint applied */
#define CONSTRAINT_CYL 3
/** Rod-like charge. */
#define CONSTRAINT_ROD 4
/** Plate-like charge. */
#define CONSTRAINT_PLATE 5
/** maze-like constraint applied */
#define CONSTRAINT_MAZE 6
/** pore constraint applied */
#define CONSTRAINT_PORE 7
//ER
/** External magnetic field constraint applied */
#define CONSTRAINT_EXT_MAGN_FIELD 8
//end ER
/** Constraint for tunable-lsip boundary conditions */
#define CONSTRAINT_PLANE 9
/address@hidden/
/* Data Types */
/************************************************************/
/** field containing the interaction parameters for
* nonbonded interactions. Access via
* get_ia_param(i, j), i,j < n_particle_types */
typedef struct {
#ifdef LENNARD_JONES
/** \name Lennard-Jones with shift */
/address@hidden/
double LJ_eps;
double LJ_sig;
double LJ_cut;
double LJ_shift;
double LJ_offset;
double LJ_capradius;
double LJ_min;
/address@hidden/
#endif
#ifdef LENNARD_JONES_GENERIC
/** \name Generic Lennard-Jones with shift */
/address@hidden/
double LJGEN_eps;
double LJGEN_sig;
double LJGEN_cut;
double LJGEN_shift;
double LJGEN_offset;
double LJGEN_capradius;
int LJGEN_a1;
int LJGEN_a2;
double LJGEN_b1;
double LJGEN_b2;
/address@hidden/
#endif
#ifdef LJ_ANGLE
/** \name Directional Lennard-Jones */
/address@hidden/
double LJANGLE_eps;
double LJANGLE_sig;
double LJANGLE_cut;
/* Locate bonded partners */
int LJANGLE_bonded1type;
int LJANGLE_bonded1pos;
int LJANGLE_bonded1neg;
int LJANGLE_bonded2pos;
int LJANGLE_bonded2neg;
/* Cap */
double LJANGLE_capradius;
/address@hidden/
#endif
#ifdef SMOOTH_STEP
/** \name smooth step potential */
/address@hidden/
double SmSt_eps;
double SmSt_sig;
double SmSt_cut;
double SmSt_d;
int SmSt_n;
double SmSt_k0;
/address@hidden/
#endif
#ifdef BMHTF_NACL
/** \name BMHTF NaCl potential */
/address@hidden/
double BMHTF_A;
double BMHTF_B;
double BMHTF_C;
double BMHTF_D;
double BMHTF_sig;
double BMHTF_cut;
double BMHTF_computed_shift;
/address@hidden/
#endif
#ifdef MORSE
/** \name Morse potential */
/address@hidden/
double MORSE_eps;
double MORSE_alpha;
double MORSE_rmin;
double MORSE_cut;
double MORSE_rest;
double MORSE_capradius;
/address@hidden/
#endif
#ifdef BUCKINGHAM
/** \name Buckingham potential */
/address@hidden/
double BUCK_A;
double BUCK_B;
double BUCK_C;
double BUCK_D;
double BUCK_cut;
double BUCK_discont;
double BUCK_shift;
double BUCK_capradius;
double BUCK_F1;
double BUCK_F2;
/address@hidden/
#endif
#ifdef SOFT_SPHERE
/** \name soft-sphere potential */
/address@hidden/
double soft_a;
double soft_n;
double soft_cut;
double soft_offset;
/address@hidden/
#endif
#ifdef LJCOS
/** \name Lennard-Jones+Cos potential */
/address@hidden/
double LJCOS_eps;
double LJCOS_sig;
double LJCOS_cut;
double LJCOS_offset;
double LJCOS_alfa;
double LJCOS_beta;
double LJCOS_rmin;
/address@hidden/
#endif
#ifdef LJCOS2
/** \name Lennard-Jones with a different Cos potential */
/address@hidden/
double LJCOS2_eps;
double LJCOS2_sig;
double LJCOS2_cut;
double LJCOS2_offset;
double LJCOS2_w;
double LJCOS2_rchange;
double LJCOS2_capradius;
/address@hidden/
#endif
#ifdef ROTATION
/** \name Gay-Berne potential */
/address@hidden/
double GB_eps;
double GB_sig;
double GB_cut;
double GB_k1;
double GB_k2;
double GB_mu;
double GB_nu;
double GB_chi1;
double GB_chi2;
double GB_alpha1;
double GB_alpha2;
/address@hidden/
#endif
#ifdef TABULATED
/** \name Tabulated potential */
/address@hidden/
int TAB_npoints;
int TAB_startindex;
double TAB_minval;
double TAB_minval2;
double TAB_maxval;
double TAB_maxval2;
double TAB_stepsize;
/** The maximum allowable filename length for a tabulated potential file*/
#define MAXLENGTH_TABFILE_NAME 256
char TAB_filename[MAXLENGTH_TABFILE_NAME];
/address@hidden/
#endif
#ifdef COMFORCE
/** \name center of mass directed force */
/address@hidden/
int COMFORCE_flag;
int COMFORCE_dir;
double COMFORCE_force;
double COMFORCE_fratio;
/address@hidden/
#endif
#ifdef COMFIXED
/** \name center of mass directed force */
/address@hidden/
int COMFIXED_flag;
/address@hidden/
#endif
#ifdef INTER_DPD
/** \name DPD as interaction */
/address@hidden/
double dpd_gamma;
double dpd_r_cut;
int dpd_wf;
double dpd_pref1;
double dpd_pref2;
double dpd_tgamma;
double dpd_tr_cut;
int dpd_twf;
double dpd_pref3;
double dpd_pref4;
/address@hidden/
#endif
#ifdef INTER_RF
int rf_on;
#endif
#ifdef MOL_CUT
int mol_cut_type;
double mol_cut_cutoff;
#endif
#ifdef TUNABLE_SLIP
double TUNABLE_SLIP_temp;
double TUNABLE_SLIP_gamma;
double TUNABLE_SLIP_r_cut;
double TUNABLE_SLIP_time;
double TUNABLE_SLIP_vx;
double TUNABLE_SLIP_vy;
double TUNABLE_SLIP_vz;
#endif
} IA_parameters;
/** \name Compounds for Coulomb interactions */
/address@hidden/
/** field containing the interaction parameters for
* the coulomb interaction. */
typedef struct {
#ifdef ELECTROSTATICS
/** Bjerrum length. */
double bjerrum;
/** bjerrum length times temperature. */
double prefactor;
/** Method to treat coulomb interaction. See \ref COULOMB_NONE "Type codes for Coulomb" */
int method;
#endif
#ifdef MAGNETOSTATICS
double Dbjerrum;
double Dprefactor;
int Dmethod;
#endif
} Coulomb_parameters;
/address@hidden/
/** Defines parameters for a bonded interaction. */
typedef struct {
/** bonded interaction type. See \ref BONDED_IA_FENE "Type code for bonded" */
int type;
/** (Number of particles - 1) interacting for that type */
int num;
/** union to store the different bonded interaction parameters. */
union {
/** Parameters for FENE bond Potential.
k - spring constant.
drmax - maximal bond streching.
r0 - equilibrium bond length.
drmax2 - square of drmax (internal parameter).
*/
struct {
double k;
double drmax;
double r0;
double drmax2;
double drmax2i;
} fene;
/** Parameters for harmonic bond Potential */
struct {
double k;
double r;
double r_cut;
} harmonic;
/** Parameters for three body angular potential (bond-angle potentials).
ATTENTION: Note that there are different implementations of the bond angle
potential which you may chose with a compiler flag in the file \ref config.h !
bend - bending constant.
phi0 - equilibrium angle (default is 180 degrees / Pi) */
struct {
double bend;
double phi0;
#ifdef BOND_ANGLE_COSINE
double cos_phi0;
double sin_phi0;
#endif
#ifdef BOND_ANGLE_COSSQUARE
double cos_phi0;
#endif
} angle;
/** Parameters for four body angular potential (dihedral-angle potentials). */
struct {
double mult;
double bend;
double phase;
} dihedral;
#ifdef GB_DIHED
struct {
double gb_mult[5];
double gb_bend[5];
double gb_phase[5];
} gb_dihedral;
#endif
#ifdef TABULATED
/** Parameters for n-body tabulated potential (n=2,3,4). */
struct {
char *filename;
int type;
int npoints;
double minval;
double maxval;
double invstepsize;
double *f;
double *e;
} tab;
#endif
/** Dummy parameters for -LJ Potential */
struct {
double k;
double r;
double r2;
} subt_lj;
/** Dummy parameters for -GB Potential */
struct {
double k;
double r;
double r2;
} subt_gb;
/**Parameters for the rigid_bond/SHAKE/RATTLE ALGORITHM*/
struct {
/**Length of rigid bond/Constrained Bond*/
//double d;
/**Square of the length of Constrained Bond*/
double d2;
/**Positional Tolerance/Accuracy value for termination of RATTLE/SHAKE iterations during position corrections*/
double p_tol;
/**Velocity Tolerance/Accuracy for termination of RATTLE/SHAKE iterations during velocity corrections */
double v_tol;
} rigid_bond;
} p;
} Bonded_ia_parameters;
#ifdef CONSTRAINTS
/** \name Compounds for constraints */
/address@hidden/
/** Parameters for a WALL constraint (or a plane if you like that more). */
typedef struct {
/** normal vector on the plane. */
double n[3];
/** distance of the wall from the origin. */
double d;
} Constraint_wall;
/** Parameters for a SPHERE constraint. */
typedef struct {
/** sphere center. */
double pos[3];
/** sphere radius. */
double rad;
/** sphere direction. (+1 outside -1 inside interaction direction)*/
double direction;
} Constraint_sphere;
/** Parameters for a CYLINDER constraint. */
typedef struct {
/** center of the cylinder. */
double pos[3];
/** Axis of the cylinder .*/
double axis[3];
/** cylinder radius. */
double rad;
/** cylinder length. (!!!NOTE this is only the half length of the cylinder.)*/
double length;
/** cylinder direction. (+1 outside -1 inside interaction direction)*/
double direction;
} Constraint_cylinder;
/** Parameters for a PORE constraint. */
typedef struct {
/** center of the cylinder. */
double pos[3];
/** Axis of the cylinder .*/
double axis[3];
/** cylinder radius. */
double rad;
/** cylinder length. (!!!NOTE this is only the half length of the cylinder.)*/
double length;
} Constraint_pore;
/** Parameters for a ROD constraint. */
typedef struct {
/** center of the cylinder in the x-y plane. */
double pos[2];
/** line charge density. Only makes sense if the axis along the rod is
periodically replicated and only with MMM1D. */
double lambda;
} Constraint_rod;
/** Parameters for a PLATE constraint. */
typedef struct {
/** height of plane in z-axis. */
double pos;
/** charge density. Only makes sense if the axis along the rod is
periodically replicated and only with MMM2D. */
double sigma;
} Constraint_plate;
/** Parameters for a MAZE constraint. */
typedef struct {
/** number of spheres. */
double nsphere;
/** dimension of the maze. */
double dim;
/** sphere radius. */
double sphrad;
/** cylinder (connecting the spheres) radius*/
double cylrad;
} Constraint_maze;
//ER
/** Parameters for a EXTERNAL MAGNETIC FIELD constraint */
typedef struct{
/** vector (direction and magnitude) of the external magnetic field */
double ext_magn_field[3];
} Constraint_ext_magn_field;
//end ER
/** Parameters for a plane constraint which is needed for tunable-slip boundary conditions. */
typedef struct {
/** Position of the plain. Negative values mean non-existing in that direction. */
double pos[3];
} Constraint_plane;
/** Structure to specify a constraint. */
typedef struct {
/** type of the constraint. */
int type;
union {
Constraint_wall wal;
Constraint_sphere sph;
Constraint_cylinder cyl;
Constraint_rod rod;
Constraint_plate plate;
Constraint_maze maze;
Constraint_pore pore;
//ER
Constraint_ext_magn_field emfield;
//end ER
Constraint_plane plane;
} c;
/** particle representation of this constraint. Actually needed are only the identity,
the type and the force. */
Particle part_rep;
} Constraint;
/address@hidden/
#endif
/************************************************
* exported variables
************************************************/
/** Maximal particle type seen so far. */
extern int n_particle_types;
/* Number of nonbonded (short range) interactions. Not used so far.*/
extern int n_interaction_types;
/** Structure containing the coulomb parameters. */
extern Coulomb_parameters coulomb;
/** number of bonded interactions. Not used so far. */
extern int n_bonded_ia;
/** Field containing the paramters of the bonded ia types */
extern Bonded_ia_parameters *bonded_ia_params;
/** Array containing all tabulated forces*/
extern DoubleList tabulated_forces;
/** Array containing all tabulated energies*/
extern DoubleList tabulated_energies;
/** Maximal interaction cutoff (real space/short range interactions). */
extern double max_cut;
/** Maximal interaction cutoff (real space/short range non-bonded interactions). */
extern double max_cut_non_bonded;
/** For the warmup you can cap the singularity of the Lennard-Jones
potential at r=0. look into the warmup documentation for more
details (who wants to wite that?).*/
extern double lj_force_cap;
/** For the warmup you can cap the singularity of the directionnal LJ
potential at r=0. look into the warmup documentation for more
details (who wants to write that?).*/
extern double ljangle_force_cap;
/** For the warmup you can cap the singularity of the Morse
potential at r=0. look into the warmup documentation for more
details (who wants to wite that?).*/
extern double morse_force_cap;
/** For warm up integration, the maximum force between any two particles
interacting via Buckingham potential can be set and this magnitude of max
force is stored in buck_force_cap*/
extern double buck_force_cap;
/** For the warmup you can cap any tabulated potential at the value
tab_force_cap. This works for most common potentials where a
singularity in the force occurs at small separations. If you have
more specific requirements calculate a separate lookup table for
each stage of the warm up.
\note If the maximum value of the tabulated force at small
separations is less than the force cap then a warning will be
issued since the user should provide tabulated values in the range
where particle interactions are expected. Even so the program
will still run and a linear extrapolation will be used at small
separations until the force reaches the capped value or until zero
separation */
extern double tab_force_cap;
#ifdef CONSTRAINTS
/** numnber of constraints. */
extern int n_constraints;
/** field containing constraints. */
extern Constraint *constraints;
#endif
/**Switch for nonbonded interaction exclusion*/
extern int ia_excl;
/************************************************
* exported functions
************************************************/
/** Function for initializing force and energy tables */
void force_and_energy_tables_init();
/** Implementation of the tcl command \ref tcl_inter. This function
allows the interaction parameters to be modified.
*/
int inter(ClientData data, Tcl_Interp *interp,
int argc, char **argv);
/** Implementation of the Tcl function constraint. This function
allows to set and delete constraints.
*/
int constraint(ClientData _data, Tcl_Interp *interp,
int argc, char **argv);
/** Callback for setmd niatypes. */
int niatypes_callback(Tcl_Interp *interp, void *data);
/** get interaction parameters between particle sorts i and j */
MDINLINE IA_parameters *get_ia_param(int i, int j) {
extern IA_parameters *ia_params;
extern int n_particle_types;
return &ia_params[i*n_particle_types + j];
}
/** Makes sure that ia_params is large enough to cover interactions
for this particle type. The interactions are initialized with values
such that no physical interaction occurs. */
void make_particle_type_exist(int type);
/** Makes sure that \ref bonded_ia_params is large enough to cover the
parameters for the bonded interaction type. Attention: 1: There is
no initialization done here. 2: Use only in connection with
creating new or overwriting old bond types*/
void make_bond_type_exist(int type);
/** This function increases the LOCAL ia_params field
to the given size. Better use
\ref make_particle_type_exist since it takes care of
the other nodes. */
void realloc_ia_params(int nsize);
/** calculates the maximal cutoff of all real space
interactions. these are: bonded, non bonded + real space
electrostatics. The result is stored in the global variable
max_cut. The maximal cutoff of the non-bonded + real space
electrostatic interactions is stored in max_cut_non_bonded. This
value is used in the verlet pair list algorithm (see \ref
verlet.h). */
void calc_maximal_cutoff();
/** check whether all force calculation routines are properly initialized. */
int check_obs_calc_initialized();
/** check if a non bonded interaction is defined */
int checkIfInteraction(IA_parameters *data);
/** check if the types of particles i and j have any non bonded
interaction defined. */
MDINLINE int checkIfParticlesInteract(int i, int j) {
return checkIfInteraction(get_ia_param(i, j));
}
char *get_name_of_bonded_ia(int i);
#endif