ClassicalDynamics ================= ClassicalDynamics is a program (and with it a library) completely written in C++ offering general classical dynamics for students and researchers in science. The main program is an example application of the more general library that can be used in your own scientific software. The library makes use of the object oriented features of C++. The use of the example main program is fairly easy and allows for undergraduates in science to get acquainted with the basics of computer simulations and dynamics. The main program also offers the use of BOINC API. It can make extensive use of OpenGL to render the dynamics, but it can also be compiled without the rendering capabilities and run in text mode. The main program and library have been tested on Linux, Windows and Apple Mac using the GCC, Microsoft Visual studio C++ 5.0, and the devC++ compilers. This document describes the program, with it's use and input format, in section A. In section B the library part of the package will be covered. The library is covered by the LGPL license and the main program is covered by the GPL. Look for more details on the license in section C. Copyright (C) 2005 M.F. Somers, Theoretical Chemistry Department, Leiden University ******************************************************************************** SECTION A: Using the ClassicalDynamics Program ******************************************************************************** 1) Compiling. ------------- Before you can use the program, the code should be compiled. Just go to the directory of the platform you want to compile for and invoke make (Linux and Apple OS-X) or use the solution file to startup the Microsoft Visual Studio. The executable will be put in the main directory tree. For your convenience the supplied binary freeglut libraries are included during the build. If you want to recompile them, the sources of the freeglut library are included in the glut_libs directory. The makefiles are by default configured to compile using GLUT and include the needed OpenGL rendering code. 2) Invoking the executable -- your first run. --------------------------------------------- The program can be started with none or several input options. If no parameters are supplied to the command line, the program automatically seeks the "classical.in" input file within the current directory and will output to standard output and the file "classical.out". If parameters are given, the input and output files can be specified: ClassicalDynamics.exe [inputfile [outputfile [stdoutfile [stderrfile]]]] As an example try the command: ClassicalDynamics.exe examples/hydrogen_gas testoutput This will invoke the program and run the hydrogen_gas example simulation, generating statistical information output to stdout and restart/checkpoint data in the testoutput file. 3) Re-Invoking the executable -- your second run. ------------------------------------------------- If during a run a key is pressed the dynamics is paused and extra statistical information is output to stdout. Pressing a key again continues the simulation. During such a pause a 'snapshot' of the system is output into the output or checkpoint file. One can use the mouse to rotate (by pressing the left button at the same time) the viewpoint of the system. If you use the right mouse button and move your mouse up/down you can zoom in/out the rendered scene! As an example restart the hydrogen_gas simulation and press a key to make a snapshot. Now try to rotate the scene and zoom in a bit. Press a key again to continue the simulation and notice that you can use the mouse during the simulation as well. After that stop the run. Inspect the outputfile and notice that the snapshots/checkpoints are appended to the file each time in a format that is directly readable as input by the program. These snapshots can also be auto-generated in a dynamical run. Read more about that in the documentation about the 'dynamics' statement below. 4) Creating input -- the inputfile format. ------------------------------------------ Inputfiles are plain ASCII text files and can include comments (# character) or blank lines. The inputfiles have the "single statement per line" convention. Morever, in the input file, arbitrary units are used. One is free to use whatever units you like, but you have to be consistent with them yourself! In general an inputfile consist of 6 sections; 1 - first create a small lists of colors you'd like to use for rendering. 2 - then specify what sort of box/cell the dynamics should run in. 3 - then specify one or several particles. 4 - then specify the interactions between particles. 5 - if you want, specify distance constraints between particles. 6 - then specify the dynamics or conformation searcher parameters. section 1 - colors: - - - - - - - - - - For the scene to be rendered, different colors can be used for different particles. In general each particle will have a color index assigned to it. The actual colors corresponding to such an index can be specified by the 'color' statement: color id=## rgb=[#.#,#.#,#.#] in which a color specified through it's Red-Green-Blue components, is assigned to a given index 'id'. The following statements will setup a color table with 4 colors (id 1 to 4) with the colors Red, Green, Blue and White: # start basic color table color id=1 rgb=[1.0,0.0,0.0] # this is red color id=2 rgb=[0.0,1.0,0.0] # this is green color id=3 rgb=[0.0,0.0,1.0] # this is blue color id=4 rgb=[1.0,1.0,1.0] # this will be white section 2 - the box: - - - - - - - - - - When the colors have been specified, the box should be specified with the 'box' statement: box cell=[#.#,#.#,#.#] or box periodic=[#.#,#.#,#.#] If the cell=[#.#,#.#,#.#] format is used, a box of [X,Y,Z]=[#.#,#.#,#.#] dimensions will be used during the dynamics. If one of the particles now moves out-of-the box, the dynamical simulation will automatically stop. Example: box cell=[10,20,30] # use a box of size 10 by 20 by 30 If the periodic=[#.#,#.#,#.#] format was used, a box of [X,Y,Z]=[#.#,#.#,#.#] dimensions will be used in the simulation. However, now periodicity of the system is used according to the minimum-system-image convention. Essentially this means that particles moving out of the box are replicated and moved into the box at the other side. It also implies that particles only interact with the images and particles within a similar box, of equal size, surrounding each particle. Example: box periodic=[1.5,3.75,5.95] # use a minimum-system-image box of size # 1.5 by 3.75 by 5.95 section 3 - the particles: - - - - - - - - - - - - - After having specified the box, particles, interactions and constraints can be specified. It is the convention to first specify particles and afterwards the interactions, and then if needed constraints, corresponding to these particles. Specifying a particle can be done using the 'particle' statement: particle id=## c=# m=#.# q=#.# r=#.# p=[#.#,#.#,#.#] x=[#.#,#.#,#.#] This will create a particle with id number 'id', radius 'r', mass 'm', charge 'q', momentum 'p' at position 'x' with color as specified by color index 'c'. Both the momenta and the position are relative to the center of the dynamical box and are given in cartesian coordinates. Specifying all those parameter per particle is not so efficient. Luckily there is a shorthand convention in the program. The last specified 'c', 'r', 'm', 'p' fields become automatically default values for subsequent particle statements, and the 'id' value is automatically increased. For example, specifying: particle id=1 c=1 m=1.0 q=0.0 r=1.0 p=[0.0,0.0,0.0] x=[1.0,1.0,1.0] particle m=2.0 x=[2.0,1.0,-5.0] c=2 particle q=-3.0 r=2.0 x=[3.0,1.0,-2.0] c=0 particle x=[-3.0,-4.0,3.0] id=10 particle x=[3.1,14.0,3.0] p=[1.0,0.0,0.0] is equivalent to specifying: particle id=1 c=1 m=1.0 q=0.0 r=1.0 p=[0.0,0.0,0.0] x=[1.0,1.0,1.0] particle id=2 c=2 m=2.0 q=0.0 r=1.0 p=[0.0,0.0,0.0] x=[2.0,1.0,-5.0] particle id=3 c=0 m=2.0 q=-3.0 r=2.0 p=[0.0,0.0,0.0] x=[3.0,1.0,-2.0] particle id=10 c=0 m=2.0 q=-3.0 r=2.0 p=[0.0,0.0,0.0] x=[-3.0,-4.0,3.0] particle id=11 c=0 m=2.0 q=-3.0 r=2.0 p=[1.0,0.0,0.0] x=[3.1,14.0,3.0] One can also use spherical coordinates to specify the position or the momentum of the particle. Specify s=[r,theta,phi] instead of the x=[X,Y,Z] to use spherical coordinates in setting up the system. Internally the program uses still cartesian coordinates and the s=[r,theta,phi] coordinates are therefore simply transformed to cartesian coordinates. The same applies to specifying the momenta of particles. Instead of p=[Px,Py,Pz] one can specify ps=[Pr,Ptheta,Pphi]. It is obviously not allowed to have a x=[X,Y,Z] and a s=[r,theta,phi] field in the same particle statement. This also applies for the p=[Px,Py,Pz] and the ps=[Pr,Ptheta,Pphi] fields. section 4 - the interactions: - - - - - - - - - - - - - - - After having specified some particles, interactions between those particles can be defined using the 'interaction' statement. Several types of interactions can be defined between particles: 1) gravitational interactions between two or more particles 2) coulomb interactions between two or more particles 3) 12-6 lennard-jones interactions between two or more particles 4) morse interaction between two or more particles 5) rydberg interaction between two or more particles 6) harmonic stretch interactions between two or more particles 7) harmonic bending interactions between three or more particles 8) periodic torsional interactions between four or more particles In the current program interactions automatically include the minimum- system-image convention when using periodicity through the box statement. It is not yet possible to introduce 'cut-offs' in the interactions within a minimum-system-image-box. This is planned to be implemented at a later stage. When describing an interaction between particles, a so-called 'particle-list' is used and several parameters can be supplied, depending on the specific type of interaction. The particle list is the list of defined particle id's enclosed in {} brackets: { 1 2 3 } corresponds to particles 1 2 and 3. { 1 2 5 6 to 10 } corresponds to particles 1, 2, 5, and 6 to 10. { all } corresponds to all particles defined so far. 1) Gravitational interaction syntax: interaction gravity f=#.# { 1 2 } Creates a mass dependant gravitational interaction between particles 1 and 2 with the gravitational constant equal to 'f'; V(r12)= f * m1 * m2 / r12, with m1, the mass of particle 1, m2 the mass of particle 2, r12 the distance between particle 1 and 2. If more than two particles are specified in the list, all possible combinations are automatically included and a corresponding 'pair force-field' is used. 2) Coulomb interaction syntax: interaction coulomb f=#.# { 1 2 } Creates a charge dependant coulombic interaction between particles 1 and 2 with prefactor 'f'; V(r12)= f * q1 * q2 / r12, with q1, the charge of particle 1, q2 the charge of particle 2, r12 the distance between particle 1 and 2. If more than two particles are specified in the list, again all possible combinations are automatically included and a corresponding 'pair force-field' is used. 3) 12-6 Lennard-Jones interaction syntax: interaction lennardjones f=#.# r0=#.# { 1 2 } Creates a 12-6 Lennard-Jones interaction between particles 1 and 2; V(r12)= f * ( (r0/r12)^12 - 2 * (r0/r12)^6 ), with r12 the distance between particle 1 and 2, 'f' the prefactor and strength of the interaction (the well- depth) which has an equilibrium distance of 'r0' (well-position). Again if more than two particles are specified in the list, all combinations are considered. 4) Morse interaction syntax: interaction morse f=#.# r0=#.# exp=#.# { 1 2 } Creates a Morse interaction, V(r12) = f * ( 1 - e^(-exp*(r12-r0)) )^2 - f between particles 1 and 2. 'f' is the well-depth, 'r0' the equilibrium distance, 'exp' the exponential parameter and r12 the distance between particle 1 and 2. If more than two particles are specified, all possible combinations are automatically included. 5) Rydberg interaction syntax: interaction rydberg f=#.# r0=#.# exp=#.# A1=#.# A2=#.# A3=#.# { 1 2 } Creates a Rydberg type of interaction between particle 1 and 2: V(r12)=f-f*(1+A1*(r12-r0)+A2*(r12-r0)^2+A3*(r12-r0)^3) * e^(-exp*(r12-r0)). Here 'A1' to 'A3' are parameters used in fitting, exp the exponential parameter, 'r0' the equilibrium distance, r12 the distance between particle 1 and 2 and 'f' the well-depth with it's minimum V(r12=r0)=0. 6) Harmonic stretch interaction syntax: interaction harmonic f=#.# r0=#.# { 1 2 } Creates a harmonic interaction between particles 1 and 2; V(r12)= 0.5 * f * ( r12 - r0 )^2, with r12 the distance between particle 1 and 2, 'f' the corresponding force-constant and 'r0' the equilibrium distance. If more than two particles are specified in the list, all combinations are considered. 7) Harmonic bending interaction syntax: interaction bending f=#.# deg=#.# { 1 2 3 } Creates a bending interaction between particles 1, 2 and 3; V(angle123)= 0.5 * fcos * ( cos(angle123) - cos(deg) )^2. The equilibrium angle 1-2-3 angle123 is specified with 'deg' in degrees, or in radians when 'rad' is used instead of 'deg'. The strength is either given by 'f' or by 'fcos'. Internally the interaction is defined as a function of the cosine of the angel 1-2-3 and uses the 'fcos' parameter. However, often only parameters of a quadratic function of the angle directly is known. In such a case 'f' can be used and internally the 'f' parameter in translated into an 'fcos' parameter in such a way that the taylor-expansions of both possible dependencies is equivalent to second order. If more than three particles are given, i.e. { 1 2 3 4 5 ... }, the interaction is automatically defined for the subsequent angles, 1-2-3, 2-3-4, 3-4-5 and so on. 8) Periodic torsional interaction syntax: interaction torsional n=# f=#.# deg=#.# { 1 2 3 4 } Creates a torsional angle interaction as a function of the angle being the angle between the two planes defined by particles (1,2,3) and (2,3,4); V(angle1234)= f * ( 1 - cos( N * (angle1234-deg) ) ). The interaction is periodic with a period of 'N', an equilibrium angle of 'deg' degrees or 'rad' radians and with an amplitude of 'f'. If more than four particles are specified, i.e. { 1 2 3 4 5 6 ... }, subsequent torsional angle interactions, (1,2,3)-(2,3,4), (2,3,4)-(3,4,5), (3,4,5)-(4,5,6) and so on, are automatically included. section 5 - the constraints: - - - - - - - - - - - - - - When the particles and their interactions have been specified, constrainst can be supplied. This is done through the "constrain" statement. Currently only distance constraints have been implemented, so the only possible constraint type is "distance" and therefore the constraint statement has the following format: constrain distance [r=#.#] [maxstep=#] [error=#.#] { 1 2 3 ... } This statement defines a constraint on the distance between particle 1 and particle 2, between particle 2 and 3 and so forth. All these distances can be set fixed to a distance of 'r'. If 'r' has not been specified, the constraint will be ignored unless also a steepest descent conformation search will be perfomed in the same run. In the later case, the constraint distance will be set to the found distance between the first two specified particles of the particle list. The particle list has very much the same format as is used in interactions, except for the 'all' specifier which is not allowed in the list for constraints. The 'maxstep' parameter sets the upperbound on the number of steps in the constraint algorithm SHAKE (steepest descent) or RATTLE (velocity verlet). By default the maximum has been set to 1000 itterations. The algorithm tries to find the correct positions of the particles within a given error. By default the error is set to 1.0E-12, but that can be changed using the 'error' field. section 6 - the dynamics: - - - - - - - - - - - - - After having specified particles, the corresponding interactions between them, and the possible constraints, several statements can be used to start the real dynamics on the system: 1) conformation searcher 2) temperature initializing and control 3) running a trajectory These statements do not need to be present in the inputfile, but leaving all of them would make an inputfile correspond to a rather boring dynamical simulation ;-). 1) The conformation searcher; The 'conformation' statement is used to start the steepest descent searcher. This is an algorithm that tries to find the stable conformation of the system. It will move the particles from their current positions to a situation (conformation) corresponding to a lower potential energy of the system. If the potential energy is at it's lowest and if subsequent movements only increase the potential energy of the system the algorithm stops. The syntax: conformation n=### error=#.# [maxstep=#.#] 'n' equals the maximum number of subsequent steps, with a maximum step size of 'maxstep' per particle, the algorithm will do. If the minimum of the potential energy has been found within the given relative 'error' parameter, within 'n' steps, the algorithm also stops. If the 'conformation' statement is present in the inputfile, the steepest descent algorithm is performed on the system, as is specified in the inputfile, but before any temperature scalings or trajectory calculations are performed. Moreover, when the conformation searcher algorithm has stopped, a snapshot / checkpoint will be output to the outputfile and some extra information will be displayed on stdout. If no conformation search is needed before a trajectory is to be performed, one can leave out the statement. When distance constraints have been specified with fixed distances the steepest descent algorithm applies the SHAKE algorithm after each step to meet the supplied constraints. If no fixed distance has been specified in the distance constraint statement, the particular constraint distance will be set to the distance between the particles found after minimization with the steepest descent algorithm. The constraint is thus not active during a conformation search, but will be initialized afterwards, before doing any dynamics. 2) The temperature control statement; After having performed a conformational search, or not, the particles can be assigned an initial temperature, or the temperature control can be specified. Syntax: temperature k=#.# initial=#.# or temperature k=#.# constant=#.# [tau=#.#] [rmf=#.#] [cmrf=#.#] or temperature k=#.# constant=#.# gamma=#.# [tau=#.#] [cmrf=#.#] Here 'k' is the Boltzmann constant, 'initial' the temperature to which the system should be initialized to the start of a simulation, or 'constant' the temperature the system should be kept at during the simulation. If no temperature init or control needed, one can just leave out the statement. Only when the start time of the simulation (see the dynamics statement below) is zero, the initial temperature is set. The velocities of the particles will be randomly chosen, the center of mass velocity will be removed and finallt the velocities will be rescalled to the desired initial temperature. When doing a constant temperature calculation, the velocities will be rescaled each timestep or a Langevin type of dynamics is performed. See below. When the 'tau' coupling parameter is supplied, the rescaling will be according to the 'weak-coupling' method of H.J.C Berendsen et. al. [JCP 81 (1984) 3684-3690]. In an update (april 2007), the 'rmf' option was added. The 'rmf' is the so-called 'random momentum factor'. When this factor is set to a non-zero positive number, at each temperature rescaling step, a random momentum vector is added to the current momentum of each of the individual particles. The size of this random momentum vector is given as a fraction of the original momentum of the particle. The 'rmf' option cannot be specified if the 'gamma' parameter was also set. In an update (may 2007), the 'gamma' and 'cmrf' options were added. If the 'gamma' Langevin friction coefficient has been supplied, the Langevin equation of motion is used and the system will be coupled to a thermal bath accordingly. The 'gamma' value should be positive and smaller than 1. The Langevin equation of motion includes a friction term and a random force Fr: F = m * a = -dV/dr - gamma*v*m + sqrt(6*m*gamma*k*T/dt)*Fr When doing Langevin dynamics, the 'rmf' parameter is not allowed. It is possible to specify a 'tau' parameter to allow for an extra Berendsen type of velocity rescaling during the Langevin dynamics. By setting the the 'tau' smaller or equal to the dynamics timestep, the normal velocity rescaling algorithm is regained. The new 'cmrf' option specifies what fraction of the center of mass velocity should be removed at each time step during NVT simulations. This parameter should be between zero and one. 3) The start of a trajectory; After the temperature control has been specified, if needed, the trajectory can be started with the 'dynamics' statement: dynamics dt=#.# tend=#.# t=#.# [error=#.#] [snapshots=#] or dynamics dt=#.# t=#.# n=# [error=#.#] [snapshots=#] The trajectory will run from start time 't', with a maximum timestep of 'dt', until 'tend' has been reached or, when using a single cell dynamical box, one of the particles moves out of the box. If the 'n' instead of 'tend' was specified, only 'n' timesteps will be propagated. If the temperature is not rescaled each timestep (no constant temperature dynamics), the timestep 'dt' is automatically adjusted to keep the total energy conservation within the given relative 'error'. The timestep will however, never be taken bigger than 'dt'. Currently the 'velocity verlet' integrator is used. This integrator calculates the positions and momenta of the particles to 4th order accuracy for each time. The velocities and positions do not lag half a timestep like in the 'leap-frog' method and does not calculate the momenta only to 2nd order like the normal 'verlet' method does. The 'snapshots' parameter specifies howmany automatic snapshots the program should generate within the 't' to 'tend' range. At each snapshot the full system is dumped into the output file and extra information is printed to stdout. The output generated by a snapshot can be used to restart the run from that point onwards directly. The output to the output file is in the same format as the input file should be. When distance constraints have been set, either by the user through specifying fixed distances in the constraint statement, or by the steepest descent searcher, if a 'conformation' statement was present in the input, the RATTLE algorithm is used during the calculation of the trajectory. The RATTLE algorithm makes sure that the nice properties of the velocity verlet method are retained, despite the inclusion of constraints in the dynamics. section 7 - extra statements for convenience: - - - - - - - - - - - - - - - - - - - - - - - When dealing with larger collections of particles, or blocks of particles, it can be convenient to be able to translate, rotate or introduce a scale, for a given block of position vectors. This can be accomplished through the 'scale', 'translate' and 'rotate' statements: scale [X=#.#] [Y=#.#] [Z=#.#] translate [X=#.#] [Y=#.#] [Z=#.#] rotate deg=#.# axis=[#.#,#.#,#.#] At the start of the program the scale values are set to X=Y=Z=1, the translate values to X=Y=Z=0 and the rotate values to deg=0 with axis=[0,0,0]. By specifying different values, transformations to all subsequent position vectors (not momenta!) are applied in the order of first rotating then translating and finally rescaling. Look in the example directory for examples of the use of these statement. Specifically look at the water simulation in which several molecules are defined and placed by these statements. ******************************************************************************** SECTION B: Using the ClassicalDynamics C++ Library ******************************************************************************** The ClassicalDynamics C++ library offers a framework to implement extendable dynamics into your own program. The library consists of several .h and .cpp files located in the sources subdirectory. The library is covered by the LGPL license, in contrary to the main example program, described in the previous section, which is covered by the GPL. 1) The five base classes. ------------------------- The library is built onto five basic classes, the 'Particle' class, the 'Interaction' class, the 'Constraint' class, the 'ClassicalStates' class and derived from the later, the 'ClassicalDynamics' class. All of these classes can be found in the ClassicalDynamics.h header, the Particle.cpp, the Interaction.cpp, the ClassicalStates.cpp and the ClassicalDynamics.cpp files located in the sources subdirectory. section 1 - the 'Particle' class: - - - - - - - - - - - - - - - - - The Particle class consists of several datamembers representing a particle. Apart form the obvious members, like 'ID', 'Mass', 'Charge', 'Position' and so on, it also contains two protected members 'Interactions' and 'DynamicsContext'. The 'Interactions' member is a STL vector of pointers to 'Interaction' classes. This vector represents all the interactions the particle feels. This list is used to calculate for example the force on a particle through the 'TheTotalForce' method, or the potential energy of the particle through the 'ThePotentialEnergy' method. One can include an 'Interaction' through the 'IncludeInteraction' method. One can also query the list with the 'IsInteractionIncluded' and 'NrOfInteractionsIncluded' methods. With the 'ExcludeInteraction' method an entry in the list is removed, but the corresponding 'Interaction' class, to which the 'Interactions' vector entry points, is not destructed. Use the 'DeleteInteraction' method to exclude an interaction from the list and destruct (free the memory of it) the pointed 'Interaction' class. The 'Interaction' class also contains an STL vector with pointers to 'Particle' classes included into the 'Interaction'. If an 'Interaction' is excluded from a 'Particle's' 'Interactions' list (using the ExcludeInteraction method), the correct 'Particle' class is also removed from the 'Interaction's' 'Particles' list. Therefore the 'Interactions' list is a protected data member and should only be accessed through the methods the 'Particle' class offers. Another consequence of this philosophy is clear now, 'Particles' can only interact with each other through an 'Interaction' class. So that when calculating the force on a given particle, only those interactions that are associated with a 'Particle' are used to calculate the force on the 'Particle'. Other 'Particles' within a dynamical contex or system, are not needed for this! The 'DynamicsContext' member of the 'Particle' class, is a pointer to a 'ClassicalStates' class. If the 'Particle' is a member of a classical system, a 'ClassicalStates' class, a reference is used to that system. Why this is needed can be understood when considering a dynamical simulation. When a trajectory is calculated, the system evolves from one 'state' into the other 'state'. Probably two or more states are needed, depending on the integrator one would like to use. However, this implies that, within a system, for each state of the context, copy 'Particle' exists in the 'ClassicalStates' class. If now the 'Particle' class 'ExcludeInteraction' method is called within this framework, the 'Interaction' should also be excluded from all other copy 'Particles', belonging to the different states of the 'ClassicalStates' class of the system. section 2 - the 'Interaction' class: - - - - - - - - - - - - - - - - - - The 'Interaction' class also consist of data members and methods. Very similar to the 'Particle' class, the data member 'Particles' of this class is an STL vector of pointers to 'Particle' classes. The 'Particles' list correspond to all the particles that interact with each other through this 'Interaction'. Methods like 'NrOfParticlesIncluded', 'IncludeParticle', 'IsParticleIncluded', 'ExcludeParticle' and 'DeleteParticle' offer the same functionality as the very similar methods described in the previous section. Again, when a 'Particle' is excluded from the 'Interaction's' 'Particles' list, the 'Interaction' is also excluded from the corresponding 'Particle's' 'Interactions' vector. Moreover, when a 'Particle' is excluded from an 'Interaction', within a 'ClassicalStates' framework, the data member 'DynamicsContext' is used like in the 'Particle' case. That means, copy particles within the same dynamical 'ClassicalStates' context, but belonging to different states of that dynamical context, are automatically also excluded form the copy 'Interaction's' 'Particles' list. At first, it all seems rather complex, but after carefull study, you will understand why this is setup like it is and why this is a very powerfull way of doing things. Any 'Particle' can automatically work with any 'Interaction' in whatever classical context the 'ClassicalStates' class corresponds to! section 3 - the 'Constraint' class: - - - - - - - - - - - - - - - - - - This base class serves as a stub for constraints. It offers a framework to list particles together that are part of a particular constraint. It has the methods 'IncludeParticle', 'ExcludeParticle', 'IsParticleIncluded', 'GetParticleId', 'ExcludeAllParticles', 'NrOfParticlesIncluded' and 'GetParticleListSize' to manage a list of particles (through their IDs) for constraints in general. The class also offers members like 'Type', 'ParticleIDs' and 'ParticleIndices' for other classes, capable of dealing with constraints, to interact. Derived from this base class is the 'DistanceConstraint' class, which in turn implements a distance constraint between several particles. The DistanceConstraint class therefore has the extra members 'Distance', 'Tolerance' and 'MaxItterations'. For example a class like 'NewtonDynamics' or 'SteepestDescent', allow for a list of constraints to be given as an argument to some of their methods. More information on this will given below. section 4 - the 'ClassicalStates' class: - - - - - - - - - - - - - - - - - - - - The 'ClassicalStates' class offers the possibility to define a classical context for 'Particles' and 'Interactions'. The 'ClassicalStates' class therefore has a list of 'ParticleList's' pointers and a list of 'InteractionList's' pointers. Again these two data members, 'TheParticleListList' and 'TheInteractionListList', are protected and should only be accessed through the 'ClassicalStates' methods acting upon them. The 'ClassicalStates' class offers the 'IncludeParticle', 'IncludeInteraction', 'IsParticleIncluded', 'IsInteractionIncluded', 'DeleteParticle', 'DeleteInteraction', 'ExcludeParticle', 'ExcludeInteraction', 'NrOfParticlesIncluded, 'NrOfInteractionsIncluded', 'DeleteAllParticles' and 'DeleteAllInteractions' for this. There is however a difference; when including a 'Particle', copies are made of the 'Particle' to be included in the different states. The same applies to 'Interactions' being included. Be sure to destruct included 'Particles' and 'Interactions' yourself, or suffer from a memory leak! Moreover, when a 'Particle', which is included into an 'Interaction', is included, all the needed 'Interactions' are also copied and included into the 'ClassicalStates' class states. Because the 'ClassicalStates' class makes copies, for each state of the dynamical context, of 'Particles' and 'Interactions', the inclusion of a 'Particle' into an 'Interaction' that's part of the 'ClassicalStates' needs to be done through the 'ClassicalStates' methods. The same applies to including an 'Interaction' into a 'Particle' that's already part of a 'ClassicalStates'. The 'ClassicalStates' class offers therefore the following methods; 'IncludeParticleInInteraction', 'IncludeInteractionInParticle', 'ExcludeInteractionFromParticle', 'ExcludeParticleFromInteraction', 'ExcludeAllInteractionsFromParticle', 'ExcludeAllParticlesFromInteraction', 'DeleteAllInteractionsOfParticle' and 'DeleteAllParticlesOfInteraction'. Apart from methods for handling the various lists within the 'ClassicalStates' class, the class also offers some methods to cycle through the different states of the dynamical context, to retrieve a pointer to a 'Particle' within a specific state, to retrieve a pointer to an 'Interaction' within a specific state, to retrieve a pointer to the list of pointers to 'Particles' of a state, to retrieve a pointer to a list of pointers 'Interactions' of a state, the sizes of these lists and the actual number of states present within the 'ClassicalStates' class; 'CycleTheStates', 'GetParticlePointer', 'GetInteractionPointer', 'GetStateParticlesPointer', 'GetStateInteractionsPointer', 'GetParticleListSize', 'GetInteractionListSize' and 'GetStateListSize'. section 5 - the 'ClassicalDynamics' class: - - - - - - - - - - - - - - - - - - - - - The last base class of the library is the 'ClassicalDynamics' class, which is derived from the 'ClassicalStates' class. This class offers the use of different integrators to be implemented and to run a trajectory. The class has no extra data members, and offers the methods 'GetFinalStateParticlesPointer', 'GetFinalStateInteractionsPointer', 'Integrator', 'IntegratorConservationCheck', 'RunTrajectory' and 'TrajectoryContinuationCheck'. The 'ClassicalDynamics' class uses the state corresponding to number 1, to be the final state and therefore has at least two (0 and 1) states within it. It can contain more states, but it has at least 2! The state associated with number 0 is the 'future' state. The state with number 1 is the 'current' state. The 'RunTrajectory' method used the 'Integrator' to compute the 'future' (0) state from the 'current' and 'past' (1 and upward to n) states. If the computation of the 'future' (0) state was correct, according to the overloaded 'IntegratorConservationCheck' method, the states are cycled so that the 'future' (0) state becomes the 'current' (1) state. The 'past' states are also cycled through in the same way so that the now obsolete 'past' (n) state will be used as the new 'future' (0) state. When the 'RunTrajectory' method has successfully integrated such a timestep, the overloaded 'TrajectoryContinuationCheck' method is called to check if the trajectory should continue or not. With this implementation, merely the overloading of the 'Integrator', the 'IntegratorConservationCheck' and the 'TrajectoryContinuationCheck' methods, is enough to allow classical dynamics with different sorts of integrators and or conditions for the trajectories. 2) The derived interaction classes. ----------------------------------- The 'Interaction' class is used a base class for essentially four different types of other extendable interaction classes; the 'FiniteDifferenceInteraction' class, the 'PointPairInteraction' class, the 'BendingInteraction' class and the 'TorsionalInteraction' class. The 'FiniteDifferenceInteraction' class is a stub for a general interaction, for which no analytical force can be computed. It uses the three-point finite difference method in three dimensions on the 'InteractionPotential' method, inherited from the 'Interaction' class, which should be overloaded. The 'PointPairInteraction' class is a derived class from 'Interaction', from which again other more specific interaction classes are derived. The base 'PointPairInteraction' class is a stub for pair potentials. Only the distance between two particles is needed to construct the full potential or force field. The base class offers two new methods to be overloaded; 'PointPairPotential' and 'PointPairForce'. Both new methods accept two particles and should return the potential of the first particle, or the force excerted on the first particle, due to the second particle. From this class another base class is derived, 'FiniteDifferencePointPairInteraction', in order to implement point pair potentials for which an analytical expression of the force is unavailable. Now a three point finite difference method is used to calculate the point pair force on a 'Particle'. Other classes that are also derived from 'PointPairInteraction' include the 'CoulombInteraction' class, the 'GravitationalInteraction' class, the 'HarmonicInteraction' class, the 'LennardJonesInteraction' class, the 'MorseInteraction' class and the 'RydbergInteraction' class. These derived classes all overload the 'PointPairPotential' and 'PointPairForce' methods of 'PointPairInteraction', for for these interactions, an analytical expression for the force is available. The third class that is derived from the base 'Interaction' class, is the 'BendingInteraction'. This class is again a base class which offers the 'BendingPotential', the 'BendingForce' and the 'DerivativeOfBendingPotential' methods to be overloaded. The default 'BendingForce' method uses the 'DerivativeOfBendingPotential' to calculate the force as a function of the cosine of the angle betwen three or more 'Particles'. If more 'Particles' are included into the class, the subsequent angles between the 'Particles', as described in section A of this document, are included into the total potential of the interaction and in the force on one of the 'Particles'. The derived class 'FiniteDifferenceBendingInteraction' overloads the 'DerivativeOfBendingPotential' method to implement the three point finite difference method if no analytical expression of the derivative of the bending potential is known. This derived class is however a stub. The 'HarmonicBendingInteraction' class is also derived from the 'BendingInteraction' base class and implements, through the overloading of the 'BendingPotential' and 'DerivativeOfBendingPotential' methods, a harmonic potential as function of the cosine of the angle defined by three 'Particles'. A similar construction is used for the fourth class derived form the base 'Interaction' class, the 'TorsionalInteraction' base class. This base class offers a framework for torsional interaction potentials between four or more 'Particles', through the derived 'PeriodicTorsionalInteraction' and the 'FiniteDifferenceTorsionalInteraction' class. The later class overloads the 'DerivativeOfTorsionalPotential' method of the 'TorsionalInteraction' class to use a three point finite difference algorithm. The other method the base class 'TorsionalInteraction' offers to overload are the 'TorsionalPotential' and the 'TorsionalForce' methods. Again the default 'TorsionalForce' uses the 'DerivativeOfTorsionalPotential' method to compute the force on one of four 'Particles' as a function of the torsional angle between them. Here is a schema of how all classes are related: Interaction (stub) | +-- FiniteDifferenceInteraction (stub) | +-- PointPairInteraction (stub) | | | +-- FiniteDifferencePointPairInteraction (stub) | | | +-- CoulombInteraction | | | +-- GravitationalInteraction | | | +-- HarmonicInteraction | | | +-- LennardJonesInteraction | | | +-- MorseInteraction | | | +-- RydbergInteraction | +-- BendingInteraction (stub) | | | +-- FiniteDifferenceBendingInteraction (stub) | | | +-- HarmonicBendingInteraction | +-- TorsionalInteraction (stub) | +-- FiniteDifferenceTorsionalInteraction (stub) | +-- PeriodicTorsionalInteraction Clearly the framework allows for even the most complicated types of interactions. Just pick the right base class to derive the special case from! 3) The derived NewtonDynamics class and the SteepestDescent class. ------------------------------------------------------------------ The 'ClassicalDynamics' base class is used to derive the 'NewtonDynamics' base class. The 'NewtonDynamics' class implements the velocity verlet integrator through overloading the base 'Integrator', 'IntegratorConservationCheck' and 'TrajectoryContinuationCheck' methods. The default 'IntegratorConservationCheck' and 'TrajectoryContinuationCheck' methods of 'NewtonDynamics' return 1, signaling the integration went successful and the trajectory should continue. Overload these method to change that behavior. The second class derived from the 'ClassicalStates' class is the 'SteepestDescent' class. The 'SteepestDecent' class implements the steepest descent algorithm to be applied to a collection of particles in a given forcefield, specified by 'Interactions' between them. The class offers the 'FindMinima', the 'StepSteepestDescent' and the 'RelativeTotalInteractionPotentialChange' methods which can be overloaded. The default behavior of these methods, however, are for 'FindMinima' to use the 'StepSteepestDescent' to step towards the local minimum in the multi-dimensional potential energy surface, using 'RelativeTotalInteractionPotentialChange' to check for convergence. The 'FindMinima' method uses smaller steps when the forces are big. It also automatically makes the steps bigger, when the forces are smaller. Here is a schema of how all these classes are related: ClassicalStates (stub) | +-- SteepestDescent | +-- ClassicalDynamics (stub) | +-- NewtonDynamics | +-- LangevinDynamics The main program example, again uses a class derived from the 'NewtonDynamics' class to implement the calculation of statistical properties during the simulation. It can however use any class derived from the ClassicalDynamics class, implementing different integrators, for example! The NewtonDynamics and the SteepestDescent classes also are capable of dealing with distance constraints. A list of constraints can be supplied to the 'Integrator', 'RunTrajectory', 'FindMinima' and 'StepSteepestDescent' methods. The distance constraints in NewtonDynamics are implemented through the RATTLE algorithm (which is a velocity verlet verion of SHAKE), whereas in the 'SteepestDescent' class, normal SHAKE is implemented. Other constraints than 'DistanceConstraint', that might be added in the future and can therefore be present in a constraint list, are ignored by these two classes. In the LangevinDynamics class, the Newton equation of motion is used with frictional terms and random forces added. The LangevinDynamics class is therefore similar in setup and implementation as NewtonDynamics is. If a friction factor 'Gamma' of 0.0 is used, the LangevinDynamics is equivalent to NewtonDynamics class. The LangevinDynamics class was added in May 2007. Information on the RATTLE and SHAKE algorithm's can be found in the sources and in "H. C. Andersen, Journal of Computational Physics 52, 24-34 (1983)". 4) The other functions and constants. ------------------------------------- The library offers a few functions to retrieve basic information on 'Particles' and 'Interactions': 'CentreOfMass', 'DipoleMoment', 'TotalPotentialEnergy', 'TotalKineticEnergy', 'TotalMass', 'TotalCharge' and 'Virial'. Their meaning are obvious and can be found in ClassicalDynamics.h. Apart from those functions, the Vector.h header offers not only the use of 3D cartesian vectors, but also the functions; 'Angle', 'Cross', 'SphericalToCartesian', 'CartesianToSpherical', 'RotateVectorAroundAxis', 'MimimumSystemImageConventionVector', 'ARandomNormalizedVector', the correct operators to calculate with Vectors and to perform IO on streams for Vectors. Furthermore, the ReadInput.h and ReadInput.cpp files offer the framework of creating, reading and writing inputfiles in a format that was described in the first section of this document. The ReadInput.h header defines three new data types; 'Color', 'ColorList' and 'InputFileDataType'. The later is a data structure that's built when input is read, or used when one wants to create an inputfile. Apart form these data types, the functions 'ReadDataFromFile', 'WriteDataToFile', 'DeleteInteractionList', 'DeleteParticleList', 'DeleteColorList', 'DeleteTheData', 'AddParticleToList', 'AddInteractionToList', 'AddColorToList', 'DeleteInteractionFromList', 'DeleteParticleFromList' and 'DeleteColorFromList'. The ReadInput routines have been written in a very extendable fashion. Internally tables, which contain pointers to specialized functions to read, parse and write strings, according to the described format, are used. Extending the ReadInput routines is therefore rather straightforward, for new types of interactions or other statements. Constants are found in the Constants.h and Global.h files. The Global.h file specifies the different types of classes defined in the library. Whenever the library is extended with an interaction, a new type is recorded into this file. ******************************************************************************** SECTION C: License stuff ******************************************************************************** The library software is licensed according to the LGPL: http://www.gnu.org/licenses/lgpl.txt a copy of that text is located in LIBRARY_LICENSE.txt The example main program is licensed according to the GPL: http://www.gnu.org/licenses/gpl.txt Copyright (C) 2006 M.F. Somers, Theoretical Chemistry Department, Leiden University -------------------------------------------------------------------------------- The sources directory contains both the code of the library and the code of an example main program. One is allowed to used the library part according to the LGPL. The main program consists of Main.cpp and *.inc and can be used to test the library. The main program, and the corresponding compile scripts or make files, are all covered by the GPL license. --------------------------------------------------------------------------------