Home

Moldy User's Manual

1.
2. Reciprocal space term 2 R la 2 24 iJlMDi m 8 2 Molecular Virial Correction Charged system term The true internal stress is the ensemble average Tim where T7 and nh are the short range force and kinetic contributions respectively enters into the Parrinello Rahman equations of motion see section 2 7 The term marked Molecular Virial Correction in equation 2 24 is the difference between the site site virial Y f ri and the molecular virial F and is subtracted after all of the site forces and the molecular forces have been calculated including the short range potential components which are not included in the equations above Though it is not apparent reference B7 Appendix A this term has exactly the same form for all parts of the stress the short range potential and the real and reciprocal space Coulombic parts 2 5 Periodic Boundaries the Link Cell Method The real space part of the Ewald equation 19 and short range potential energy is a sum of contributions over pairs of sites In both cases the interaction decays rapidly with separation which means that only site pairs closer then some cutoff distance need be considered Several methods are available to enumerate site pairs and choose those within the cutoff Most simple MD codes simply loop over all pairs of particles in the MD box and compute the separation r for each If r l
3. off on text mode save write restart print config write save file write text restart Figure 5 3 The main timestep loop in file main c showing the outermost control structures Periodic analysis trajectory dump or output tasks are all shown here including the accumulation of the radial distribution function data which is integrated into the link cell force calculation constraints of equations D 6 b call force calc and ewald to evaluate the potential energy and the forces on all atomic sites The site co ordinates must themselves be computed from the centre of mass co ordinates and the quaternions by function make sites c calculate the molecular centre of mass forces and torques using equations and which are implemented in functions mol force and mol torque At this point the accelerations from the previous timestep are shuffled down from the acc qddot etc arrays in system and species to acco qddoto to make way for the new accelerations computed in this timestep d calculate the internal stress equation 34 and the accelerations of the MD cell matrix from the Parrinello Rahman equations of motion Function rahman implements equation 72 33 calculate the new centre of mass and extended system variable accelerations including the velocity dependent parts which arise in a constant pressure temperature calculation Iterate steps iii v of equations 2 13 u
4. Set up sink and cosk arrays for k 0 a b c Get sin for k ha kb Ic using recursion loop over k vectors Calc pre factors for energy force qsincos Evaluate sin cosk for this k and all sites Eval struct factor and energy term Eval forces on all sites amp add to total Figure 5 6 Evaluation of the reciprocal space parts of the ewald sum 5 3 3 Site Forces Calculation The evaluation of the forces on atomic sites comprises almost all of the run time of any reasonably sized simulation over 95 for systems of only a few hundred atoms Efficiency of execution is therefore a primary consideration in this part of the code It is frequently the case that a straightforward implementation of an algorithm is not optimally fast which regrettably means that the optimized code is not as transparent to read as for other less critical parts of the program The link cell algorithm for the short ranged part of the forces section D 5 is one which scales very well to large systems In its original form 5 the MD cell is divided into a number of subcells and a linked list structure BO lists the molecules or atoms belonging to each subcell The inner loops of the force calculation must therefore use an indirect indexing operation to obtain the co ordinates on
5. the remainder of a line following a 4 symbol and keywords may be entered in upper or lower case For example title Moldy example This is a comment 1 The above blank line is ignored nsteps 1000 step 0 0005 restart file RESTART DAT end The control file ends here On VMS moldy may be defined as a foreign command by moldy mydisk mydir moldy 21 is important that the dumpext utility is found in a directory in the shell search path as it is invoked by several of the other utilities to read dump files 3Some operating systems Unix and MS DOS allow file redirection whereby the standard input is associated with some file This may also be used to supply the control file provided that no command line parameter is given 19 sets the title of the simulation to Moldy example the number of steps to execute to 1000 the timestep to 0 0005 and supplies the name of a restart file to read from It is not necessary to specify all of the parameters on each run Unless it is explicitly assigned a value in the control file each parameter has a default value This is either the value listed in table B I or the case where the simulation is continuing from a restart file the value it had on the previous run see section Parameters are read in sequence from the control file and if one appears more than once only the final instance is used The two most important parameters are step which sets the size of th
6. mode for dumps as well as restart files As yet XDR is not available on every machine Therefore a program called dumpconv is provided which converts dump files to a portable text file format which may be easily moved between machines and back again It is described in appendix 3 8 Constant Stress or Pressure Simulation The control parameter const pressure switches from a constant volume simulation to one of several constant stress or constant pressure ensembles see section D 7 according to its value Available values are 1 The constant stress method of Parrinello and Rahman B9 2 constant pressure method which is the uniform dilation limit of the Parrinello and Rahman dynamics This is convenient if switching between the two but like P R the dynamics is not modularly invariant and does not preserve the virial 3 The improved constant stress method of Cleveland IT and Wentzcovitch b6 4 The Andersen constant pressure 8 which is the uniform dilation limit of Wentzcovitch dynamics The value of the MD cell mass parameter W is given by the control parameter w and the external pressure by pressure Note that because of the different ficticious kinetic energy term in the Lagrangians typical values of W for the Wentz covitch and Andersen ensembles are smaller by a factor of the order of L where L is the MD cell edge than those useful in the Parrinello Rahman and compatible uniform dilation cases At pr
7. 56 As OPS W ot OR SAG A qub E PEA eR BE 56 Aqueous MgCb 5 56 57 Prog 58 58 58 58 59 60 61 62 63 63 64 65 i List of Tables List of Figures 21 The Skew Start method D ample Mo GUU Be eee et td AP c oS D xample output of radial distribution functions llle b torage layout of a Z dimensional pointer based D 2 a Block diagram of initialization function start up d D emaintmestep loop 4s vos m Lao A 9 y E AS RU Per Y 4 4 ow diagram function do step which performsasingletimestep D 4 b ow diagram of function do step which performs a single timestep D e Link Cell short range force calculation valuation of the reciprocal space parts the ewald sum iv Chapter 1 Introduction Moldy is a computer program for performing molecular dynamics simulations of condensed matter It can handle any assembly of rigid polyatomic molecules atoms or ions and any mixture thereof It uses the link cell method to calculate short range forces and the Ewald sum technique to handle long range electrostatic f
8. BT who also managed to permute the components which means that the parameters do not form an ordered number quartet which obeys quaternion algebra Like Allen and Tildesley 2 page 88 we follow the definition of Goldstein I9 pages 143 and 155 Experiments conducted while developing Moldy show that enforcing this constraint significantly decreases the fluctua tions in the total energy Linear molecules are a slightly special case as the moment of inertia about the molecular axis is zero Though there are unique methods to represent this situation Z page 90 Moldy uses a minor modification of the quaternion algorithm that is necessary is a little special case code to avoid dividing by the zero component of inertia in the solution of equation 2 4 and to hold the components of angular velocity and acceleration about the molecular axis to zero This has the considerable advantage of uniform treatment of all kinds of molecules which is convenient when dealing with heterogeneous mixtures 2 22 Integration Algorithms The dynamical equations 2 3 and are integrated using this author s modification 47 of the Beeman algorithm H For atomic systems the accuracy is of the same order as the commonly used Verlet algorithm 55 The accuracy of common integration algorithms was discussed by Rodger 9 who showed that the Beeman algo rithm is the most accurate of the supposedly Verlet equivalent algorithms The Beeman algorithm is the on
9. bit 1 centre of mass velocities quaternion and unit cell matrix derivatives bit 2 centre of mass accelerations quaternion and unit cell matrix second derivatives bit 3 forces torques and stress tensor Items selected are written in the order laid out above Within each set of variables values are ordered primarily by species in the order they appeared in the system specification Within a species ordering is by molecule or atom and at the finest level by x y or z component qo q3 for quaternions Therefore if n is the total number of molecules and is the number with rotational freedom the size of each record is 3n 4n 9 1 if bit 0 is set 3n 4n 9 if bit 1 is set 3n 4n 9 if bit 2 is set 3n 3n 9 if bit 3 is set single precision floating point numbers The header is a copy of a struct dump_t see source file structs h for the format It contains the simulation title and version number the timestep at the beginning of the file the control parameters dump interval and dump level the maximum and actual number of dump records in the file a unique marker actually a timestamp common to all the files in a dump run and the timestamp of any restart file used to start the run It is not possible to dump directly to magnetic tape Moldy must rewind to the beginning of a file to keep the header up to date with the number of dumps in the file as well as extend existing files Neither operation is allowed on a tap
10. for machines lacking specialist libraries as well as for scalar computers 4 2 3 Optimization for Vector Architectures The program should of course be compiled with options specifying vectorization Since highly optimizing and vector izing compilers frequently contain bugs and since some options generate unsafe optimizations it may be necessary to restrict the highest optimization level to those modules which contain critical code To allow the compiler to generate vector code it must be instructed to ignore apparent vector recurrences The reason is that the run time dimensioned arrays necessary to implement such a flexible program must use pointers as their base See any C textbook e g Kernighan and Ritchie D8 Chapter 5 for an explanation of C pointers and arrays Unfortunately this means that the compiler can not determine that each iteration of the loop is independent of the preceding iterations In the jargon of vectorizing compilers there may be a vector dependency or recurrence compiler can be notified that these are not genuine recurrences either globally by use of a command line directive or on a per loop basis using machine specific compiler directives inserted into the source C uses the opposite convention to FORTRAN in storage layout of multidimensional arrays 5Technically the C standard treats arrays passed as formal function parameters as pointers which are permitted to refer to overlapping areas of memory Th
11. is provided in defs h for this purpose This also contains other xalloc macros customized for various data types The complementary function t ree calls the library function ree but also allows tracing of allocation and freeing for debugging purposes All of Moldy s memory management functions are in the source file alloc c Like FORTRAN C only permits the declaration of arrays of size fixed at compile time Unlike FORTRAN C lacks any method of declaring arrays with adjustable innermost dimensions as function formal parameters However through the use of pointer array mapping and dynamic memory allocation variable sized multidimensional arrays may be emulated Z8l 1071 4 pp 20 23 Multidimensional array emulation is done by the function arralloc see Fig ure This takes the size of the atomic data type the number of dimensions and the lower and upper bounds for each dimension as arguments and returns a pointer to the allocated pseudo array This is an array of pointers to an array of pointers to to the data area The C mapping of array notation onto pointer syntax allows this construct to be referenced exactly as if it was a true multidimensional array For example double r r arralloc sizeof double 2 0 m 1 0 n 1 defines 2D array m x n of double so that r i j is a reference to the i j th element in the usual manner This pseudo array may be passed directly to a function e g double doit
12. o r 4 Bexp Cr This is the most general 6 exp form and includes potential of the Born Mayer and Bucking ham forms MCY Aexp Br Cexp Dr generic r Aexp Br C r D r E r F r This is a composite which contains terms of several inverse powers which may be useful for ionic solution simulations In addition to the short range potential electric charges may be specified for each atomic site which interact according to Coulomb s Law See the following section for details 2 4 Ewald Sum The long range Coulomb interactions are handled using the Ewald Sum technique in three dimensions 6 2 p 156 The electrostatic potential of a system of charges is expressed as a sum of short range and long range contributions Each term is written as a series the first real space and the second obtained by Fourier transformation using the periodicity of the MD cell in reciprocal space The expression for the Coulomb energy U is where the daggered summation indicates omission of site pairs i j belonging to the same molecule if n 0 The meaning of the symbols is 1 N erfc a r n qiq j 4 j i 1 rij Real space term i 8 N 2 N 2 y 12 4a Y aicos k ri aisin k ri E0V 2 Reciprocal space term N a 2 1 NS Oral 3
13. p Wh V Wh h 2 37 which replace equation 7 33 Both of the above equations cause the MD cell to dilate or contract uniformly and apply whether or not the initial h matrix corresponds to a cubic cell It is therefore possible to switch between a constant stress and the corresponding constant pressure ensemble during a run This is the only reason to choose equation 2 37 Equation 2 36 selected by const pressure 4 is normally the best option for liquid simulations 2 8 Radial Distribution Functions The radial distribution function or RDF is one of the most important structural quantities characterizing a system particularly for liquids For one component system the RDF is defined as ZU p445 g r EEA 2 38 T ji V r2 0 where we use the angle brackets to denote a spherical average as well as the usual configurational average Allen and Tildesley 2 pp184 185 describe how to evaluate g r from a histogram of pair distances accumulated during a simulation run With the notation that N is the total number of particles b is the number of the histogram bins r is the bin width so that r bdr nnis b is the accumulated number per bin is the number of steps when binning was carried out pond In the case of a molecular system the partial RDF for atoms and f is defined as 20 p 445 1 gog 7 E N N 1 8 r rio rog 2 40 which may be
14. selector format Default is O to half the number of dump records selected with an increment of 1 range of species to be included given in selector format where 0 represents the first species Defaults to all species calculate the VTF instead of the VACF name of optional output file Defaults to standard output If the executable file is named mdvaf then the default calculation is the VACF If it is renamed to mdvtf then the default is the VTF In both cases the output is in the form of a single column of data headed by the name of the species in question The lines down each column correspond to increasing time intervals Similar to msd the time intervals selected using m are taken relative to the dump slices selected using t 62 8 Mdavpos Mdavpos is a utility for calculating the average positions of particles from dump slices recorded during a Moldy run The utility is similar to Mdshak with command line options mdavpos s system specification r restart file d dump file format t dump range h p x o output file where the arguments have the following meanings 5 read a system specification file read a restart file Only one of s or r may be given read configurational data from a dump file given as a prototype name containing a printf format string see section range of records in the dump file specified in selector format to average positions over give the output in SCHA
15. specified in selector format i the increment the initial time fg is increased by during calculation of msd for a given time interval Defaults to 1 g range of species to be included given in selector format where 0 represents the first species Defaults to all species u option to by pass msd calculation and output the connected particle trajectories in the format selected by w w selects output format for trajectory data Set to 0 for format readable by GNUPlot default or 1 for generic format simple matrix X y Z X Y indicate that limits in are to be specified in the given cartesian direction If the direction is in lowercase only particles whose initial positions lie within these limits are included in the msd calculation or trajectory output If the direction is in uppercase the program will only include particles lying outside these limits limits can be imposed in one two or all three directions in any possible combination However if both options for a single direction e g x and X are used priority will be given to the latter option entered on the command line The user is prompted for the limits at run time The default is to include all particles o of optional output file Defaults to standard output The output format for msd calculations consists of four columns corresponding to the x y z components and total msd with increasing time intervals down the columns The time interv
16. to multi million particle systems would be desirable most practical simulations of solids or liquids can be accomplished using systems of a few tens of thousands of particles or fewer 5 4 2 Implementation Two of the design goals are that the serial and parallel versions of Moldy can be built from the same code and that the code is easily portable to a number of communications libraries The C preprocessor macro SPMD is used to conditionally include the parallel code calls to the parallel library interface are protected by conditional compilation and if SPMD is undefined at compile time then the serial version is built These calls are limited to a very few modules main c for the parallel set up do step in accel c for the global summation of the forces virial and potential energy and message in output c for error handling Furthermore the communications library is not called directly but through interface functions defined in file parallel c these are listed in section 3 2 This means that only parallel c need ever be modified if a port to a different communications library is needed as the rest of the code sees a common interface As supplied parallel c contains implementations for the MPI library the TCGMSG library the Oxford BSP library and the Cray SHMEM library for the T3D T3E series machines see section 1 One of these is selected by conditional compilation of parallel c with one of the preprocessor symbols MPI TCGM
17. 8888213e 14 MCY Water Mg2 Cl solution Water 200 1 0 0 0 16 0 0 2 0 7569503 0 0 5858822 1 0 717484 H 2 0 7569503 0 0 5858822 3 0 0 0 2677 0 1 434968 M Magnesium 4 4 0 0 0 24 31 2 Mg2t Chloride 8 56 1 2 2 4 4 111 5 5 S o p 0 1088213 2 455 427 666 3373 47750 0 198855 0 1857 0 28325 5 4 Quartz 152712 961895 760844 836 06 910 408 65 0 35 45 1 0 213 5954 0 546 3 0 0 77 94 Cem co qp NS 233264 253 New values of Mg potl 0 369 Quartz parameters from Van Beest Kramer and Van Santen Physical Review Letters 64 16 1955 1990 Units are eV el chg so time unit 1 0181e 14 Oxygen 384 1 0 0 Silicon 192 2 0 0 end buckingham del 175 0000 172 133 5381 22 0 0 end 4 903 4 903 5 393 90 90 Oxygen 0 415000 0 Oxygen 0 857000 0 Oxygen 0 728000 0 Oxygen 0 143000 0 Oxygen 0 272000 0 Oxygen 0 585000 0 Silicon 0 465000 0 Silicon 0 535000 0 Silicon 0 0 end 0 16 ds 0 28 0855 2 4 1388 7730 2 76000 18003 7572 4 87318 0 0 0 0 120 4 4 272000 0 120000 585000 0 453300 143000 0 453300 728000 0 880000 415000 0 546700 857000 0 213300 0 535000 0 333300 465000 0 666700 57 0 Appendix Utility Programs Seven utility programs are included in the Moldy package mostly for manip
18. B I integer parameters Parameters such as dump level take a numeric value which should be an integer timestep related parameters Several parameters govern when some calculation begins and ends and how frequently it is performed in between These are known as begin end and interval parameters but are really a special case of integer parameters For example begin average dump interval and scale end The calculation begins on the timestep specified on the begin parameter occurs every interval timesteps thereafter and ends after the timestep specified by the end parameter Setting the interval parameter to zero is the usual method of turning that calculation off The begin and end parameters behave in a special fashion when the simulation is continued from a restart file they are interpreted relative to the current timestep Notice especially that nsteps the number of timesteps is treated in this way A complete list of the parameters their meanings and default values appears in table B I 3 2 Setting up the System 3 2 1 System Specification The information which describes to Moldy the system to be simulated and the interaction potentials is contained in a file known as the system specification file This may be presented to Moldy in either of two ways If the control file parameter sys spec file is null or absent it is assumed to be appended to the end of the control file Otherwise it is read from the file whose name is th
19. Cambridge MA 02139 USA Goodness knows who will own the Unix registered trademark by the time you read this Chapter 2 Algorithms and Equations This chapter describes the implementation of the molecular dynamics technique as used in Moldy It is not intended as an introduction to MD simulations and does assume some familiarity with the basic concepts of microscopic models and simulations thereof The book by Allen and Tildesley Z is a very good introductory text covering both the theory and the practicalities and is highly recommended It also contains comprehensive references to the scientific literature of microscopic computer simulation Moldy is designed to simulate a common class of models of atomic or molecular systems The assumptions are that the system is an assembly of rigid molecules atoms or ions that the forces of interaction are derived from continuous potential functions acting between usually atomic sites on each molecule that the dynamics are governed by the classical Newton Euler equations of motion A major aim of Moldy has been to allow the most general of models within that class and to impose as few restrictions as possible In particular arbitrary mixtures of different molecules are allowed and several popular forms of potential functions are catered for 2 4 The Equations of Motion If we denote the force exerted by atom 04 of molecule i on atom B of molecule j by fia sl then the total force acting on m
20. It is therefore handled using the parallel link cell method exactly as for the site forces calculation by functions rdf inner and rdf accum However the contributions from different processors do not need to be summed every time inner is called but may be left to accumulate separately The data is globally summed by function par isum only when rd out is called to calculate and print the RDFs or just before a restart file is written 55 Appendix Example System Specifications Argon LJ Argon about as simple as you can get Parameters from Allen and Tildesley Table 1 1 Argon 108 1 0 0 0 39 948 0 Ar end Lennard Jones 113 984 3 41 2 TIPS2 Water This is the four site water model of Jorgensen Z5 Only the oxygen site interacts via the Lennard Jones potential and the charge site M is displaced 0 15A from the Oxygen Modified TIPS2 water Water 64 1 0 0 0 16 0 0 2 0 7569503 0 0 5858822 1 0 535 H 2 0 7569503 0 0 5858822 3 0 0 0 15 0 1 07 M end lennard jones 1 1 0 51799 3 2407 Aqueous MgCl Solution This is a three component system consisting of MCY water B3 magnesium and chloride ions The Mg potential was fitted to the SCF calculations of Dietz and Heinzinger II7 and the to the calculations of Kistenmacher Popkie and Clementi 29 Note that the potential parameters are expressed in kcal mol and the control file must set the parameter time unit 4
21. Xi Yi Zi species name Xn Yn Zn 40 qni dm qm end Here a b c and o D y are the crystal unit cell parameters and are the number of unit cells in each direction which comprise the MD cell The next lines describe the molecules of the basis which will be replicated to form the full configuration Molecules may appear in any order but of course the total number of each multiplied by the number of unit cells 1 n n must agree with that given in the system specification file 23 Each molecule is identified by its name as given in the system specification file 7 are fractional ordinates between 0 and 1 giving the location of the molecular centres of mass in the crystal unit cell The orientation is given by the four quaternions 40 41 42 43 Which specify a rotation about the centre of mass relative to the orientation of the prototype molecule in the system specification file Notice the slight inconsistency with the positions which are of the centres of mass not the zeroes of co ordinates in the system specification file This may be fixed in future releases Quaternions need only be included for polyatomic species that is molecules 1 and above and omitted for the monatomic species i After the molecular positions and orientations have been set up their velocities and angular velocities if appropriate are initialized Their values are sampled at random from the Maxwell Boltzmann distribution for the tem
22. a multiprocessor machine 4 2 Portability A major goal in writing Moldy was that it be as portable as possible between different computers and operating systems It is written in the Kernighan and Ritchie Z7 compatible subset of ANSI C and assumes the library calls and header files defined for a hosted implementation of the standard It should therefore be possible to compile and run Moldy on any computer which has a good C compiler The configure script immensely eases portability on unix like systems by detecting the capabilities of the compiler and the installed header files and library functions The details are given below for completeness or in case of a configure failure or for a system where configure does not work Though hosted ANSI standard C environments are now commonly available it may still be necessary to compile Moldy using a pre ANSI compiler Replacement ANSI library functions are supplied in ansi c for the VMS and unix operating systems For ease of portability all other system dependent functions are in the module auxil c and all prepro cessor conditionals are in the header file defs h If the target machine has ANSI conformant C libraries all that must be done is to define the preprocessor symbol ANSI LIBS either in defs h or by using a compiler option e g DANSI LIBS 4 2 1 System Dependencies In this section details of system dependent functions are described for the major operating systems Replacement ANSI heade
23. and shape in response to the imbalance between the internal and externally applied pressure For an exposition of the method the reader should refer to the papers of Parrinello and Rahman B9 Wentzcovitch 56 Cleveland IT and to Nos and Klein s extension to rigid molecular systems B7 The equation of motion for the reduced centre of mass co ordinates 8 h R 18 MS MiG GS 2 32 replacing the straightforward Newton equation 2 3 h denotes the 3 x 3 MD cell matrix whose columns are the MD cell vectors a b and F is the centre of mass force and h h Equation 4 the Euler equation governs the angular motion exactly as in the constant volume case In the Parrinello Rahman method the dynamics of the unit cell matrix h are governed by the equation Wh x p o 2 33 where W is a fictitious mass parameter O and is the external pressure The instantaneous internal stress 7 is given by 1 his S T 2 34 B iSi 1 00 2 34 with the short ranged and electrostatic components given by equations 2 15 and 2 24 respectively Wentzcovitch 56 and Cleveland II T showed that the Parrinello Rahman equations suffer from at least one signifi cant flaw namely that they are not modularly invariant The dynamics depends on the choice of MD cell which can give rise to unphysical symmetry breaking effects They offer an improved but more complicated expression for t
24. be disabled by setting specific values for control alpha The exact behaviour is If control alpha gt 1077 is set to that value and used in the Ewald sum e If control alpha 0 then the value of amp is chosen automatically as above e If control alpha lt 0 the reciprocal space part of the Ewald sum is not invoked and the computation of coulom bic forces in the real space sum is turned off This happens automatically if all charges in the system specification file are zero e If 0 lt control alpha lt 1077 then the reciprocal space part of the Ewald sum is not called but coulombic interactions are included in the real space part and truncated at re The other relevant parameter is the switch surface dipole which includes the dipole surface energy term of De Leeuw Perram and Smith BZ See the note in section 2 4 for an explanation of why this term should not be used for an ionic as opposed to dipolar system The two adjustable parameters which control the link cell force calculation see section 5 are subcell and strict cutoff The former specifies the length in of the side of a link cell and determines the number of cells the MD cell is divided into In fact the MD cell is divided into a whole number of subcells whose side in each of the three directions is nearest to the value of subcell The default of zero though is special and sets subcell to one fifth of the cutoff radius In general the smaller the link cell the
25. c Primary initialization function start up plus subsidiary functions allocate dynamics to create dy namic variable arrays and initialise sysdef to compute molecular and whole system properties Also con tains functions to create initial configuration for system values c Function init averages allocates and sets up pointers to averages database Other functions store and retrieve instantaneous values of thermodynamic quantities and compute the averages thereof xdr c Functions to read and write Moldy s own struct types using the External Data Representation XDR library 5 3 Flow of Control 5 3 1 Input and Initialization The initialization and input file reading stage of Moldy is rather more complicated than is usual in a molecular dynamics simulation program because Adjust timestep Find Md 1 related control rom sys spec file parameters 0 rom backup replace conv potentials system spec check sysdef abort if failed Skew start interpolate_ timeste derivatives changed Figure 5 2 b Block diagram of the initialization function start up and a list of the functions called Continued from Figure The paths beginning at A B and are for a new run a restart from a save file and a restart from a backup file respectively of the use of dynamic data structures The size of the arrays required can only be dete
26. described hereafter 2 10 1 Implementation The framework is in many respects exactly like an ordinary molecule It should be defined to exactly fill the MD cell so that the periodic repeat of the cell generates the periodicity of its crystal structure The distinctive features of a framework molecule are The framework is fixed in space and is not allowed to rotate Any rotation would of course destroy the 3D structure For most purposes it is convenient to regard the framework as being at rest hence translational motion is forbidden also No interactions between sites on a framework molecule and on itself or any of its periodic images are evalu ated That is framework framework interactions both potentials and forces are systematically excluded from the real space and reciprocal space parts of the Ewald sum including the point and intra molecular self terms of equation 2 19 The exact modifications to equations 19 20 etc are left as an exercise for the reader thereal space force calculation sites are treated as being independent atoms rather than belonging to a molecule In particular the cutoff is applied on the basis of the framework site to fluid molecule distance By virtue of the all image convention all the requisite molecule framework interactions are correctly included When assigning sites to sub cells each site is therefore placed in the sub cell which contains its co ordinate By contrast sites bel
27. deviations 1 8214 71 893 56 441 19942 23 90 1 43 4 0 00 0 00 0 00 1 32 03 750 51 0 013 0 17 173 0 1 0 0 0 00 000 0 00 750 19 552 0 00 0 00 0 00 51 552 49 Figure 3 1 Sample Moldy output from a simulation of a two component mixture The first component is a polyatomic molecule and the second is atomic There are three frames for the instantaneous values the rolling averages and their associated standard deviations Within a frame each row has the following meaning for translational and rotational kinetic energies and temperatures it is the per species value for the potential energy it is the direct and reciprocal space components and the MD cell matrix h and the stress are laid out as 3 x 3 matrices However in the traditional implementation of scaling velocities are multiplied by a factor of J desired KE instantaneous KE Thus the scaling factor is wrong by the ratio of the instantaneous to average KE s which means that the temperature can not be set more accurately than the relative size of the fluctuations in the KE The option selected by bit 2 goes some way towards the ideal scaling factor by using the rolling average KE instead of the instantaneous value The fluctuations in this short term average should be much lower than in the instantaneous value allowing more accurate temperature control However it will almost always be easier to use a true thermostat to achieve this goal 3 5 Output At the beginning of each run Mol
28. double see section 5 1 1 code for par rsum must handle either case which may be tested using the sizeof operator For example the preprocessor macro define M REAL sizeof real sizeof double MPI DOUBLE MPI FLOAT is used to determine which constant to pass to the MPI sum function 38 Chapter 5 Program Structure The source of Moldy consists of 31 different C source files amounting to 167 functions and 9000 lines of code A complete and detailed description would be a compendious volume of questionable value Instead much of the detailed documentation is contained in the source code in the form of comments This chapter concentrates on describing the organization of the calculation and the data structures used throughout the code describes some of the more complicated and less obvious algorithms and provides call graphs and brief function descriptions to act as a map of the program structure Moldy is written in modular fashion in adherence to the principles of structured programming to as great a degree as practical This does not mean merely the avoidance of gotos but instead the organization of the program into modules and functions which are independent of each other and of global environmental variables Functions are kept to a limited size and as far as practical serve a single well defined purpose This ensures that the internal workings of a function are unaffected by changes elsewhere in the program and do not have any
29. e e ed ee Ae a ROS ed Rs UN SIR e E Ge 8 D Potentials and Short Range Forces 2 0 ol 2 10 arged FramewoOrks bq aan ao a mue dike x eX ge BR g Moldy 3 pup the System ee PAL 5 2 stem pecificatiOD B 2 elnitialConieurationy 2 2 u oso e oom RR o X RO E XN ee D Restart ng from 5 D P dic Ds s eo ox Ro Fus X Sox Xo 3 4 etting the lemperature s es e sa ma reaa deee e e e e a ss se D sni D Dutput units 25222 55 E CE Gu US S CR VOCE Y a E RU GN O 4 HNN B ke ER doe amp BA we SE 27 B8 Constant Stress or Pressure Simulation 29 VV 29 BIO Framework Simulations 2222s s sss 30 D Messages and ETTOT9S x09 po gt vox sce cr wow Xo We RO Wo VOR M ADAC UR UE E RU 31 4 ompiling and Modifying Moldy 32 ompilation eR as 32 T T Parallel Version Message Passing 34 4 ared Memory Parallel 35 4 Qua ML
30. for the centre of mass acceleration depends indirectly on the velocities through 7 as well as explicitly The iteration on steps iii v of equations 2 215 therefore essential to integrate these equations consistently There are two criteria to be considered when choosing the value of Q The coupling to the heat bath introduces non physical oscillations of period to 214 Q 2gkgT which may be easily detected in the total energy B6 Firstly it must be arranged that there are sufficiently many oscillations during a simulation run that the computed thermodynamic values represent averages over many periods This ensures that the configurations used to calculate the averages sample the whole of the phase space generated by the fluctuations in G to represent the canonical ensemble Secondly Q should be chosen so that fo is large compared to the characteristic decay time of dynamical correlation functions This is to ensure that the fictitious dynamics of the heat bath are decoupled from the real molecular dynamics and is particularly important in a simulation of dynamic amp such as velocity correlation functions Since the first criteria favours a small Q and the second a large one it may be necessary to increase the total simulation time in order to satisfy them both 2 6 3 Gaussian Thermostat An alternative approach is known as constrained dynamics whereby the equations of motion are modified to generate trajectories whi
31. influence on any other part of it except through the defined interface In Moldy functions are grouped into different files according to a rough classification of their purpose The other primary consideration in the design of a large computer program is the provision and organization of storage for the data Structured programming advocates that definitions of data objects be restricted as closely as possible to the code that uses them and that the code be organized in a modular fashion to encourage data locality This minimizes the risk of a programming error in one part of the code modifying a variable used in a completely different part and producing difficult to locate side effects With two exceptions f global data is avoided in Moldy and all arrays are passed as function arguments where necessary Heavy use is made of C structures or structs to further group related data so that it may be manipulated in toto 5 1 Data Structures 5 1 1 Types and Typedefs A number of derived types are defined in the header file defs h for reasons of economy of expression portability and ease of customization Most of these derived types are named using a suffix _mt to mark them as such These are real This is the standard precision floating point type used for all the dynamic variables and most other internal variables It is set to double by default Changing this to float will give a single precision version of Moldy with a consequent memory s
32. is a restart then those control parameters whose value is interpreted relative to the current timestep see section B I must have its value added to theirs This is only done if the parameter was explicitly specified in the new control file otherwise its saved value is untouched Function re re sysdef reads the system specification variables and arrays system species potpar and site info from the restart file These were stored in their final form after being set up by the previous run so a call to initialise sysdef is unnecessary Alternatively a new system specification file may be used in which case a similar sequence to that for a new run is followed Memory for the dynamic variable arrays and the RDF and averages databases must be allocated before their saved values are restored by read restart The simplest startup mode is that from a backup file Once it has been determined that a backup file exists setup proceeds much as in the restart case except that the parameters and system specification are restored unchanged Thus the run continues exactly from the point the backup file was written A restart file is still opened but only to obtain the name of a backup file to search for The final act of start is to call banner page to print out a nicely formatted page of information to record the simulation details in the output file Control is then transferred back to the main program where everything is now ready for the mole
33. molecule The computation of the real space potential and forces is distributed over processors on the basis of link cells For the reciprocal space part of the Ewald sum the k vectors are distributed among processors This is an extremely efficient way of implementing parallelism since the forces etc must be summed over processors only once per timestep thus minimizing interprocessor communication costs It is therefore possible to get considerable speedup for a small number of workstations coupled by a fast network such as an Ethernet The biggest disadvantage of the replicated data strategy is that every processor must maintain a copy of all of the data and therefore that the memory requirement per processor increases with the size of the system In many cases this is not a severe problem as MD memory requirements are not large compared with memory sizes of modern computers However the poor scaling will eventually limit the number of processors which may be used On a shared memory multiprocessor the alternative parallel version in section may provide a solution if it can be ported to that machine The memory limitation will be most acute when the goal is to simulate an extremely large system on a massively parallel distributed memory computer where it is desirable to scale the system size with the number of processors In such architectures the available memory per processor is usually a constant independent of the number of processors But the me
34. more accurately the cutoff radius is implemented but too many of them reduces the efficiency of the program In the default cutoff mode st rict cutoff false the list of neighbour cells is constructed to include all cells whose centre centre distance is less than the cutoff This means that some molecule pairs separated by more than the cutoff will be included and some by less will be omitted Setting strict cutoff to true generates a larger cell neighbour list which is guaranteed to include all appropriate molecule pairs Furthermore molecules separated by more than the cutoff are excluded from the force calculation by setting their separation to a large number 100 times the cutoff at which distance it is assumed the potential is very close to zero This is therefore the mode of choice for liquid simulations where any artificial anisotropy is undesirable See section 2 5 for a full explanation It is worth noting that it is unnecessary to recompile the program or change anything else when the cutoffs are modified Unlike most MD codes Moldy employs dynamic array allocation and automatically sets up arrays of the correct size and no more for any given ke 3 10 Framework Simulations There has recently been much interest in simulations of systems of molecules interacting with some rigid framework such as zeolites clays and other surfaces Moldy has the capability to include such a framework in a simulation by defining it as a special kind of mole
35. potential params into neighbour list arrays Calculate squared pair site distances in neighbour list force inner test for close site approaches Print warning kernel Evaluate potential amp force on neighbour list Calculate stress amp forces on sites in neighbour list Distribute neighbour list forces onto force array M ML Figure 5 5 The Link Cell short range force calculation velocities have converged f calculate the new angular and rotational thermostat accelerations The Euler equations 2 4 and equation 2 10 are used by function euler to evaluate the second time derivatives of the molecular quaternions Iterate until the corresponding angular velocities have converged g call step 2 to apply the final step iv of equations 2 13 h write requested co ordinates etc to the dump file The loops of steps f iterate until the molecular centre of mass velocities and quaternion time derivatives have adequately converged in practice about 3 or 4 iterations No additional iteration is applied to converge the extended system dynamic variables since in any practical system these should be very slowly varying compared to molecular motion and any error introduced will be very small 50 1 Calculate self and sheet energy terms Precompute list of h k l for k vectors within cutoff
36. r i j double r rnt d J r i j 0 0 since the underlying pointer arrays contain all the necessary information to access indivual elements Function arralloc is implemented so as to lay out the individual pointer and data arrays as part of a single block of memory which is allo cated by a single call to 11 It then sets up the values of the pointers to emulate the requested array The memory can be freed with a single call to t ree The arralloc mechanism is unnecessary for 2D arrays where the innermost dimensions are fixed such as for an array of position co ordinates which has dimensions n 3 In such cases one of the xalloc macros is used to allocate an array with fixed innermost dimensions and whose type is a pointer to a fixed size array double x 3 rather than a pointer to a pointer double x 3At the time Moldy was being planned in 1988 C was the only widely available language offering dynamic memory allocation Fortran 90 which also offers dynamically declared arrays was standardized in 1991 and compilers only became common in the mid nineties 40 0 1 2 3 2 n l 0 1 221 2 L m 1 r r i r i j pointer array of m pointers m arrays of n data objects Figure 5 1 Storage layout of a 2 dimensional pointer b
37. t 8 F px r 81 4 8t i 1 n 2 13 iv O t 8r 3 Dx r 4 8r S x r r v Replace x P with x and goto iii Iterate to convergence The predictor corrector cycle of steps iii to v is iterated until the predicted and corrected velocities have converged to a relative precision of better than 1 part in 1077 which in practice takes 2 or 3 cycles This iteration is not as inefficient as it might at first appear as it does not include the expensive part of the calculation the recalculation of the site forces Only the angular accelerations and quaternion second derivatives must be evaluated at each iteration using equations 2 4 and 2 10 and this operation is relatively cheap 2 3 Potentials and Short Range Forces The forces determining the dynamics of the system are derived from the potential function denoted by jp The indices i and j run over all molecules in the system and o and p over sites on the respective molecule The total potential energy of the system is U tia 2 14 i where fiaj jB F io 38 is the force acting on site B of molecule j from site of molecule i The total force and torque acting on any particular molecule are calculated using equations e and 4 Since and therefore are short ranged forces i e they decay faster than r one can define a cutoff radius re Interactions between sites whose separation r is greater
38. the calculation of dynamic properties such as time correlation functions and additional static averages not normally calculated by Moldy During a run dump files are produced which contain a record of the simulation dynamic variables positions quaternions etc at varying degrees of detail Any property of interest dynamic or static may then be evaluated using the data in the dump A dump consists of a sequence of files since the amount of data generated in a run can be very large indeed and it is usually more convenient to manipulate a series of smaller files rather than one large and unwieldy one Moldy takes considerable pains to ensure that a contiguous sequence of dump files is maintained and also ensures that dumps from different runs are not accidentally intermixed There is no requirement that a dump file be produced by a single run of Moldy which extends an existing file or starts a new one as appropriate simulation may stop and be restarted many times without disturbing the dump file sequence The sequence should in most cases even survive a system crash and a restart from a backup file see section B 3 1 previous versions of Moldy the calculation of the interatomic distances was done on the basis of the minimum image convention Conse quently the calculated value of tailed off for re gt L 2 This restriction is now lifted 27 Each dump file in a sequence is a binary file consisting of a dump header wh
39. the potential and the forces on the neighbour list sites The actual potential and scalar force evaluation is delegated to function kernel for ease of comprehension and modification This takes as its arguments an array of squared pair distances and corresponding arrays of potential parameters and site charges and returns the potential energy and an array containing the scalar part of the force rij rij The structure of kernel is designed to evaluate all the different kinds of potential functions see section 2 3 3 as efficiently as possible It contains separate versions of a vector loop for each distinct potential type both with and without electrostatic charges For the charged case the complementary error function is approximated to an accuracy of approximately 1 part in 1076 by the formula given in Abramowitz and Stegun IT section 7 1 26 Pair distances for the calculation of radial distribution functions are evaluated by rdf_inner whose structure resembles a simplified force inner It does not call kernel or compute forces but instead calls accum which bins the computed pair distances The separate function allows the use of a larger cutoff than for the force and potential evaluation The reciprocal space part of the Ewald sum is evaluated in function ewald shown in figure 9 This is a fairly straightforward implementation of the k space terms of equations 2 19 2 20 and 2 24 and differs from common imple menta
40. to move it to any other machine with a good compiler To be of real use a simulation program must run efficiently on modern high speed computers which are increasingly of vector or parallel architectures Moldy is written so as to be highly vectorizable and has been tested on a range of vector machines from manufacturers including Cray Convex Stardent and Alliant On the cray XMP 48 its performance can exceed 100 MFlop sec on a suitably large system Moldy is also able to run on a parallel computer of either shared or distributed memory architectures and has been tested on multiprocessors from Stardent Convex Cray Research SGI SP1 and Cray Research T3D massively parallel machines 1 1 Terms and Conditions Permission to compile and run Moldy is granted to any individual or organization without restriction This program is distributed in the hope that it will be useful but WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE See the GNU General Public License for more details Moldy is free software you may redistribute it and or modify it under the terms of the GNU General Public License as published by the Free Software Foundation either version 2 or at your option any later version The terms of the GNU General Public License are described in the file COPYING which is included with the source distribution or from the Free Software Foundation 675 Mass Ave
41. was often necessary to write programs in machine specific languages HI use proprietary and obscure library calls IJ and be designed around a specific communications topology Nevertheless molecular dynamics simulation was one of the earliest applications of parallel computers and many simulations have been undertaken despite the diffi culties 21 ITO The emergence and widespread adoption of the single program multiple data SPMD programming model and the standardization of parallel communications libraries in the 1990s has much improved the situation In this model a parallel program is treated as a set of copies of a single program executing asynchronously and independently on different processors and which communicate using library calls The program is written in a standard language and arranged to divide the calculation among the processors in some way so that each handles a subset of the operations The model is completed by the addition of a library of communications routines which allow the separate parts of the calculation to be assembled to give the full desired result A significant advance was the availability of parallel libraries such as TCGMSG ZI 8 BSP B4 and the new standard MPI T7 Together these developments allow a program to be 52 both parallel and sufficiently independent of architectural details to be portable to a wide range of parallel environments These may include shared memory and distribut
42. with very low communications latency run this version effectively Practical tests show that the communications overhead completely negate any parallel gain on systems of networked workstations and most multiprocessors However a significant speedup is obtained on a Cray T3D which is exactly the case where this version is needed 4 11 2 Shared Memory Parallel Version An alternative parallel version is available for shared memory multiprocessors with parallelizing compilers This relies on the compiler handling the multi threading synchronization and allocation of local memory stacks for inner function calls It requires compiler specific directives to be inserted in the code and is therefore less portable than the distributed memory version of the previous section Note that that version works on this class of machines too under the message passing interface Nevertheless it works and has been run on Stardent Convex and Cray computers It consists of replacements for files force c and ewald c called force parallel c and ewald parallel c Then the program should be compiled with the preprocessor macro PARALLEL defined not SPMD The distributed memory parallel version section 1 is generally recommended over this one However because the parallel sections reference a single global copy of most of the arrays the shared memory version uses much less memory This version may therefore be of use if memory limits the size of the system on
43. 000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000481 035710 0 183334 0 86 0 613992 1 024402 1 046396 0 964906 0 830174 0 660035 0 693341 0 615902 0 593192 0 510595 0 530697 030 0 535959 0 57 0 523221 0 466219 0 496028 0 438487 0 456500 0 410547 0 443861 0 457956 0 446822 0 452202 0 419768 0 439333 0 465509 0 486887 0 461970 0 475745 0 478883 0 480854 0 509090 0 533728 0 552747 0 552555 0 575402 0 547278 0 544836 0 493597 0 488168 0 520727 0 508073 0 479948 0 501159 0 484000 0 485378 0 489160 0 464448 0 466791 0 476508 0 446576 0 470948 0 474468 0 449340 0 462169 0 501220 0 519107 0 513338 0 510192 0 499766 0 525963 0 504663 0 517673 0 498359 0 512156 0 507061 0 466390 0 464342 0 445886 0 417555 0 407778 0 387220 0 374041 0 RDF 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 0 000000 26 976688 0 000000 0 000000 0 000000 0 000000 0 000000 0 016214 0 061257 0 304082 0 647342 0 847404 0 757188 0 0 478273 0 462682 0 449614 0 450424 0 518998 0 572242 0 689704 0 914269 1 184674 1 441772 1 570390 1 609068 1 600392 1 430457 1 322722 1 183606 1 103701 1 061788 0 980018 0 960570 0 924390 0 908883 0 877591 0 857668 0 890761 0 852463 0 815447 0 824963 0 841255 0 890416 0 929030 0 960589 0 984145 1 020650 1 028199 1 047496 1 064600 1 099812 1 095715 1 073793 1 078131 1 049212 1 052160 1 1 020737 1 010782 0 979748 0 983158 0 988946 0 967620 0 955655 0 944384 0 952145 0 948509 0 946692 0 9
44. 2 20 if used performs the calculations under different periodic boundary conditions It was suggested by De Leeuw Perram and Smith B2 in order to accurately model dipolar systems and is necessary in any calculation of a dielectric constant The distinction arises from considerations of how the imaginary set of infinite replicas is constructed from a single copy of the MD box Z pp 156 159 Consider a near spherical cluster of MD cells The infinite result for any property is the limit of its cluster value as the size of the cluster tends to infinity However this value is non unique and depends on the dielectric constant of the physical medium surrounding the cluster If this medium is conductive s the dipole moment of the cluster is neutralized by image charges whereas in a vacuum 1 it remains It is trivial to show that in that case the dipole moment per unit volume or per MD cell does not decrease with the size of the cluster The final term in equation 2 is just the dipole energy and ought to be used in any calculation of the dielectric constant of a dipolar molecular system It is switched on by Moldy s control parameter surface dipole Note that as it represents the dipole at the surface of the cluster the system is no longer truly periodic Conversely it must not be used if the simulated system contains mobile ions Consider an ion crossing a periodic boundary and jumping from one side of the MD cell to a
45. 35 17 stem Dependencie eA 35 F Z Optimization and 36 T2 Optimization for Vector Architecture 2 222 les 36 4 uunc 37 Adding a New Potential 37 4 3 Porting the Parallel Version 37 5 Prog 39 Data Structured 2 222 22 0 E 39 pes and Typedefg ko pomo woe wee E ACE NDA ROBUR EUR AE SUR REOR X 39 Memory Management 40 ee EE es ee 41 E2 HlesandFunctiond 42 BZI SoureFled iieegg heec twee m MG S RA d S E E pr RU e ed UR d 42 E3 Flow 44 5 3 nput and Initialization 44 5 3 Main timestep ss ee HA o ER ER 46 6 3 3 Site Forces Calculation 51 b 4 Parallelizatiotl 222255222 ede RA AeA Re weg eee ERR ay 52 5 4 Parallel StraleBy ee RE EA ee SR ee 53 aa ee ee wee od at bee eee ote 53 5 4 put Outp d Duck oe BS ep o eR RE c d d CUR 54 54 22 xample System Specification 56 A I Es Re ety ae N
46. 60097 0 959299 0 0 969219 0 972704 0 998504 1 027791 1 041576 1 037637 1 039961 1 016804 1 004726 1 026805 1 030903 1 006268 0 972421 0 948140 0 908959 0 877089 0 849855 0 817964 0 776986 0 721485 Figure 3 2 Example output of radial distribution functions After the header line consisting of underscores there is an indication of the bin width b that is the distance between points at which the RDF is tabulated Then for each site type pair there is a line listing which pair e g 0 0 RDF followed by nbins values of ggg i 1 2 b 3 6 Radial Distribution Functions Radial distribution functions are calculated by binning site pair distances periodically throughout the simulation see section 2 8 As this process is expensive in computer time the binning subroutine is invoked only every few timesteps as set by the control parameter rdf interval 20 by default Since the pair distances only change a little on each timestep very little statistical information is lost Collection of binning data may also be turned off during an equilibration period specify when binning is to start by means of the parameter begin rdf The parameters rdf limit and nbins control the details of binning giving respectively the largest distance counted and the number of bins that interval is divided into The calculation of the interatomic distances is done separately from that used in the evaluation of the forces using the same link cell scheme This ensures that al
47. C convention This constant should be used as the label to a new case in the switch statement of kernel and this is where the code to evaluate the potential goes The existing cases may be used as a model especially for the evaluation of the electrostatic term erfc ar r which is evaluated by the polynomial expansion of Abramowitz and Stegun M section 7 1 26 There are currently two versions of each loop the second omitting this term for efficiency when all the electric charges are zero which case is flagged by a negative value of e Finally the distant potential correction for the new potential should be added as a new case to function dist pot The code should evaluate T rU dr fe for the potential U r 4 3 2 Porting the Parallel Version It should be relatively straightforward to port the distributed memory parallel version to a new message passing library Section 4 1 describes the parallel implementation All of the interface code is contained in parallel c and it will only 6A mechanism is provided to insert appropriate directives using the C preprocessor The text VECTORIZE has been placed before each loop which ought to be vectorized and the file defs h contains machine conditional 4 defines to replace it with the appropriate directive Currently directives for the CRAY Stellar and Convex compilers are included and null text is substituted for other machines Notice that in each case the substituted
48. Comp Phys 13 1975 430 432 67 46 D C Rapaport Large scale molecular dynamics simulation using vector and parallel computers Comp Phys 9 1988 no 1 1 53 47 Refson Molecular dynamics simulation of solid n butane Physica 131B 1985 256 266 48 K Refson and G S Pawley Molecular dynamics studies of the condensed phases of n butane and their transitions I Techniques and model results Mol Phys 61 1987 no 3 669 692 49 P M Rodger On the accuracy of some common molecular dynamics algorithms Mol Simul 3 1989 263 269 50 W Smith Molecular dynamics on hypercube parallel computers Comp Phys Commun 62 1991 229 248 51 A replicated data molecular dynamics strategy for the parallel Ewald sum Comp Phys Commun 67 1992 392 406 52 W Smith and D Fincham The program MDMPOL SERC Daresbury Laboratory 1982 CCP5 Program Library 53 R Sonnenschein An improved algorithm for molecular dynamics simulation of rigid molecules J Comp Phys 59 1985 no 2 347 350 54 P Du Val Homographies quaternions and rotations Oxford Mathematical Monograph 1964 55 L Verlet Computer experiments on classical fluids I thermodynamical properties of lennard jones molecules Phys Rev 165 1967 201 214 56 Renata M Wentzcovitch Invariant molecular dynamics approach to structural phase transitions Phys Rev 44 1991 no 5 2358 2361 68
49. DCL command file called moldy com and is unpacked by the command moldy MS Windows 3 Windows 95 NT The shareware program winzip may be used to unpack the compressed tar archive moldy tar gz There is also a precompiled binary version for Windows 95 98 NT which is available as moldy zip Winzip may also be used to uncompress this distribution Two versions of each executable are provided in di rectories cygwin and mingw One of these sets should be installed in a suitable directory which should be added to the user path The Cygwin compiled binaries rely on the shared library cygwin1 dll and incorporate support see section If Cygwin is not already installed the DLL included must also be installed in a suitable directory on the system The Mingw32 compiled executables do not require a dll and should run on any Windows machine but do not support the creation or reading of XDR format restart files or dumps MS DOS files must be unpacked from the tar archive on another host and transferred to the PC by disk or ftp 4 1 Compilation The source code of Moldy consists of 23 files of C programs which have suffix c and 8 header files suffix h To build Moldy all of the c files must be compiled and linked together to form an executable The method of doing this depends on the operating system of the target computer unix This version of Moldy uses the GNU autoconf system to configure for compilation In mo
50. It remains to choose the initial velocities of the molecules to complete the specification of the initial configuration The recipe is the same irrespective of whether the molecules are started on a lattice or from a skew start The initial centre of mass velocities are chosen from the Maxwell Boltzmann distribution at the temperature specified for the simulation Z pp 170 That is the velocities are chosen from a set of random numbers with a Gaussian distribution and normalized so that the probability density p v of the xyz component of the velocity vj of molecule k is 1 2 2 my 2 47 P vit zziz Expl This is easily accomplished given a random number generator which samples from a Gaussian distribution with unit variance Given a random number Rx each component of velocity is set to Vik 4 Rik 2 48 Each component of the angular velocity expressed in the principal frame has a probability distribution E W Ij l eik 2 4 pon exp 2 49 which is ensured by assigning a random velocity gie e 2 50 Since the dynamical variables used by Moldy for the angular co ordinates are in fact the quaternions and their derivatives we must set the quaternion derivatives and accelerations to the corresponding values Using equations 9 and LIQ we have q 0 2 2 51 and 4 2 52 Finally we note that if a molecule has less than three degrees of f
51. KAL format default give the output in Brookhaven Protein Data Bank PDB format give the output in XYZ format suitable for Xmol or RasMol name of optional output file Defaults to standard output B 9 Mdbond and bdist Mdbond is a utility for calculating the distances between the centres of mass of pairs of molecules as well as angles between triplets of molecules within minimum and maximum values specified using options and a for bonds and angles respectively The utility calculates values for ALL pairs triplets within the cutoffs including those involving mirror atoms so that maximum cutoffs greater than the simulation box lengths are possible Bond lengths are output in columnar format in increasing order with the number and species type of each pair of molecules given in the first two columns The angle data are output after the bond data also in columnar format Each row contains the molecule numbers and species types for triplets i j k followed by the angle formed about molecule i and the distances i j and i k The command line options are mdbond s system specification r restart file dump file format t dump range b bond limits a angle limits g species range c o output file where the arguments have the following meanings 5 0 read a system specification file read a restart file Only one of s or r may be given read configurational data from a dump file given as a prototy
52. M Meyer and V Pontikis eds vol E205 Kluwer Dordrecht 1991 NATO ASI pp 21 42 37 S Nos and M L Klein Constant pressure molecular dynamics for molecular systems Mol Phys 50 5 1983 1055 1076 38 Rafael Rodriguez Pappalardo and Enrique Sanchez Marcos Recovering the concept of the hydrated ion for mod elling ionic solutions A monte carlo study of 21 in water J Phys Chem 97 1993 4500 4504 39 M Parrinello and A Rahman Polymorphic transitions in single crystals A new molecular dynamics method J App Phys 52 1981 no 12 7182 7190 40 S Pawley Molecular dynamics simulation of the plastic phase a model for 8 6 Mol Phys 43 1981 no 6 1321 1330 41 Computer simulation of the plastic to crystalline phase transition in SFg Phys Rev Lett 48 1982 410 413 42 G S Pawley and M T Dove Quaternion based reorientation conditions for molecular dynamics analyses Mol Phys 55 1985 no 5 1147 1157 43 J Powles W A B Evans E McGrath E Gubbins and S Murad computer simulation for a simple model of liquid hydrogen chloride Mol Phys 38 1979 893 908 44 William H Press Saul A Teukolsky William T Vetterling and Brian P Flannery Numerical recipes in C the art of scientific computing 2nd ed Cambridge University Press 1992 45 B Quentrec and C Brot New methods for searching for neighbours in molecular dynamics computations J
53. Moldy User s Manual Revision 2 25 2 6 for release 2 16 Keith Refson Department of Earth Sciences Parks Road Oxford OX1 3PR Keith Refson earth ox ac uk June 5 2003 Copyright 1988 1993 1996 Keith Refson Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies Permission is granted to copy and distribute modified versions of this manual under the conditions for verbatim copying provided also that the entire resulting derived work is distributed under the terms of a permission notice identical to this one Permission is granted to copy and distribute translations of this manual into another language under the above conditions for modified versions Acknowledgements The initial implementation of the Nos Hoover and Gaussian thermostats was by Vladimir Murashov at Dalhousie University Rafael R Pappalardo contributed the modifications to mdshak to write XYZ format molecular graphics files Craig Fisher contributed the msd mdavpos bdist mdbond and ransub programs Thanks are due to Craig Fisher and Jos Manuel Martinez Fernandez for diligent proofreading of several versions of the manual and many helpful suggestions on how to improve it Any mistakes that remain are my own Contents Introduction 2 Algorithms and Equations D e Equations of Motion Rm ors 2 Aleorithms e
54. SG BSP or SHMEM defined Tt is essential that the particle co ordinates are identical on every processor and remain so throughout the simulation run Since the equations of motion are chaotic the tiniest numerical difference in the forces between processors will cause the trajectories to diverge It is not difficult to arrange that a global sum returns results identical to the last decimal place on all processors and this is guaranteed by the BSP and cray SHMEM implementations and strongly recommended by the MPI standard Moldy relies on this behaviour to ensure that no divergence occurs If this were not the case it would be necessary to add a step to synchronize the co ordinates periodically This may also be an issue for running on networks of workstations it is assumed that the numerical properties of all the processors are identical The equality of the separate copies of the co ordinates is checked periodically and the simulation exits if divergence is detected 53 5 4 3 Input Output and startup Since each processor runs an identical copy of the program executable steps must be taken to ensure that input and output are handled by a single processor Otherwise a parallel run would print P copies of every line if running on P processors Indeed there are some parallel systems on which only one processor is allowed to perform input or output It follows that all parts of the code which perform I O must be aware of which processor this instance
55. a timestep matrix c Functions for manipulating 3 x 3 matrices and applying to 3 x n arrays of co ordinates etc output c Main output functions responsible for regular output sent to main output file Also contains error handling function message which may terminate simulation parallel c Interface between Moldy and various parallel message passing libraries Also contains functions to copy Moldy s data structures from the input output node to other nodes using the parallel interface quaterns c Functions for manipulating arrays of quaternions rdf c Contains init rdf which sets up RDF database rdf accum which periodically bins pair distances and rdf out which calculates and prints all the radial distribution functions restart c Functions for reading and writing the restart and backup files 43 start up default control re re header get saved params read_control conv_contro to input uni conv control read control get new params conv contro to prog units p backup lockfiles exist Create backup amp dump lock files backup file exists banner page write output file Functions called by start up default control Initialize struct cont rol with default parameter values read control Read parameters from control file and store in struct control convert control Convert physical values in cont rol betw
56. alculation seed integer 1234567 Seed for random number generator page width integer 132 Number of columns on output paper page length integer 44 Number of lines on a page of output mass unit real 1 6605402e 27 Unit of mass for system specification file length unit real le 10 Unit of length for system specification file time unit real le 13 Unit of time for system specification file charge unit real 1 60217733e 19 Unit of charge for system specification file where species name is the name of the molecule and N is the number of molecules of that type in the system Each molecule has n atoms one for each line in that group and each kind of atom is identified by a number id the site id which will be used to specify the appropriate potential parameters Its co ordinates are x z its mass is mj its charge is q and its name is name See Appendix A for some sample system specification files If there is more than one atom of any type in the system not just the same molecule it is sufficient to identify it by its id and the site co ordinates If 4 or name are given they must agree exactly with the previous values or an error will be signalled Site ids masses and charges are all checked for reasonableness and impossible values cause an error set of site ids does not have to start at 1 or be contiguous but since this may indicate a mistake a warning is issued Following the physical specificat
57. alculation see section Consequently the serial part becomes an increasingly insignificant fraction of the computation as the system size increases Since one might reasonably expect that in practice N is scaled with the number of processors the serial part does not seriously limit the parallel performance The advantages of the replicated data strategy are firstly simplicity of programming Most of the code is identical to the serial program which greatly eases maintenance Secondly load balancing is much more straightforward than with domain decomposition methods This will be discussed below Thirdly the communications granularity is coarse consisting of one global summation per timestep This makes it very suitable for loosely coupled processors such as a workstation cluster where the communications latency is high Finally it is very efficient for small and medium sized systems The big disadvantage is that the scaling to very large numbers of processors is poor which will eventually limit the size of system which can be simulated This is because both the amount of memory used per processor and the total amount of data to be communicated per processor increase linearly with N the number of particles in the system However in comparison with other supercomputing applications molecular dynamics simulation uses relatively little memory and most supercomputers can easily accommodate simulations of many thousands of particles Though scaling
58. als specified by m and initial time increment i are taken relative to the dump slices selected with t rather than the complete dump records stored in the dump files For example if dump records 1 3 5 19 are extracted using t then time intervals of m 1 5 correspond to real intervals of 2 4 10 relative to the original dump file data It should also be pointed out that the msds for each time interval are only averaged over as many initial time slices as the longest time interval Therefore it is usual to specify the largest value of m to be half that of the total number of dump slices in order to optimise the msd time interval and number of calculations over which each msd is averaged Selecting option u circumvents calculation of the msds and instead sends the connected trajectory coordinates to the output The default format for trajectory data is species name followed by a row of xyz for a single particle 61 for each increasing time interval Trajectory data for each particle are separated by two blank lines This is the format favoured by GNUPlot Specifying w 1 selects the generic matrix format in which a single row contains consecutive xyz coordinates for all particles with time increasing down the columns B 7 Mdvaf Mdvaf is a utility for calculating the Velocity Autocorrelation Function VACF from moldy dump files of at least level 2 see SectionB 7 It can also optionally calculate the Velocity Total Corre
59. array hk1 Each element of hk1 is a struct containing the values of l and the vector components of k The loop over k vectors then takes the form of a single loop over this list which is easily sliced over processors exactly as the loop over subcells above for phkl hkl4ithread phkl lt hkltnhkl phkl nthreads Compute the contribution of k vector in phk1 At the end of this loop the force arrays contain the contributions to the force on each site from the particular set of k vectors handled by that processor The total forces are obtained by globally summing these arrays over processors In fact since the same arrays are used to store both the short ranged and long ranged contributions to the force then only one global sum is needed for the forces Since the contribution at each k point is a sum of fixed length and takes exactly the same time to evaluate as any other this algorithm is always optimally load balanced There is also an implementation of the alternative RIL parallelization strategy in file ewald RIL c which is a direct replacement for ewald c In this case each processor handles a subset of the particles and all k vectors This involves an extra global sum within the loop over k vectors to evaluate the total structure factor for each 5 4 5 Radial Distribution Functions The accumulation of distances for the radial distribution function calculation is potentially costly since it is a site pair property
60. ase available to the restart file reading and writing functions in restart c Radial distribution functions are calculated between all pairs of site types defined by the site identifier given in the system specification file Storage for the accumulating binned pair distances is provided by a 3D pseudo array rdf idi idj bin with two indices for the site identifiers of the pair and one for the bin This does not have the geometry of a true 3D array but is set up by function init rdf to economise on memory by making rdf idi idj bin an alias for rdf idj idi bin Both the pointers and the actual data area are of course dy namically allocated The latter is in a single block of memory whose address is stored in rd base so that it may be easily stored in and retrieved from a restart file The function ptr is provided to make the size and base address of the RDF database available to the restart file reading and writing functions in restart c 5 22 Files and Functions 5 21 Source Files The source files consist of 8 header files and 23 code files which contain functions grouped a modular fashion A detailed description of the purpose and interface to each function is given in the source code in the form of comments which should be regarded as the primary documentation for the function An annotated list of files follows structs h Definitions of structure types used throughout the program 42 d
61. ased pseudo array The base pointer r is of type pointer to pointer to double and declared double r This points at a 1 dimensional array with length of pointers which in turn point to the rows of data arrays This structure contains all the information needed to access element i j using pointer indirection rather than an indexing computation and therefore without any reference to the values of m or n Higher dimensional arrays are set out in a similar fashion with more levels of pointer arrays n 1 for an n dimensional array There is one important instance of a more sophisticated use of arrays of pointers the allocation of the atomic site site force centre of mass force and torque arrays in the main timestep loop function do step These syntactically resemble 3D arrays declared e g site nspecies nsites 3 but because the number of sites differs between molec ular species they do not map onto rectangular 3D arrays However the array of pointers construct does allow for rows of different lengths and this is easily set up These pseudo arrays can again be passed straightforwardly as function arguments and used with the obvious and intuitive indexing scheme at the cost of a little extra code to allocate and initialize the pointers 5 1 3 System Variables and Data Structures Moldy makes consistent use of C language structures or structs to combine related pieces of data into a single variable This may then be passed a
62. ather than from the beginning The related control parameters are backup file which specifies the file name and backup interval which gives the frequency of backups It should not normally be necessary to change the name but the optimum interval will depend on the size of the simulated system and the speed of the computer By default it is 500 At the successful end of a run the backup file is deleted so that only if there is an abnormal termination does one remain The restart from a backup is entirely automatic If a backup file exists when a run is started it is read in and the run continues from it In contrast to a normal restart all of the control parameters are taken from the backup file and the control file and a restart file if one is specified is ignored f In consequence if a run is aborted or stops abnormally for some reason the backup file must be removed manually otherwise next time a run starts the unwanted simulation will continue instead If a run terminates abnormally there may also be a lock file called MDBACKUP Ick which ought to be removed Moldy attempts to prevent two runs from overwriting each other s backup files by creating a lock file whose name is formed Whether the old file is lost depends on the operating system Under systems such as VMS which have version numbers a new version is created and the old one remains Under Unix the old file is renamed by the addition of a character and thus is saved On other s
63. ature to the correct value while the system approaches its equilibrium state the kinetic energy is made to approach its desired value Scaling may be performed every timestep or every few depending on the simulation conditions An MD run with scaling does not generate a valid statistical ensemble and it must therefore be switched off before any calculation of thermodynamic averages is performed Moldy incorporates two refinements to the basic scaling algorithm which are selected by the control parameter scale options Linear and angular velocities can be scaled independently either for the whole system or for each species individually In this way one does not rely on the interactions of these degrees of freedom for convergence to equilibrium In many systems these degrees of freedom are loosely coupled and the exchange of energy between them is slow In these cases individual scaling can speed up the approach to equilibrium considerably The other refinement addresses the problem of setting the temperature accurately At equilibrium the system s kinetic energy fluctuates with mean square amplitude 8K 1 gt which corresponds to rms fluctuation in the instantaneous temperature 4 72 2 gT The difficulty with applying equation is the instantaneous kinetic energy K in the denominator Strictly scaling ought to use the average kinetic energy K as in equation 2 27 but this quantity is not known until after the run is co
64. aving However additional work must be done to create a fully single precision version since the standard double precision maths library will still be used The C standard does not require a single precision maths library but most systems do make the functions available However there is no standardization of the interface so this can not be done portably boolean Used for storing logical values Typed to int gptr Generic pointer type set to void in the case of ANSI C Pre ANSI compilers do not support void so char is used instead time_mt This is used to store integer format times and dates as returned by the time function It is declared as unsigned long and is used because pre ANSI compiling systems may not define the t ime t type for this pur pose l The term function in C corresponds to both functions and subroutines in FORTRAN These are the struct cont rol section and the integers ithread and nthreads holding the parallelization parameters section 5 4 3 39 size mt Used for storage of C object sizes as returned sizeof Like time mt this would be unnecessary if we were guaranteed an ANSI C compilation system vec mt Array of 3 reals for holding a vector type quat mt Array of four reals for holding a quaternion type mat mt 3 3 array of reals for storing a matrix 5 11 2 Memory Management The allocation of memory for storage of the multitude of data required in a molecular dynamics simulation is one of
65. bject OPT is used for less performance critical modules and may be set to a lower level of optimization than OPT2 There is just one other option to configure enable parallel which attempts to configure for building the parallel version The syntax is enable parallel parsys where parsys is one of mpi bsp shmem or tcgmsg If parsys is omitted the default is mpi See the section compiling for parallelism for a description of these libraries This is all that is required for a parallel build if the include files for the parallel system are in a directory searched by the C compiler and the libraries can be found using 1 The Cray T3D T3E and SGI Power Challenge and Origin 2000 systems fall into this category Otherwise it may be necessary to specify these by setting the environment variables CPPFLAGS and LIBS eg env CPPFLAGS I usr local mpi include LIBS L usr local mpi bin lmpi configure Some parallel systems such as the MPICH implementation of MPI from Argonne National Laboratories supply a compiler script called in this case mpicc If you are using MPICH it is sufficient to do env CC mpicc configure enable parallel Configure assumes that mpicc invokes the system C compiler and adds the appropriate optimization flags Other cases where an alternate compiler specification is used are for the SP2 where the command is mpcc and for the Oxford bsp library which provides bspcc VMS Si
66. ch exactly satisfy J T at all times One such is the Gaussian constraint see Allen and Tildesley 2 pp 230 231 The equations of motion are R irk v EF Cr Lvi Fy 2 31 1 1 0 bs Y Qj Nj 2 0 1 0 These bear superficial resemblance to equations but now is not itself a dynamical variable but a function of the other dynamical variables Note that does not enter explicitly into equations 2 31 since 7 is now a conserved quantity equal to its initial value at all subsequent times Perhaps surprisingly these equations of motion generate con figurational averages from the canonical ensemble but this does not hold true for the momenta and therefore dynamics nor for fluctuations 27 Constant Stress pressure It is frequently desirable to conduct simulations under conditions of constant pressure or stress rather than constant volume For example this allows the simulation of a solid state phase transition with a change of symmetry or unit cell size Moldy incorporates two constant stress methods which allow the MD cell to change size and shape and two variants which constrains the MD cell to uniform dilation to give an isobaric ensemble There has been some discussion in the literature of the validity of representing a heat bath by a single dynamical variable which the more diligent reader may wish to explore D BG 12 In a constant stress simulation the MD cell changes in size
67. cially FORTRAN ones Dumpext is a program which processes dump files extracts a subset of the information and outputs it in form more suitable for reading by other programs It is invoked with the command dumpext R nmols Q nquats c components t timeslices m molecules output file 1 0 dump filel dump file2 The meanings of its arguments are R the total number of molecules This argument is compulsory Q the number of polyatomic molecules This argument is compulsory c which piece component of the information in dump record to extract Cof M positions quaternions unit cell matrix potential energy Cof M velocities quaternion derivatives unit cell velocities C of M accelerations quaternion accelerations 10 unit cell accelerations 11 C of M forces 12 torques 13 stress tensor This argument is compulsory t which timeslices or dump records to extract This is a selector format list and is the index in the whole dump sequence not just a single file Timeslices or dump records are sequentially numbered in the dump files from 1 Defaults to all timeslices m extract information for selected molecules This is a selector list specifying the molecule index Defaults to all molecules o of optional output file Defaults to standard output b selects binary output in single precision floating point numbers Defaults to ASCII formatt
68. control file to specify energy units other than kJ mol See section 3 2 1 At present ransub only introduce monatomic solute species but polyatomic species can easily be added manually by modifying the system specification file created by ransub It is also possible to introduce randomly distributed vacancies by prefixing the solute name with 64 B 11 Plotrdf Plotrdf is a perl script which reads Moldy output files and extracts RDF information skipping page breaks It produces output suitable for reading into most graphics programs but most notably gnuplot and xmgrac It can also write a command script for either of these two programs if requested Obviously you must have perl installed on your computer to use it plotrdf a gp xg ps eps np m REGEXP moldy output file a average of multiple RDF blocks gp Write a script and invoke GNUPLOT to plot data xg Write a script and invoke GRACE to plot data ps Invoke GNUPLOT or GRACE to plot data and write as a PostScript file eps Invoke GNUPLOT or GRACE to plot data and write as an encapsulated PostScript EPS file np Do not plot data write a GNUPLOT or GRACE script m output plot only columns which match the perl regular expression REGEXP no needed The output format is in columns with different components laid out across the page r rdf1 rdf2 rdf3 Multiple RDF blocks in one or multiple files are processed and outpu
69. cular dynamics to begin 5 3 2 Main timestep loop The initial entry point into Moldy is into function main which calls start up section 5 3 1 and then proceeds into the main loop over molecular dynamics timesteps Figure D 3 This loop calls do step which computes the forces and advances the co ordinates velocities and their angular counterparts by one timestep Most of the periodic tasks not directly associated with progressing the simulation such as the accumulation of running averages velocity rescaling and writing the output files are called directly from main The exceptions are the binning of site pair distances for computing the radial distribution functions which is integrated into the site forces evaluation and the writing of trajectory data to the dump file which is called from do step since the force arrays are local to that function The timestep loop continues until either the requested number of steps have been completed a SIGTERM or SIGXCPU signal is received or one more step would exceed the CPU time limit set in cont rol cpu limit Only the normal termi nation case is shown in Figure 5 3 In the case of an early termination write restart is called to write out a backup file so the simulation may be restarted later The signal handling functions shutdown for SIGTERM and SIGXCPU and siglock for other signals are installed in main following the call to start up Function shutdown simply sets a flag whic
70. cule The system specification should contain an entry similar to that for a normal molecule which describes the atomic sites belonging to one MD cell s worth of the framework Its periodic images should fill space to construct the required infinite framework This is notified to the program by modifying the first line of the specification of that molecule to read species name 1 Framework compare with section B Z T The effect of the special keyword framework is 1 to remove the rotational freedom of the molecule This preserves the infinite structure over MD cell repeats by disallowing relative motion of its parts Linear motion does not destroy the structure and is allowed 2 to modify the effect of the periodic boundary conditions Normally a molecule is assumed to be small and periodic relocations are applied to all of its atoms depending on its centre of mass co ordinates relative to some interacting molecule In contrast the atoms of a framework are independently relocated This ensures that each molecule sees all framework atoms from any unit cell which are within the cut off distance 30 In the present version of the program only one framework molecule is allowed though more may be permitted in future versions Consequently the configuration given as a lattice start must fill the entire MD box A skew start is not sensible under these circumstances since the orientation of the framework must be explicitly specifie
71. d from the standard input which may be a terminal or a job command file depending on the operating system and circumstances and the output is written to standard output which may be a terminal or batch job logfile Here are examples for VAX VMS and Unix tm which assume that in each case the command has been set up to invoke Moldy Under VMS the commands moldy control dat output lis moldy control dat will start Moldy which will read its input from control dat The output will be directed to the file output lis in the first case and written to the terminal or batch log in the second Under UNIX any of moldy lt control gt output lis moldy control output lis moldy control oe will cause moldy to read from the file called control and in the first two examples to write its output to output lis The command line interface for Windows and MS DOS is very similar 3 1 The Control File The information needed to initiate and control a run of Moldy is specified in a file known as the control file This contains the parameters governing the run e g the number of timesteps to be executed or the frequency of output and the names of files to be used e g for reading a restart configuration from or for writing the output to Parameters in the control file are specified by entries of the form keyword value which appear one to a line terminated by the special keyword end Spaces and blank lines are ignored as are comments
72. d to construct a good space filling structure 3 11 Messages and Errors Apart from the periodic output there are occasional once off messages which Moldy writes to the usual output file Such messages begin with the characters I W E or F denoting the classes information warning error or fatal respectively Their meanings are I information These are often produced by subroutines to give useful information on their particular calculations For example when temperature scaling is turned off a message to that effect is recorded in the output file Various routines which calculate internal quantities such as the Ewald sum self energies and distant potential corrections also record their values using an information message W warning When the system specification is suspicious but not clearly wrong or some untoward condition is detected such as two atoms approaching too closely a warning message is issued E error Error messages are issued when a mistake is detected reading any of the input files To make correction easier processing continues until the end of that file so that all of the errors are found The simulation is then stopped with a fatal error F fatal The simulation is terminated immediately Faulty input files generate fatal errors after they have been com pletely processed There are many other conditions which also generate fatal errors for example if the simulation algorithms violate some criterion such a
73. di Fro QnkGnX 47280 i 1 n 1k 1A K 1 ral Point self energy Intra molecular self energy 1 Oe AD i Zo ifi o2 6 9V ha Charged system term Surface dipole term lattice vector of periodic array of MD cell images reciprocal lattice vector of periodic array of MD cell images modulus of k absolute indices of all charged sites index of molecules indices of sites within a single molecule total number of charged sites total number of molecules number of sites on molecule m co ord of site i relative to molecular centre of mass charge on absolute site i charge on site of molecule m Cartesian co ordinate of site i Tj Fi real reciprocal space partition parameter instantaneous stress tensor Kronecker delta symbol xyz tensor indices volume of MD cell and the force on charge i is given by fi P ya Q rijtn 20 dia rijtn ruta 17 rij n Sa Real space term bol 2 k zy l4 P k rj cos k r 2 2 20 j l Reciprocal space term Surface dipole term The molecular forces and torques F and N are evaluated from the site forces f using equations E T and Notice that the equation 2 T9for the energy contains a correction for the intra molecular self energy whose derivative is absent from the equation for the forces equation 2 20 This term corrects for interactions between charges on the sam
74. dient using the Einstein relation lt r t r 0 gt 6Dt B 1 see Allen and Tildesley 2 p60 For a species of N particles the mean square displacement is calculated as 1 N lt r t r 0 gt NN Y Yn 0 r9 B 2 t n 1 to where r t is the position of particle n at time t In interactive mode you are prompted for the source of the system specification sys spec or restart file and the dump file for reading the configurational information You will also be prompted to supply the dump range limits and msd time interval limits before the calculation can commence Be careful to either redirect the output using gt on unix or specify an output file with o otherwise the output will go to the screen More options can be specified in command line mode msd s system specification r restart file d dump file format t dump range m msd time range i initial time increment g species range u w trajectory output format x X yl Y 2 21 7o output file where meanings of the options are s read the species data from a system specification file r read the species data from a restart file Only one of s or r may be given d read the centre of mass coordinates from a dump file given as a prototype name containing a printf format string see section 7 t range of records in dump file specified in selector format m time intervals to calculate the msd values for
75. dy writes a banner page containing a summary of the system being simulated and details of the important control parameters The bulk of the output file is the periodic output which contains the instanta neous values of various thermodynamic variables their rolling averages and associated standard deviations The rolling average is just the mean over the preceding n timesteps where n is set by the control parameter roll interval an notated example is given in figure B I The frequency of periodic output may be altered by setting the control parameter print interval to the interval required This may be necessary to constrain the size of the output file which can grow to be very large indeed with the default interval of only 10 As well as the short term rolling averages long term averages are calculated and printed out at regular but usually infrequent intervals Accumulation starts on the timestep given by the control parameter begin average and every average interval timesteps thereafter the means and standard deviations are calculated and printed This output is interspersed with the periodic output and is formatted with one variable to a line in the form mean s d Where a variable has more than one component such as multiple species for the translational temperature or Cartesian compo nents for the mean square forces the components are printed across the page In addition to those variables printed as part of the periodic output the pres
76. e drive Large disk stores are now very cheap so this should not be a problem in practice If disk store is limited then the simulation may be divided into multiple Moldy runs interspersed with copying of dump files to tape Notice that Moldy must sometimes read an existing but complete dump file to propagate the unique marker to all of the files in a sequence Therefore when continuing a simulation and a dump run at least the immediately preceding dump file must still be accessible This should be borne in mind when copying dumps to tape Moldy is careful to ensure that existing files are not overwritten especially necessary since dump records are added to the end of an existing dump file Whenever Moldy prepares to start a new dump file it checks to see if one of that name is already present If so a new name is chosen by mutating the old one and a warning message to that effect is written to the output file On the other hand if the first file of a new dump run including one initiated because of some error in continuing an old one already exists the prototype file name is mutated as above and the whole dump run is written to files based on the mutated name When a run is restarted checks are made to ensure that the values of the dump control parameters have not been altered If they have it is not possible to continue an existing dump sequence and a new one will be started If existing dump files are present the new sequence will have muta
77. e molecule which are implicitly included in the reciprocal space sum but are not required in the rigid molecule model Though the site forces f do therefore include unwanted terms these sum to zero in the evaluation of the molecular centre of mass forces and torques equations Z T and by the conservation laws for linear and angular momentum 2 41 Parameter Values Both the real and reciprocal space series the sums over n and k converge fairly rapidly so that only a few terms need be evaluated We define the cut off distances and ke so that only terms with r 4 n lt re and k lt ke are included The parameter determines how rapidly the terms decrease and the values of and needed to achieve a given accuracy For a fixed and accuracy the number of terms in the real space sum is proportional to the total number of sites N but the cost of the reciprocal space sum increases as N An overall scaling of 3 be achieved if varies with This is discussed in detail in an excellent article by David Fincham I6 The optimal value of is I2 Q 21 where fg and fr are the execution times needed to evaluate a single term in the real and reciprocal space sums respec tively If we require that the sums converge to an accuracy of exp p the cutoffs are then given by te 2 22 ke 2 4 2 23 representative value of specific to Moldy has been established as 5 5 Though thi
78. e compiler must therefore assume that if an array element is written in a loop then elements of any other arrays may also be changed It must therefore reload from memory even though it already has a copy of the value in a register But if all loads are completed before any stores then the compiler is at liberty to re use the values and save memory accesses 36 Most compilers also have option which directs it to ignore recurrences throughout the whole program e g va on the Convex va and h ivdep on the Cray compilers It should normally be safe to use these options Each manu facturer s compiler has its own peculiar set of inline directives For example the CRAY compilers use a pragma ivdep statement whereas the convex and Stellar compilers use a significant comment dir no recurrence E 4 5 Modifying Moldy 4 3 1 Adding a New Potential By default Moldy supports potential functions of the Lennard Jones six exp and MCY forms However it should be very easy to add further types The program is written in a highly modular fashion so that the only code which need be altered is in file kernel c and occasionally in defs h The calculation of the potential and forces is performed entirely in the function kernel This function is called repeatedly with a vector of squared distances between some reference site and its neighbour sites Vectors of potential parameters and charges are supplied which bear a one to one correspondence wit
79. e in modern heterogeneous networks where diverse machines frequently share a common file space The format is based on Sun Microsystems XDR protocol 24 The XDR routines are available on almost every modern Unix machine and are simple enough to implement on any other system If the control parameter xdr is set to 1 then all files will be written using this format Moldy automatically determines whether a restart file was written using XDR by examining the file header and reads it in the appropriate fashion irrespective of the value of xdr 3 4 Setting the Temperature Moldy implements three different methods to control the temperature of the simulated system These are the veloc ity rescaling technique described in section 6 1 the Nos Hoover thermostat and constrained equations of motion section D 6 2 Scaling is selected by the parameters scale interval etc Every scale interval timesteps un til scale end the velocities are adjusted so that the kinetic energy corresponds exactly to the desired temperature the value of control parameter temperature The Nos Hoover and constrained thermostat are selected by setting const temp equal to 1 or 2 respectively The control parameter scale options selects refinements to the basic scaling or thermostat algorithms This is an integer parameter interpreted as a set of bit flags with the following meanings bit 0 perform scaling or thermostatting for each molecular species individually bit 1 scal
80. e simulation timestep in ps and nsteps which specifies the number of steps to perform Together these control the length of time to be simulated It is also possible to specify that a run should be terminated after a certain amount of computer time has been used given by parameter cpu limit This will be particularly useful in batch mode systems where the run is killed after some specified CPU time has elapsed Setting cpu limit to the maximum time allowed will cause Moldy to terminate the run before the limit is reached and write out a backup file see section 3 3 1 There are several kinds of parameters character strings Apart from t itle these are just file names e g sys spec file No checks are performed on the validity of the name because Moldy has to work on many different computer systems so if you make a mistake you are likely to get an error message to the effect that Moldy couldn t find the file you asked for To remove a default value just specify a null string e g save file booleans These are just switches which turn a feature off or on 0 means off or false and 1 means on or true The parameters text mode save new sys spec surface dipole and lattice start are booleans real parameters Several parameters are real numbers e g step which specifies the timestep They may be entered in the usual floating point or scientific notation e g step 0 0005 or step 5e 3 and are taken to be in the units given in table
81. e thermostat the rotational and translational components of the kinetic energy separately bit 2 use the rolling averages of kinetic energy to calculate the scale factor rather than the instantaneous values bit 3 discard all existing velocities and accelerations and re initialize from the Maxwell Boltzmann distribution The bits may be set in any combination so for example scale options 6 sets bits 1 and 2 6 2 22 and scales separately for rotation translation using the rolling averages If bit 3 is set the others are ignored Only bits and 1 have any meaning in the case of a thermostat and signify that each species or the translational and rotational degrees of freedom are isolated from each other and coupled to their own individual heat baths The options for scaling separately rotation and translation and per species may be useful for achieving equilibration in difficult systems where mode coupling is ineffective In those situations it is otherwise possible for all the energy to be transferred into the rotational modes of a particular species halting any progress to equilibrium for other degrees of freedom These options ensure that all degrees of freedom have some thermal energy The option controlled by bit 3 to discard all existing information and start from a random set of velocities may be of use when starting from far from equilibrium situations In such cases the forces are frequently so large that the velocities and acceleratio
82. e value of sys spec file This file is divided into two sections First is the description of the molecules atoms or ions which is followed by the potential functions As for the control file the input is case independent and free format but line structured Blank lines spacing and comments are ignored The physical description consists of a series of entries one for each molecular species terminated by the keyword end The entry for species i should have the form species name id x yi zi my qj name id 9 yo zo mo q name ida yn Zn qn namen 20 Table 3 1 Control Parameters name type default function title character Test Simulation A title to be printed on all output nsteps integer 0 Number of MD steps to execute cpu limit real 1 20 Terminate run if excessive CPU time used step real 0 005 Size of timestep in ps sys spec file character null Name of system specification file Appended to control file if null lattice start boolean false Switch for crystalline initial configuration save file character null File to save restart configuration in restart file character null File to read restart configuration from new sys spec boolean false Read restart configuration with changed system specification text mode save boolean false Write a portable restart file consisting of control system specification and lattice start files density real 1 0 I
83. ed in the header file defs h If a new potential with more than 7 parameters is to be added to Moldy it will be necessary to increase its value The above primary storage of the potential parameters is convenient for reading in and storing and rereading restart files However it can not be indexed efficiently to retrieve the potential parameter values in the critical innermost loops The data is therefore copied into an expanded 3 dimensional array max id NPOTP nsites in the force evaluation function forces before being used The masses and charges of the various sites are stored in another array of structs site info simply indexed by the site identifier As with the potential parameters this is not convenient for access within the inner force loop so the data are expanded in do step to fill an array chg nsites to allow direct indexing by the loop counter The two final data structures of importance are those containing accumulated data for computing the usual averages and the radial distribution functions Both databases are considered private to their respective source modules values c and rdf c and are only accessed using the function interfaces provided in those modules The averages database provides a general and extensible scheme for storing partial sums and computing rolling and time averages of instantaneous values such as temperature kinetic energy etc It consists of two parts a linear array av of structs of type av
84. ed memory multiprocessors and even networks of workstations Of course formal portability is no guarantee of good performance on hardware with widely differing processor and communications performance That is a function of the granularity of the problem with respect to communications latency and volume of data with respect to communications bandwidth Both of these depend on the nature of the problem and the parallelization strategy 5 4 1 Parallel Strategy Moldy uses the replicated data approach to parallel molecular dynamics 10 50 whereby each processor maintains a complete set of dynamical variable arrays for all particles The short ranged and reciprocal space force calculations are divided among the processors and a global sum operation is used to add the separate contributions to the forces and propagate a copy of the complete force array to every processor This parallel part of the calculation dominates the execution time and shows a close to linear speed up It takes over 9596 of the CPU time even for quite small systems of a few hundred particles The remainder of the calculation in particular the integration of the equations of motion is not parallelized but is performed redundantly on every processor simultaneously This serial part of the code will therefore limit the parallel speedup on large numbers of processors by Amdahl s law However this part of the computational work only grows as N compared with N for the forces c
85. ed numbers If any of the compulsory arguments are omitted you will be prompted to supply a value In particular you must always supply the number of molecules and polyatomic molecules This information is not recorded in the dump header and is needed to determine the number of co ordinates and quaternions in each dump record You must specify which pieces of information to extract using the c option The dump files must of course be part of the same dump sequence this is carefully checked They may be supplied as arguments in any order dumpext automatically determines the sequence from the information in the headers This is not as pointless as it sounds since the list may be generated using unix shell wild cards which arrange them in alphabetical rather than numeric order The output is arranged in columns one line per time slice So if for example you wish to extract positions and quaternions there will be 7 columns corresponding to 2 40 41 42 43 Multiple components are printed in the in teger order of the component not the order specified with c If binary output is asked for with the b option the 59 order is the same The numbers written sequentially as single precision floating point numbers in machine tive format without record delimiters The records may be easily read in C with an read call with the appropri ate number of bytes or from FORTRAN77 using direct access unformatted read of the appropriate l
86. een input and program units read sysdef Read the system specification file Allocate storage for arrays species site info and potpar and copy in values from file initialise sysdef Finish set up of system and species structs allocate dynamics Allocate memory for the dynamic variable arrays and set up pointers in system and species Skew start Setup skew cyclic initial state lattice start Read crystal structure and set up co ordinates quaternions and cell vectors init cutoffs Determine optimum parameters for Ewald sum re re header Readthe restart header and control structs from the restart file re re sysdef Read the system specification structs arrays system species site info and potpar from the restart file read restart Read dynamic variables and the RDF and averages databases from the restart file check sysdef Check new system specification is consistent with old thermalise Setup Maxwell Boltzmann velocity distribution init rdf Prepare to bin RDFs Allocate memory and pointers init averages Allocate space for and initialize the averages database convert averages Update averages database if roll interval changed or if old style restart file conv potentials Convert potential params between input and program units Figure 5 2 a Block diagram of the initialization function start up and a list of the functions called Continued in Figure 5 206 startup
87. efs h Main Moldy header file to be included in all files It contains machine and operating system configuration macros physical constants unit definitions and global non structure type definitions string h time h stddef h stdlib h These files are replacements for the ANSI C library header files of the same name included here for portability to pre ANSI environments messages h Error messages file containing messages in the form of preprocessor macros This may allow for non English language customization xdr h Includes the system headers for the usual External Data Representation XDR functions and prototypes of addi tional functions for XDR reading and writing of Moldy specific structure types accel c Contains do step which implements the major part of the main MD timestep procedure plus functions for velocity rescaling and thermostatting algorith c Contains functions to implement algorithms and computations related to the dynamical equations of motion This includes the computation of molecular centre of mass forces torques the Newton Euler equations and the constant stress and temperature functions These functions all have an interface which makes no reference to the system species structs alloc c Contains functions for memory management allocation 11 freeing t free and setting up of pointer based multi dimensional arrays arralloc ansi c Replacements for certain functions r
88. ength records OPEN ACCESS DIRECT LRECL nnn In the above example there are 3 x 4 4 4 x 4 28 bytes in each record B 5 Mdshak Mdxyz Mdshak was originally written as an interface between Moldy configurations and a molecular graphics program called SCHAKAL 2G6 The current version can write the output in additional formats Brookhaven PDB Cambridge CSSR or XYZ form suitable for input to the freely available viewers such as and RasMol There is also a viewer called crystal written by the author for use with AVS the general purpose visualization program and available free from the International AVS Centre Mdshak writes a file containing the co ordinates of a single or of multiple time slices during a MD run It can read configurations from a a system specification file plus lattice start b a system specification file plus dump file c a restart file or d a restart file plus dump file It can be driven either interactively or using command line arguments If the executable file is named mdshak then the default output is in SCHAKAL file format If it is renamed to mdxyz mdcssr or mddcd then the default output format is of the corresponding type To assist in making animations of a simulation the trajectory from many timesteps can be written to a file in the CHARMM and Xplor DCD format using the v option In conjunction with a PDB file this may be used as input to the VMD visualization program To
89. epulsion in support of ionic solution models It takes the form 12 A 8 rij Aijexp Bijrij Cij rij Di Ti Ey rf with the six parameters Fij morse This is the Busing Ida Gilbert plus Morse potential as used in the Japanese MD program Mdxorto It has the form rij Bijexp Cij rij Dij Fj exp 2Gi rij Hij 2Fijexp Gi rij Hij with seven parameters B HIW This is a Lennard Jones potential with an extra fourth power term to describe ion water interactions in solution using the HIW model of R R Pappalardo and E Sanchez Marcos B8 The form is 4 6 12 rij Ai rij Bij rij Cijfrij Other types of potential may be easily added see section 3 1 It is possible to specify the units in which these quantities are given by means of the control file parameters mass unit length unit time unit and charge unit which are themselves specified in SI units poten tial parameters site masses and charges but not co ordinates are taken to be in those units Their default values are amu 0 1 and qe which means that the unit of energy is kJ mol So to read A amu and kcal mol spec ify time unit 4 8888213e 14 for amu and eV use time unit 1 0181e 14 and for Lennard Jones potentials specified in A amu and in K ie use time unit 0 548343e 12 which includes the extra factor of 4 Once the system specification has been read in all quantities a
90. equired in any ANSI C standard library but missing from pre ANSI environ ments auxil c Machine and OS dependent support functions Contains a Fast vector arithmetic functions including interface to various vector libraries as well as pure versions b OS dependent time and file manipulation functions beeman c Functions implementing the separate stages of the modified Beeman algorithm including basic steps and the updating of all the necessary dynamical variables convert c Functions for converting the potential parameters and control file parameters from input units to program units and vice versa dump c Functions for managing periodic trajectory etc dumps ewald c Evaluates the reciprocal space part of the Ewald sum for the long ranged Coulombic forces force c Implements the calculation of the short ranged forces and the real space part of the Ewald sum using the Link Cell algorithm This excludes the actual evaluation of the potential which is contained in kernel in kernel c input c Functions for reading the control and system specification files and allocating the structures to hold the data read in from them eigens c Matrix diagonalizer from the netlib archive kernel c Contains kernel which actually evaluates the potential given a vector of pair distances main c Function main is the main program which controls the set up and contains the main MD timestep loop that calls do step to do most of the work of
91. er near the top This means that there is no longer a unique molecular centre of mass and equations B T5 and 2 24 can longer be used to compute the internal stress or pressure On the other hand entire molecules may be assigned to subcells on the basis of their centre of mass co ordinates atomic sites belonging to a molecule are placed in the list belonging to the particular cell containing the molecular centre of mass This enables intramolecular interactions to be excluded from the calculation and equations 2 15 and 2 24 can be used to compute the internal stress There is however a computational disadvantage which arises because sites belonging to a large molecule may lie outside the subcell In order to correctly include site site interactions within the cutoff the list of neighbour cells must be extended to include all cells where any parts are closer than re 27m where is the molecular radius There is an inevitable increase in computer time required to process the lists and in the case of a strict cutoff many unnecessary intra site interactions must be evaluated In the case of a system containing a mixture of perhaps a few very large molecules and other small molecules the computer time penalty may be very large indeed Versions of Moldy prior to 2 16 used this molecular binning approach and suffered from the associated compu tational penalty Versions 2 16 and later adopt the former site binning approach whic
92. es to the corresponding mem ber of struct cont rol The default values were previously set by default_control In the case of a restart the saved parameters in the restart file are restored to control by function re re header overwriting all of the current values In this way the saved values become the defaults for a second call of read control which rereads the control file and assigns any new values The repeated read is necessary because the name of the restart file is supplied by the control file Note that those parameters representing physical quantities must be converted to and fro between input and program units for consistency since they are stored in program units in the restart file The three alternative routes within start up for reading the system specification and initializing the dynamic variables etc are shown in Figure The case of a new run is reasonably straightforward Memory for the system specification arrays species potpar and site info is allocated as the system specification is read in by read sysdef The raw atomic site co ordinates are then shifted to the molecular centre of mass frame and rotated to the principal frame of the inertia tensor by initialise_sysdef which also computes a number of other molecular quantities and completes the set up of the system specification Next the dynamic variable arrays are allocated and initialized and finally the RDF and averages databases are installed If the run
93. esent it is not possible to specify an anisotropic external stress though this capability may be added in future versions of the program The matrix may be constrained so as to disable motion of any or all of its components using the parameter strain mask Strain mask is a bitmap each bit of the integer freezes one of the components of h bit i freezes with i 3 k 1 1 1 The bitmask is the sum of 2 over the i s to be set so the strain mask values 1 2 4 h hy 8 16 32 constrain the corresponding components of h hn hy 64 128 256 h31 h32 h33 Thus the default constraint of ho a a4 0 0 h b bx by 0 is given by strain mask 200 84644128 Another useful value is 238 which freezes all the off diagonal components This is needed for a liquid simulation since there are no shear restoring forces acting on those components 3 9 Cutoffs and Adjustable Parameters There are four parameters related to the Ewald sum method of force evaluation see section 2 4 re and subcell1 In addition the two options strict cutoff and surface dipole select how rigorously the real space cutoff is applied and whether to include the De Leuuw surface dipole term By default and ke are chosen automatically using equations 2 23 to give a default accuracy of exp p 1075 ie p 11 5 in the Coulombic potential energy evaluation An empirically determined value of tg tr 5 5 is used If a diffe
94. estore the default par begin int argc char argv int ithread int nthreads Initialize the library and return the number of processes and the ID of this process par finish void Terminate the parallel run normally par abort int code Terminate the run abnormally Return code if possible par broadcast void buf int n size mt size int ifrom Broadcast the specified buffer from node ifrom to all nodes par_ r d i sum void buf int n Perform a global parallel sum reduction on the buffer containing n reals dou bles or ints par imax int idat Perform a global maximum reduction on the single int argument The SPMD parallel strategy updates the dynamic variables independently on all processors see section 4 1 To update the separate copies of the co ordinates and other dynamic variables synchronously the results of the floating point arithmetic must be identical Therefore the result returned by rsum and dsum must be identical to the last bit on all processors see the footnote on page 63 Another consequence is that execution on heterogeneous workstation networks is not supported the identity of floating point operations in not guaranteed even if all use IEEE arithmetic Moldy periodically tests for divergence of trajectories and will exit with the message F Trajectories on parallel threads are diverging if this condition is detected real is a typedef defined in defs h which is set either to float or
95. et electric charge Although the system as a whole may be electrically neutral the omission of framework framework interactions from the calculation also means that the k 0 term does not vanish To see this examine equation 2 19 In the non framework case the indices i and j in the terms 3 N gig k ri cos k r 1 1 M gt N qicos k r i l i run over all site pairs If k 0 the squared sum is x 2 0 i l If a framework is present the formulation becomes N 2 N 2 qicos k r qicos k r i l i M 1 assuming sites 1 N are the framework sites On setting k 0 this reduces to N 2 3 qi qi i 1 i M 1 It is therefore necessary to modify the charged system term of equation 2 19 to RD i M 1 2 55 Ya i 1 17 Ew 18 Chapter 3 Running Moldy The way Moldy is invoked depends to some extent on the operating system but usually by issuing the command moldy fj For Unix tm Windows 95 and NT and MS DOS the executable file moldy or MOLDY EXE should be placed in the shell search path e g in the current directory There is no GUI and Moldy should always be run from a command line from a terminal window under Unix or a MS DOS window under Windows 95 NT There are two optional arguments the name of the control file see section and the output file see section 3 5 If either is omitted control input is rea
96. every iteration the index of a molecule on each iteration depends on the index of the previous iteration This kind of loop is inherently un vectorizable which grossly limits the efficiency on vector computers To overcome this limitation a further step is needed as suggested by Heyes and Smith 22 Inside the loop over subcells the linked list is pre scanned and used to construct an array containing the indices of the molecules in the list This is known as the neighbour list since it points to the molecules in the region of the subcell under consideration The co ordinates are then assembled into a contiguous array using a gather operation The inner force loop then becomes a straightforward DO loop over this temporary array which is vectorizable In addition to the co ordinates a corresponding gather operation using the same neighbour list operates on the electric charge and potential parameter arrays The corresponding temporary arrays are also accessed in the inner loop over neighbour sites with a single array subscript Finally the computed forces are added to the main site force arrays using a scatter operation inverse to the gather Although this algorithm was designed for vector machines it is also highly efficient on modern RISC processors 51 which universally employ a cache based memory architecture Since adjacent iterations of the inner loop access adjacent locations in memory loading a cache line will bring the operands of the next fe
97. extract information from the instantaneous section of the output not the rolling averages or standard deviations B 2 Dumpanal Dumpanal examines the dump files given as arguments and prints out the header information to help with identification For example dumpanal dump filel dump file2 B 3 Dumpconv Dumpconv is a tool for moving binary dump files between computers of different architectures It has mostly been superseded by the portable XDR format dump files introduced in version 2 1 see section but is retained in case of machines for which no XDR implementation is available The command dumpconv binary dump file text dump file creates a straightforward ASCII text file with a representation of the dump information in the binary file including the header The command dumpconv d fext dump file binary dump file 58 converts it back Seven significant decimals are used which ought to retain almost all of the precision of the single precision binary floating point numbers You can also convert an old native format dump into format using the x option viz dumpconv x native dump file xdr dump file The x and d options may be combined if the input is a text format dump 4 Dumpext Dumps are designed so that Moldy can take care of all the bookkeeping perform data security checks and to divide a lot of data into manageable portions It is therefore not in a convenient form for reading by other programs espe
98. f molecules It also contains pointers to the array of site identification numbers and Cartesian site co ordinates expressed in the molecular principal frame for this species and also to the arrays of dynamical variables This is analogous to a COMMON block in FORTRAN 5strictly speaking molecules with rotational degrees of freedom 41 The dynamical variables are stored in dynamically allocated arrays of total length nmols and nmols r which may be regarded as the concatenation of individual arrays containing variables belonging to each molecular species The pointers in species i locate the beginning of the array of variables belonging to species i and those in system to the beginning of the entire array The variables are manipulated by either route as is most convenient but never both within a single function The potential parameters are stored separately in the array potpar since they refer to pairs of site types site id s This is an array of structs of type mt laid out as a 2 dimensional symmetric matrix of dimension system max id with rows and columns indexed by the site identifiers as given in the system specification file In fact it is stored in a one dimensional array of length system max id and the index calculation is performed explicitly when it is referenced Structure type pot mt contains an array p for the actual parameter values This is a fixed length array with NPOTP members where the constant NPOTP is defin
99. f pseudo Monte Carlo equilibration for difficult cases 2 6 2 Nos Hoover Thermostat A more sophisticated approach than rescaling is to couple the system to a heat bath Nos B5 proposed an extended system Hamiltonian to represent the degrees of freedom of the thermal reservoir Moldy uses the simpler but equivalent formulation by Hoover 23 The equations of motion equations 2 3 and 2 4 are modified thus This formula actually applies to the Canonical rather than the microcanonical ensemble but it serves for the purpose of this argument 11 Ri mR Ij Q 0j xl 0 8j 2 30 g T ke where most of the symbols have their previous meanings g is the number of degrees of freedom in the system G is a new heat bath dynamic variable and is the associated fictitious mass parameter With a suitable choice of Q these equations generate trajectories which sample the canonical 83 In other words both averages and fluctuations calculated as time averages from a simulation run tend to their canonical ensemble limits This does not necessarily guarantee the correctness of dynamical quantities It is apparent from equations that G has the dimensions of a time derivative and is analogous to the unit cell velocities in equations It is therefore treated as if it was a velocity and updated using only steps ii v of the Beeman integration procedure equations Z Z Note that the equation
100. f these are installed on your machine some public domain implementations are available for workstation clusters multiprocessors and many distributed memory parallel machines MPI The MPICH implementation can be downloaded 2 anonymous ftp from the BSP s Oxford BSP Library is available through Oxford Parallel s WWW server or by anonymous ftp SHMEM The CRI native communications library for T3D and T3E systems Alternatively a port to another parallel interface should be quite straightforward see section 3 2 Once a suitable message passing library is installed the procedure for building Moldy is as described in section ET On a non unix like system the configuration must be done by hand using whatever build mechanism is available The C preprocessor macro SPMD must be defined as well as one of MPI TCGMSG BSP or SHMEM This is usually done in the makefile by setting the Make macro PARLIBC DSPMD DMPI for example This macro should also include a I directive specifying the directory for the library s header files if these are not in the default path searched by the compiler The similar make macro PARLIBL should contain the linker directives necessary to link to the library itself Examples are provided at the top of the supplied Makefile This parallel implementation makes use of the replicated data approach BU whereby every processor has a complete copy of all the arrays containing dynamical variables and every site on every
101. gement the list of neighbouring subcells is fixed and may be precomputed Its construction takes only O N operations and only O N pair interactions need be calculated One drawback of the link cell method has been the difficulty of implementing it efficiently for vector processors The linked list of neighbour particles is not stored in the regular stride array which is required for vectorization Heyes and Smith 22 pointed out that gather operations might be used to assemble a temporary array of neighbour particle co ordinates from which the interaction potential and forces could be evaluated in vector mode A scatter operation is then used to add the resulting forces to the total force array This is the technique used in Moldy Almost all modern vector machines have scatter gather hardware which means these operations are fairly cheap 2 5 1 Molecular Aspects The treatment of molecules as rigid bodies adds an additional complication to link cell and neighbour list algorithms If atomic sites are assigned to link cells then there is no obvious way of determining in the innermost force loop whether two sites belong to the same or to different molecules It is therefore not possible to exclude intra molecular site site interactions from the force and potential energy calculation Furthermore molecules may be split by the periodic boundaries for example half of the sites of some molecule may reside near the bottom of the MD box and the remaind
102. gy Note that these messages only occur for polyatomic molecules If the system is monatomic the errors mentioned above may still be present but will not be detected and the simulation may fail in some less predictable manner for example particles may approach too closely and or acquire high velocities and fly off to infinity The message W Sites n and m closer than 0 5A gives some warning of this condition 31 Chapter 4 Compiling Modifying Moldy The Moldy distribution consists of numerous files of source code for Moldy and the utility programs command or job files to compile the source input for the manual and example control and system specification files For ease of transport these are packed into one large archive file whose format and method of unpacking depends on the operating system of the target machine At present it is available for unix The archive is usually a tar archive called moldy tar possibly compressed with the compress or gzip programs and named moldy tar Z or moldy tar gz These files may be uncompressed using the gunzip or uncompress programs viz gunzip moldy tar gz or uncompress moldy tar Z whereupon the archive is unpacked by the command tar xvf moldy tar An alternative form of archive for those unusual systems without a tar program is a shell archive a Bourne shell script called moldy shar This is unpacked by bin sh moldy shar VMS The archive is a
103. h is checked at the end of every iteration of the timestep loop This allows for an orderly shutdown if the CPU limit is exceeded or if the run must be interrupted for some other reason For conditions which require immediate program exit function siglock is called to delete the lock files before exiting The tasks of calculating the forces implementing the equations of motion and updating the dynamic variables each timestep are managed by function do step The conceptually simple calculation of forces and stepping of co ordinates is made lengthy and less transparent by the need to handle the rigid molecule equations of motion and the constant pressure and temperature extended system equations The procedure is set out in figures and The main stages of the calculation are a update all of the centre of mass co ordinates quaternions and extended system variables according to steps i and ii of equations 3 This is done by function step 1 which also applies the quaternion normalization and 46 input initialization rdf inner over timesteps link cell RDF calc istep do step 4 bin pair distances Step co ords by dt dump write trajectories every dump interval steps output periodic write rescale values log thermo values scale velocities averages calc print averages print rdf calc print RDFs write restart write backup file
104. h is much more efficient The intra molecular potential energy which is necessarily included in the site site interactions is computed separately and subtracted from the total to yield just the inter molecular part Intra molecular forces sum to zero in the calculation of the centre of mass forces and molecular torques just as with the reciprocal space contribution to the forces The internal stress is computed using a novel approach inspired by the extension of the Verlet neighbour list algo rithm developed by Bekker et 5 Bekker showed that the tensor virial may be written N N W Y n i l j i l N N rifi ng 2 25 i l j i 1 where n is the lattice vector between the central MD cell and its periodic images is also used as an index to represent a sum over all image cells and the expression j n refers to the periodic image of particle j in image cell n The quantity Sn eas Me N N 1 Y j i 1 i 1 j l fmi 2 26 is half the total force exerted by all particles in the central MD cell on all particles in image cell n This is easily evaluated in the inner loops of the link cell force calculation by summing separately the pair forces on all particles in each image cell The image cell translation vector n has already been computed in order to evaluate the site site vector rj nij Using equation 25 to compute the virial has another advantage It is not neces
105. h the elements of the distance vector 1 dU ri force calculation loop one for each kind of potential The potential type in use is passed as a parameter to kernel and is used in a switch statement to select the appropriate code To add a new potential the following modifications must be perfomed It calculates the corresponding values of which it stores in forceij There are several variants of the The array of structs called potspec must be extended The new array element should contain the name of the new potential in all lowercasd and the number of potential parameters for each site The name of a potential type given in system specification files will be matched against this name after being converted to lowercase n parallel with potspec the array dim must also be updated with a new entry which describes the dimensions of each parameter for the new potential This array is used by the input routines to convert from the input units into program units Each entry consists of triplets one triplet per parameter containing the powers of mass length and time corresponding to the dimensionality of that parameter e Define a new preprocessor symbol after the line define MCYPOT 2 to an integer which will be used as the index of the new type in the arrays and to select the code in the case statement The value must correspond to the index of the new entry in potspec starting from 0 in accordance with the usual
106. he unit cell dynamics 2 7 2V 1 pos ls s l m n p h h Tr o oh h ho h ho ch oh ho h 2 W p w V ty ht ooh oh ho 2 35 which rectifies these deficiencies The atomic dynamics are governed by equation as before The choice of fictitious mass W is governed by similar considerations to those concerning the heat bath mass Q discussed in section 2 6 2 Pl 2 7 Constrained h Matrix Nos and 1 87 describe and address the problem of the whole MD cell rotating during the course of the simulation This is because angular momentum is not conserved in a periodic system and because the h matrix has 9 degrees of freedom three more than needed to specify the position and orientation of the MD cell Their solution is to constrain the h matrix to be symmetric and involves a modification of the Parrinello Rahman equations Moldy incorporates a different constraint which is simpler to implement as it does not require modification of the equations of motion and which also has a more obvious physical interpretation The three sub diagonal elements of the h matrix are constrained to zero That is the MD cell vector a is constrained to lie along the x axis and b is constrained to lie in the xy plane Physically this may be thought of as implementing an MD box lying on a horizontal surface under the influence of a weak gravitational field The implementation is trivial at each timeste
107. ich contains information about the contents of the file followed by a number of dump records which contain the actual data Several control parameters govern dumping It starts at the timestep specified by begin dump and a dump record is written every dump interval timesteps thereafter After ndumps dump records have been written to a file it is closed and another is begun Filenames are generated from a prototype given by the parameter dump file by appending a number so that if the prototype is MDDUMP then successive files will be named MDDUMPO MDDUMP1 MDDUMP2 etc If it is not convenient for the sequence number to appear at the end of the file include the characters d at an appropriate point For example under VMS specifying dump file mddump d dat will name the files mddump0 dat mddump1 dat etc Each dump record is a sequence of single precision floating point binary numbers These are written either in native i e the machine s own format or XDR format see section depending on the value of the control parameter xdr The record s exact contents are determined by the control parameter dump level which is a bit flag i e a value of 2 means that bit is set Four bits are used and any combination may be specified but the cumulative values 1 3 7 and 15 are most useful A value of 0 disables dumping The data dumped for each bit are as follows bit 0 centre of mass co ordinates quaternions unit cell matrix and potential energy
108. ion are handled in exactly the usual manner However there are situations in which this atomic approach is impractical Because the system being modelled is essentially two phase the atoms of the framework find themselves under the influences of two distinct kinds of force There are the strong forces usually covalent or ionic which bind the atoms to form the framework itself and the weaker non bonded forces of the interaction with the fluid Ideally these could all be modelled by a single consistent set of interatomic potentials which are sufficiently transferable to yield an accurate crystal structure for the framework as well as the molecular structure of the fluid and its interaction with the surface Regrettably such potentials are hard to find Furthermore the characteristic vibrational frequencies of the solid framework will probably be much higher than those of the fluid Consequently the timestep must be chosen sufficiently small to accurately model the crystalline vibrations This will usually be far smaller than the value required for the fluid necessitating very lengthy MD runs to model both regimes properly This is of course exactly the argument used to justify rigid molecule models Moldy implements a variation on the rigid molecule model to simulate a rigid framework structure periodically extended throughout space There are a few subtleties which must be correctly handled to achieve a consistent imple mentation which are
109. ion is the specification of the potential functions This takes the form potential type i Pi se k l Pa Pu Pu m n phn cem Phm end where potential type is one of the keywords Lennard Jones Buckingham MCY generic morse or hiw to identify the kind of potentials to be used i j k m n are site ids and Pi is the potential parameter between sites i and j There should be one line for each distinct pair of site ids If any pair is omitted a warning is issued and the parameter values are set to zero The meaning of the parameters for the currently defined potentials 15 as follows Lennard Jones The potential is 6 rij 6 rij o rij and has two parameters which occur on each line in that order Note that the definition of includes the factor of 4 more usually separated out The control parameter time unit may be divided by four to read potentials specified in the standard form Buckingham This includes potentials of the Born Mayer type and has formula rij Aij r Bijexp Cijri The three parameters appear on each line in the order A B C 22 MCY This type supports potentials of the form of the water model of Matsuoka Clementi and 33 Aijexp Bijrij Cijexp Dijrij and the four parameters appear on the line in the order A B C D generic This contains a number of inverse power terms and an exponential r
110. is running on This is done using a global integer variable ithread which contains the index or rank of the processor Another variable nthreads contains the value of P These are set up in function par begin called from main and passed by external linkage to every module which needs them Processor 0 performs all input output which is arranged by testing ithread prior to any I O call e g if ithread 0 amp amp control istep control print interval 0 output which is the code in main used to print the regular output These variables are also present in the serial code where they are harmlessly redundant In that case ithread is always set to 0 and nthreads to 1 Parallel c also contains the function replicate to pass the system specification and dynamical variables from processor 0 to all of the others This is called from main after start up has returned It allocates memory for the system specification and dynamical variables on all other processors using the same functions as start up itself and calls sysdef and copy dynamics to broadcast the data from processor 0 to all of the others function par broadcast defined earlier in parallel c calls the underlying communications library to broadcast the data globally On exit from replicate every processor now contains a complete set of data and is ready to begin the run Error or warning conditions are handled by function message in o
111. l interactions between pairs of sites within the cutoff are included The neigh bouring cell list contains all pairs of cells with any parts closer than the cutoff radius This ensures that all appropriate interactions are included Furthermore all interactions between sites further apart than the cutoff are excluded by the expedient of setting their separation to a large value in the potential calculation This means that large and asymmetric molecules are handled correctly For a given cutoff radius the cell based mode is rather quicker than the strict mode since the neighbouring cell list is much smaller and fewer interactions are computed However to ensure that significant interactions are not omitted the cutoff ought to be set to a greater value than strictly required This tends to offset the potential speed gain On the other hand if strict isotropy is required in a liquid simulation for example then the strict cutoff option ought to be used 2 6 Temperature Initialization and Control From the equipartition theorem every degree of freedom in the system f has the same kinetic energy given by K 187 The effective temperature of the system is therefore given by the ensemble average of its kinetic energy 9 2 Y v 1 2 27 RATE jf Wj gi NA 7 S Here is the instantaneous kinetic energy of degree of freedom f g is the number of degrees of freedom N is the number of molecules J is an i
112. l site pairs separated by less than rdf limit are included This parameter may be varied independently of the interaction cutoff thereby allowing RDFs to be evaluated out to large distances without incurring the time penalty of increasing the cutoff 1 Every rdf out timesteps by default 5000 the RDFs are calculated from the binned distances and printed out and the counters are reset to zero to begin accumulation again Distances are binned and RDFs amp ap r calculated separately for each distinct type of atom atom or site site pair An explanation of the output format is given in figure Note that each number should be considered as the value at the centre of its bin so that entry i in each list is the value of gap i 1 2 b where b is the bin width There are a couple of tricks which may be played with the system specification if the atomic pair RDFs do not give exactly the functions required Firstly it is possible to calculate RDFs about a particular site distinguishing it from otherwise identical atoms by assigning it a new and unique site id in the system specification file This is the MD equivalent of the isotopic substitution method used in neutron diffraction Secondly if the molecular pair correlation is required this is identical to the RDF of an atom located at the molecular centre of mass A ghost site without charge mass or potentials may be added if necessary 3 7 Dumping The dump facility is provided in order to allow
113. lation Function from the same data Source For a species of N particles the Z t is defined as and the VTF as MN 20 t vu to t B 3 t 0n 1 1 N N VTF t Y qn Vy fo Y 4 t o 0 n 1 n l where 4 is the total charge of the species v t the velocity at time f The self diffusion coefficient of a species D can be obtained from the integral of its VACF 1 D Z t dt B 5 Similarly the frequency dependent ionic conductivity 6 be calculated directly from VTF via 1 l VTF te dr B 6 where V is the volume of the simulation box kg is Boltzmann s constant and T is temperature The command line options are mdvaf s system specification r restart file 1 dump file format t dump range v vaf range g species range c q output file where the arguments have the following meanings read the species data from a system specification file read the species data from a restart file Only one of s or r may be given if reading system specification file skip any preceding control parameter info read the centre of mass velocities from a dump file of dump level gt 2 given as a prototype name containing a printf format string see section 7 range of records in dump file specified in selector format time intervals to calculate the VACF or VTF values for specified in
114. lgorithms Mol Simul 8 1992 165 178 16 Optimization of the Ewald sum for large systems Mol Simul 13 1994 no 1 1 9 17 Message Passing Interface Forum MPI A message passing interface standard Int J Supercomput App 8 1994 no 3 4 165 18 Al Geist Adam Beguelin Jack Dongarra Weicheng Jiang Robert Manchek and Vaidy Sunderam PVM Parallel virtual machines a user s guide and tutorial for networked parallel computing MIT Press Cambridge MA 1994 19 H Goldstein Classical mechanics 2nd ed Addison Wesley Reading MA 1980 20 J P Hansen and I R McDonald Theory of simple liquids 2nd ed Academic Press London 1986 21 Robert J Harrison Tcgmsg send receive subroutines Battelle Pacific Northwest Laboratory 1994 Available by internet ftp from url ftp ftp tcg anl gov pub tcgmst tcgmsg 4 04 tar Z 66 22 D M Heyes Large numbers of particles on the Cray 1 Information Quarterly 5 25 1987 57 61 23 W Hoover Canonical dynamics equilibrium phase space distributions Phys Rev A 31 1985 1695 1697 24 Sun Microsystems Inc External data representation Protocol specification 1987 RFC1014 25 W L Jorgensen Revised TIPS model for simulations of liquid water and aqueous solutions J Chem Phys 77 1982 4156 63 26 E Keller Neues von SCHAKAL Chemie in unserer Zeit 20 1986 no 6 178 181 27 B W Kernighan and D Ritchie The C programmi
115. ly one accurate to O 8r in the co ordinates and O 8t in the velocities compared to the velocity Verlet algorithm which is only O 6t in the velocities This is insufficient for an accurate determination of the kinetic energy pressure and other dynamic quantities More seriously it fails badly in the case of polyatomic molecules and for constant pressure and constant temperature methods where the generalized velocities enter the dynamical equations themselves Velocity dependent forces occur in equations 2 4 and 2 10 in the Parrinello Rahman constant pressure equations section and in the Nos Hoover heat bath algorithms section D 6 7 These usually present a problem to non predictor corrector algorithms which are based on the assumption that the forces depend only on the co ordinates Fincham has devised a scheme to allow integration of the rotational equations using Verlet like algorithms II5 which is widely used despite the potential problems caused by the low accuracy of the velocities being propagated into the dynamics These cases are easily and accurately handled by the modification to Beeman s equations proposed by the These may be summarized using the symbol x to represent any dynamic variable centre of mass co ordinate quaternion or MD cell edge x and 9 to represent predicted and corrected velocities respectively i x8 x t 80 Fl AX r x r 81 ii t 0 x n
116. mory needed per processor increases with system size Mostly the scaling is linear but the reciprocal space sum uses temporary arrays whose size scales with the product of the number of sites and k vectors and hence to the 3 power of the system size An alternative version of ewald c which implements the reciprocal space term of the Ewald sum by distributing over sites rather than k vectors is included in the distribution as ewald RIL c It is based on RIL algorithm of Smith and distributes the temporary arrays containing the cos k r and sin k rj over the nodes In other words each node only stores the terms involving the r s to be considered on that node Since these arrays are by far the largest users of memory there is a substantial decrease in overall memory requirement Moreover the size per node now scales linearly with the number of k vectors and therefore assuming is optimized to the two thirds power of the number of sites These arrays will not therefore dominate the memory requirement in the limit of large numbers of processors and system size The disadvantage of the RIL scheme is that the partial sums of cos k r and sin k r must be summed over nodes separately for each k vector Though the amount of data transferred each time is small the communication and inter processor synchronization is far more frequent than for the RKL scheme and the parallelism becomes very fine grained 34 The upshot is that only machines
117. mpleted Because of this equation can only set the temperature to an accuracy of 1 1 This is often inadequate for purposes of comparison with experiment In order to allow the temperature to be set with greater accuracy Moldy allows the use of a partial average in the denominator _ 2 29 where K is the rolling average of K over some number of preceding timesteps That number is determined by control parameter roll interval This option should be used cautiously The change in K upon scaling only has a significant effect on the average after many timesteps If the subsequent scaling is performed before this change is reflected in the value of it will use an out of date value of the average kinetic energy It is therefore recommended that the number of timesteps between scalings be greater than or equal to the number used to compute the rolling average Otherwise it is possible to produce wild overshoots and oscillations in the temperature Finally there is a method for tackling really difficult cases when even individual scaling is unable to keep the temperature under control This might be a far from equilibrium configuration where the potentials are so strong that velocities rapidly become very large or when a single molecule acquires a very large velocity In that case the velocities can all be re initialized randomly from the Maxwell Boltzmann distribution periodically This provides a kind o
118. mply type compile to execute the command file compile com This will build Moldy and the utilities DOS MS Windows 3 There is a makefile for Borland Turbo C called Makefile mak in the standard distribution This must be edited to select the appropriate compiler options before executing make f Alternatively the programs may be built within the interactive environment Consult the makefile to find out which source files to link to build Moldy and the utilities Moldy has also been built using Watcom C Windows 95 98 NT The most straightforward way to build Moldy is to install one of the ports of the GNU gcc compiler such as Cygwin or Mingw32 Cygwin provides a full unix like environment while Mingw32 has the advantage that the executable files do not depend on a DLL file and are therefore more portable On the other hand only Cygwin has support for the XDR calls to produce portable restart and dump files see section B 3 2 As of version B 20 Cygwin incorporates Mingw32 which can be selected using a command line option Cygwin may be downloaded from the homepage at URL and links to Mingw32 as well as general information on GNU compilers under Windows may be found at It is also well worth installing the additional sunrpc ec es contains the XDR headers and libraries needed to produce portable restart and dump files see section B 3 2 A binary installation for Cygwin is available from URL ftp tp ranken tar gz There is also a so
119. mt defined in values c and array of classified types of data which contains pointers to av the numbers of components and to format strings and units conversion factors for output formatting and printing This is the compile time struct array info again defined in values c To compute averages of a different quantity new entry should be added to the array info and another corresponding enum type to enum n in defs h The storage retrieval and averages computation functions will then recognize the new type and may be called from within values and averages The array av itself is set up using the struct hack to allocate space for the rolling average data The structure type av mt contains as its final entry an array with one element roll 1 When storage for av is allocated the amount requested is calculated as if the member array was declared roll n with n equal to the rolling average interval The array member roll may then be used as if it had n members This does mean that the array av can not be simply indexed using pointer arithmetic or array subscripting since the compiler s notion of the size of a pointer of type av mt is wrong It is instead accessed solely via the pointers contained in array av info The general functions add average etc in values c provide the means to store and retrieve values and compute averages The function av ptr is provided to make the size and base address of the averages datab
120. nardus and H J C Berendsen An efficient box shape independent non bonded force and virial algorithm for molecular dynamics Mol Simul 14 1995 137 151 6 Berthaut L nergie lectrostatique de r seaux ioniques J de Physique Rad 13 1952 499 505 7 K C Bowler R D Kennaway G S Pawley and D Roweth An introduction to OCCAM 2 programming Chartwell Bratt 1987 8 K Cho and J D Joannopoulos Ergodicity and dynamical properties of constant temperature molecular dynamics Phys Rev A 45 1992 no 10 7089 7097 9 K Cho J D Joannopoulos and L Kleinman Constant temperature molecular dynamics with momentum conser vation Phys Rev E 47 1993 5 3145 3151 10 E Clementi G Corongiu and J H Detrich Parallelism in computations in quantum and statistical mechanics Comp Phys Commun 37 1985 287 294 11 Charles L Cleveland New equations of motion for molecular dynamics systems that change shape J Chem Phys 89 1988 no 8 4987 4993 12 W Deitz W Riede and K Heinzinger Molecular dynamics simulation of an aqueous MgCl solution Z Naturforsch 37a 1982 1038 1048 13 D J Evans On the representation of orientation space Mol Phys 34 1977 no 2 317 325 14 D J Evans and S Murad Singularity free algorithm for molecular dynamics simulation of rigid polyatomics Mol Phys 34 1977 no 2 327 331 15 D Fincham Leapfrog rotational a
121. nd may be lifted in future versions 17 2 10 2 Stress and Pressure Undefined There is no unique definition of the internal stress or pressure of a framework system The pressure of a system in a space with periodic boundary conditions is defined in terms of the molecular virial ET 32 26 Fi 2 53 i ljzi But the framework has no centre of mass and so the quantity R can not be defined site virial formulation of equation D T5 is of no assistance as the definition of the internal co ordinate pj involves the centre of mass co ordinate R Neither can one simply choose a convenient reference R Since the force exerted by the fluid acting on the framework is in general non zero the term Y Y RF i 2 54 clearly depends R The situation may also be viewed physically The definition of pressure of a system is the derivative of the free energy with respect to volume But with an infinite rigid framework the volume derivative can not be defined The useful quantity in this case is the partial pressure of the fluid Though not currently implemented this may be rectified in a future version of Moldy Finally we note that in the case of a 2D layer structure which is not rigid in the third dimension the perpendicular component if the stress tensor does have a physical meaning and is correctly calculated 2 10 3 Charged Frameworks A minor complication arises when using a framework which has a non zero n
122. nded to three dimensions and yields the results for the molecular and inter line spacings 4 and d respectively 1 VERE ven dy L 2 44 k 1 1 d Lr assuming h k is an integer The equal spacing requirements are satisfied approximately by how NM k m NIB 2 45 l 1 which when substituted into equation 2 44 yields dy d LN 1 5 2 46 15 Using this method an arbitrary number of molecules may be packed into a cubic cell with a guaranteed minimum centre of mass separation of approximately LN 1 5 n contrast to other methods such as random placement with excluded volume it will always yield a configuration no matter how high the density It is also very simple to implement It still remains to determine the initial orientations of the molecules in the case of polyatomics In the current im plementation these are assigned randomly which works well for small or near spherical molecules Further refinements which may help avoid overlaps for elongated molecules are possible such as a periodic sequence of orientations along the line but no investigation of this possibility has yet been carried out In practice the method has proved to work well for water and aqueous systems and always yields a configuration which may be evolved towards equilibrium by the MD equations of motion Any feedback on its performance in more difficult cases will be welcome 2 9 22 Initial Velocities
123. ng language 1st ed Prentice Hall Cambridge 1978 28 The C programming language second ed Prentice Hall Cambridge 1988 29 H Kistenmacher H Popkie and E Clementi Study of the structure of molecular complexes III Energy surface of a water molecule in the field of a fluorine or chlorine anion J Chem Phys 58 1973 no 4 5842 30 Donald E Knuth Fundamental algorithms second ed The Art of Computer Programming vol 1 section 1 2 Addison Wesley 1973 31 A Laakonsen Computer simulation package for liquids and solids with polar interactions I McMOLDYN H2O aqueous systems Internal Report KGN 41 IBM Corporation Kingston N Y 12401 USA 1985 32 S W De Leeuw J W Perram and E R Smith Simulation of electrostatic systems in periodic boundary conditions I Lattice sums and dielectric constant Proc Roy Soc A373 1980 27 56 33 O Matsuoka E Clementi and M Yoshimine CI study of the water dimer potential surface J Chem Phys 64 1976 no 4 1351 1361 34 Richard Miller A library for bulk synchronous parallel computing Parallel Processing Specialist Group workshop on General Purpose Parallel Computing British Computer Society 1993 35 S Nos A molecular dynamics method for simulations in the canonical ensemble Mol Phys 52 1984 255 268 36 Molecular dynamics simulations at constant temperature and pressure Computer Simulation in Materials Science
124. nitial density in gt Used by skew start only to determine initial MD cell dimensions scale interval integer 10 Number of steps between velocity scalings scale end integer 1000000 When to stop scaling const temp integer 0 1 for Nos Hoover 2 for Gaussian thermostat ttmass real 100 Translational inertia parameter for Nos Hoover thermostat kJ mol ps rtmass real 100 Rotational inertia parameter for Nos Hoover thermostat kJ mol ps scale options integer 0 Select variations on scaling or thermostat temperature real 0 Temperature of initial configuration for scaling and thermostat K const pressure integer 0 Parrinello and Rahman constant stress 2 P R compatible constant pressure ensemble 3 Wentzcovitch Cleveland constant stress 4 Andersen constant pressure w real 100 0 Extended system mass parameter amu pressure real 0 External applied pressure MPa strain mask integer 200 Bitmask controlling A matrix constraint alpha real auto parameter for Ewald sum Set any negative value to disable Ewald sum k cutoff real auto Reciprocal space cut off distance in f cutoff real auto Direct space cutoff distance in strict cutoff boolean false Flag to select rigorous or cheap but approximate cutoff algorithm surface dipole boolean false Include surface dipole term in Ewald sum roll interval integer 10 Period over which to calculate rolling averages print interval integer 10 How frequently to p
125. nother In that case the dipole moment of the MD cell changes discontinuously Because of the surface dipole term the calculation would model a discontinuous macroscopic change in the dipole moment of the whole system caused by an infinite number of ions jumping an infinite distance This is manifested in practice by a large and discontinuous change in the energy of the system and on the force on each charge within it This situation is completely non physical but is easily avoided However the problem may also arise more subtly even when there are no mobile ions if a framework is being simulated section 2 10 The framework is treated as a set of discrete but fixed atoms rather than a molecular unit If the shape of the unit cell is allowed to vary then ions constituting the framework may indeed cross MD cell boundaries causing the aforementioned problems 2 4 4 Stress The internal stress and pressure of an atomic system is given by the volume derivative of the internal energy The situation is slightly more complicated for rigid molecules since molecules do not scale with volume and only the inter molecular distances vary The resulting expression for the Coulombic part of the instantaneous stress v 15 37 Appendix A XY cement vip Enn CON C A Inj n Real space term 2 N 2 _ 1 1 N pe Y 5 2 Ez kk 249 kong _
126. ns exceed the limits of the integration algorithm and timestep which results in Moldy stopping with a quaternion normalization or quaternion constraint error Judicious use of this option every few timesteps using scale interval ought to allow the system to relax to a state sufficiently close to equilibrium for normal scaling to take over Bit 2 is intended to deal with the problem of setting the temperature accurately using scaling The ensemble average kinetic energy which characterizes the temperature of the system and the instantaneous value fluctuates about this value This may be necessary if the restart file is located on a different device or disk partition from the current directory To rename the temporary file successfully it must reside in the same partition or device as the restart file 8From version 2 1 onwards Because the XDR calls are not part of ANSI standard C however the XDR code is conditionally compiled into Moldy only if the USE XDR preprocessor symbol is defined during compilation 25 Nov 17 15 13 06 1989 Water test Trans KE Rot KE Pot Energy Tot Energy TTemp Stress Stress Stress 243 88 453 88 187 35 533 5 305 5 22 053 0 1 0401 221 0 589 464 120 464 373 901 120 901 207 Rolling averages over last 10 timesteps 240 27 319 31 82 472 533 39 301 0 400 0 342 9 12 53 0 00 0 00 1 20403 296 127 22 077 0 34 205 221 3 0 0 0 00 12 53 0 00 296 589 133 0 00 000 12 53 127 133 132 Standard
127. nstantaneous temperature 10 It is almost always desirable that a simulation be conducted so that the temperature is the supplied parameter rather than the kinetic energy This requires some mechanism to fix the average kinetic energy at thermal equilibrium The initial kinetic energy may be set approximately by choosing random velocities which sample the Maxwell Boltzmann distribution at the desired temperature and this is indeed what Moldy does on starting a new run see section 2 9 2 But because the initial configuration is usually far from equilibrium it will have too much potential energy As the run progresses this will be converted into kinetic energy raising the temperature above the desired value It is therefore necessary to have some mechanism for removing excess kinetic energy as the run progresses Moldy offers several mechanisms to control the temperature The common technique of velocity scaling is suitable for use during the equilibration period but does not generate meaningful particle trajectories The Nos Hoover method couples the system to a heat bath using a fictional dynamical variable and the Gaussian thermostat replaces the Newton Euler equations by variants of which the kinetic energy is a conserved quantity 2 6 1 Rescaling At periodic intervals linear and angular velocities are multiplied by a factor of kpT 2 28 where is the desired temperature By repeatedly setting the instantaneous temper
128. ntil the 47 do step Calculation of molecular centre of mass forces and torques mol force Get C of M forces Calc distant pot l amp pressure correction any mol_torque polyatomics Get torques Allocate arrays for sites and site forces Get atomic molecular virial Set up array of correction term site charges Store from a last timestep in step 1 old arrays etc Apply Beeman steps i ii to update dynamic variables newton New C of M acc s b Calculation of atomic site forces and potential energy d make sites Calc site co ords energy dyad kinetic stress term rahman Get unit cell acc s beeman 2 Step cell velocities force calc real space forces Y ewald 9 coulomb forces recip space forces Get surf dipole surface dipole forces amp energy switch on NL Figure 5 4 a Flow diagram of function do step which performs a single timestep continued in Figure 48 Nose Hoover nhtherm Get N H Const stress or thermostat Gaussian beeman 2 Predict N H G gtherm Calc gaussian parinello Get P amp R acc corr Const stress simulation hoo
129. olecule i is Fi Y fap 241 jipa and the torque is given by Ni ria Ri X fi 2 2 Qa where 1 M is the centre of mass of molecule i The motion is governed by the Newton Euler equations MR 2 3 0 x Ij 0 2 4 where is the angular velocity of the molecule mio 2 1 iS the inertia tensor and pj Fia is the atomic site co ordinate relative to the molecular centre of mass The orientations of the molecules are represented by quaternions as has now become common practice They are pre ferred over Euler angles for two reasons Firstly they lead to equations of motion which are free of singularities 3 which means that no special case code is required This leads to much improved numerical stability of the simulation IA Secondly molecular symmetry operations and combinations of rotations are elegantly expressed in terms of a simple quaternion 1 14 42 A quaternion is an ordered number quartet which obeys the algebra given by Pawley 40 The multiplication rule in that reference may be restated as a matrix product treating each quaternion as a 4 vector If p po p1 p2 and comment on notation is appropriate here In this chapter site quantities are denoted by lowercase letters molecular quantities by uppercase sites are indexed by greek letters and molecules by roman A missing index denotes a sum over the corresp
130. ols for manipulating DCD files and converting to other formats are also available from http www ks uiuc edu Development MDTools In interactive mode you are prompted for the source of the system specification sys spec or restart file and the configurational information restart or dump file But you must either redirect the output using gt on unix or specify an output file with o or the output will go to the screen Alternatively the command line options are mdshak s system specification r restart file d dump file format t dump range 7c x v h p g b y extra text o output file where the meanings of the options are s read a system specification file r read a restart file Only one of s or r may be given c ifreading a system specification file skip any preceding control parameter info d read the configurational information from a dump file The argument is a dump prototype name containing a printf format string see section B7 t range of records in dump file specified in selector format i this inserts its argument into the output file and is used to add extra SCHAKAL commands such as BOX x Write an output file in XYZ format suitable for Xmol or RasMol v Write co ordinates from multiple dump files to an output file in DCD format suitable for animation by h Write an output file in SCHAKAL format p Write an output file in Brookhaven PDB format suitable fo
131. onding sites or molecules so that for example is a site site vector and F the molecule molecule force q qo 1 42 d3 quaternions then p P2 P3 40 P2 41 2 5 P2 P3 po pi q2 24 2 43 The quaternion d conjugate to q is defined as d qo 41 q2 43 so that qo 44 43 45 0 0 0 Q 6 The norm is defined as q g5 97 d 43 and q is called a unit quaternion if q 1 Any possible rotation can be represented by a unit quaternion Du Val shows 54 that if cos 2 lsin 2 where we have combined the last three components to form a 3 vector and p 0 r then the operation p 0 7 2 7 corresponds to a rotation of the vector by an angle of about the axis 1 components may also be expressed terms of the Euler angles asf 40 cos nd cos 41 sin 9 V sin v 8 q2 cos 2 x sin 2 43 sin 2 8 The relationship between the time derivative of a quaternion and the principal frame angular velocity was given by Evans II3 Equation 27 and rewritten using quaternion algebra by Refson 4 amp as 24 q 0 9 2 9 The second derivative is given by 24 a 1 2 4 214 2 2 10 Equations 2 10 and 2 4 allow the simulation to be implemented using quaternions and their derivatives as the dynamic variables for rotational motion and this is the me
132. onging to an ordinary molecule are placed in the cell which contains the molecular centre of mass With these modifications Moldy is able to successfully model a fluid in contact with a 3D rigid framework 2 dimensional systems such as a fluid at a surface or between two surfaces may be represented as a 3D system by adding an artificial periodicity in the third dimension To reduce the errors so introduced the framework can be made large with space inside to fill MD cell with a large c axis In the case of a 3D framework it is clearly not sensible to allow the MD cell to change shape or size since this would destroy the internal structure of the framework However if the framework represents a 2D layer or surface then the layer spacing may be allowed to vary using the constant stress equations section 7 and applying constraints to allow only the c axis to fluctuate In that case bear in mind that the dynamics are governed by the Parrinello Rahman equations of motion for the cell rather than the Newtonian equations for the layer In particular the mass is given by the parameter W rather than the mass of the framework molecule These may of course be set equal if required Note also that no layer layer interactions are calculated and any force is the result of the externally applied pressure f Finally there are two subtle complications which arise from the framework concept This is a restriction of the current implementation a
133. orces Simulations may be performed either in the usual NVE ensemble or in NVT NOH NoT ensembles using Nos Hoover thermostat and Parrinello and Rahman constant stress methods As the MD cell need not be cubic the program is equally suitable for simulations of solids and liquids Most existing MD programs are limited in their capabilities for example to one kind of potential function or molec ular symmetry or to some restricted number of molecules Moldy is as far as possible free from such arbitrary con straints The system is specified at the beginning of each run and its size is only limited by the amount of memory available to the program if a system is too large to handle the solution is to buy some more memory The system may contain a mixture of an arbitrary number of molecular species each with an arbitrary number of atoms and an arbitrary number of molecules of each Molecules or ions may be monatomic or polyatomic linear or three dimensional in any combination The potential functions may be of the Lennard Jones Buckingham including Born Mayer or MCY types and other potential types may be easily added Such flexibility is possible because Moldy is written in the language which permits dynamic memory allocation Moldy is written to be highly portable and has been tested on a wide range of computers and operating systems including VAX VMS MS DOS and Unix both BSD and system V varieties It should be straightforward
134. p the acceleration of those components j 18 set to zero which is equivalent to adding a fictitious opposing force This constraint technique is not restricted merely to eliminating redundant degrees of freedom but can also be used for other purposes For example it may be used to allow uniaxial expansion only MD cell constraints are selected using the control parameter strain mask see section 8 A value of 238 will freeze all of the off diagonal components of h 2 7 2 Constant Pressure The above h matrix constraints not adequate for simulations of liquids at constant pressure Since a liquid can not support shear stress there is no restoring force to keep the simulation cell nearly cubic and it will therefore drift to a parallelepiped or elongated shape Moldy also supports an alternative constraint which only permits isotropic MD cell expansion The isotropic variant of equation 2 35 gives 2 zd 2W V Wh Tr n p h 5 9 2 36 5Because of the different forms of the fictitious kinetic energy term Lagrangians typical values of W for equation are smaller by a factor of the order of L where L is the MD cell edge than those useful when the dynamics is given by equation 13 and if the cell is cubic this is equivalent to the constant pressure dynamics of Andersen Bl with W 3M Andersen Applying a uniform dilation constraint to the Parrinello Rahman Lagrangian gives equations of motion Tr X
135. particles in cell icell The iterations of this loop are distributed over processors by rewriting it as for icell ithread icell nx ny nz icell nthreads Evaluate forces on particles in cell icell where the integer variables ithread and nthreads contain the rank of the processor and the total number of processors respectively Since ithread has a different value on each processor the separate copies of the force arrays contain only the contributions from the set of subcells assigned to that processor After completion of the forces calculation the individual arrays are summed in function do step Similar partial contributions to the potential energy and the stress are also globally summed Provided that the number of subcells is rather greater than the number of processors and the number of particles per subcell is fairly constant this algorithm automatically ensures that the computation is divided evenly between the processors The reciprocal space sum may be parallelized either by distributing the loop over k vectors over processors L0 54 or by partitioning the particles between processors Smith 51 called the former strategy the reduced k vector list RKL method and the latter the reduced ion list RIL method Moldy contains implementations of both methods In the method it is necessary to pre compute a list of k vectors within the cutoff sphere see figure D 6 which is stored in the
136. pe name containing a printf format string see section BZ range of records in the dump file specified in selector format to average positions over bond limits in tenths of Angstroms entered in form minimum maximum Limits be entered as integers only so smallest difference is 0 1 Defaults to 0 2 to 2A angle limits in degrees entered in form minimum maximum Limits may be entered as integers only so smallest difference is 1 Defaults to 0 to 180 range of species to be included given in selector format where 0 represents the first species Defaults to all species if reading system specification file skip any preceding control parameter info name of optional output file Defaults to standard output Species positions can be input using system specification restart or dump files If the latter case is chosen the bond and angle data for each time slice are written consecutively to the same output file with each set of data separated by a title indicating the number of the respective time slice Bdist converts the data output from mdbond into a format suitable for plotting bond and angle distribution curves The utility scans through the output file counting the number of bond lengths angles within the specified intervals for 63 different pairs triplets of species types The data is output in two columns first for bonds and then for angles with the mid point of the bond angle interval given in the fi
137. perature as given by the control parameter temperature This is done for both starting methods 3 3 Restarting from a Previous Run At the end of a simulation run it is often desirable to store the configuration of the system in a file This restart file may be used at a later date to continue the simulation from that point rather than from scratch To instruct Moldy to write a restart file simply set the control parameter save file to a suitable filename to start from a restart file set restart file to be the name of that file Each restart file is a binary file which contains enough information to reconstruct the exact state of the system and of the program It includes a copy of all the control parameters in force the current timestep number a complete system specification all the simulation dynamic variables and the intermediate data used in the calculation of averages and radial distribution functions Thus a run continued from a restart file will proceed just as if there had been no interruption and will generate identical results provided the control parameters are not changed When continuing a simulation it is only necessary to explicitly specify control parameters which are to be changed Their previous values are read from the restart file and are used as defaults when reading the control file Consequently control files for restarting tend to be rather short Caution always include a new possibly null value for save file Othe
138. r VM D Xmol or RasMol g Wirite an output file in Cambridge CSSR format suitable for reading by e g Cerius2 b Write the atom co ordinates in raw binary unformatted format y Join up trajectory fragments to eliminate jumps when a molecule exits from one side of the MD cell and re enters the opposite side Mdshak may be useful for more than just visualization purposes as it extracts atomic co ordinates from the system specification positions and quaternions It is written in a modular fashion to make it easy to add a new output routine in file molout c usual FORTRAN sequential unformatted read is not suitable as this expects record size information to be present Dumpext can not write this as its format is entirely determined by the whim of the compiler writer FORTRAN unformatted sequential files are not guaranteed to be portable even between different compilers on the same machine Be aware that the FORTRANT7 standard does not guarantee what units the record length is specified in Predictably some manufacturers compilers use bytes and others words Consult the documentation 3 Available from http www ks uiuc edu Research vmd http www iavsc org 60 B 6 Msd Msd is a utility for calculating the mean square displacements of selected species from the dump files written during the course of a Moldy run A plot of msd against time enables the diffusion coefficient of that species to be calculated from the gra
139. r files The ANSI header files string h stdlib h stddef h and time h are missing from Berkeley unix or incomplete Replacements are included which may be dispensed with on an ANSI conformant system If the symbol ANSI LIBS is defined they simply include the system version Replacements for ANSI functions The ANSI function to delete a file remove the signalling function raise and the string functions strstr and strerror are missing from pre ANSI libraries Replacements are supplied in ansi c Replacements are provided in ansi c for functions memset memcpy and strchr which are missing from Berkeley UNIX e The function vprintf is often absent from older libraries Replacements are provided which call the internal function doprnt or b implements a portable vprintf Use the preprocessor macros HAVE VPRINTF or HAVE DOPRNT to select which Timing routines The supplied clock function on some 32 bit UNIX systems resets to zero after 36 minutes Re placements called cpu for system V and Berkeley UNIXes and POSIX are supplied in auxil c The func tion rt clock is also defined and returns the elapsed time in seconds For a non unix system cpu and rt clock simply call the ANSI functions clock and time File manipulation routines Auxil c contains the functions replace and purge replace renames a file making backup of any existing file of that name purge removes the previous or backup
140. re converted to internal units a m u ps and v a m u 2 4m 9 The prototype molecule for each species is then shifted so that its zero of coordinates lies on its centre of mass and rotated into the principal frame polyatomics only 3 2 2 Initial Configuration Moldy provides two methods of setting up an initial configuration By default the skew start method of section 2 9 1 15 used to place the molecular centres of mass in a regular arrangement which ensures molecular separation If there is more than one species present molecules of each are chosen randomly for each site Molecular orientations are chosen randomly from a uniform distribution This method has been found to work well for reasonably small or fairly isotropic molecules and it is anticipated that it will be the usual method of starting a simulation of the liquid state On the other hand if the constituent molecules are sufficiently large and irregular or if it is intended to simulate the solid state then the lattice start method will be more appropriate This method is activated by setting the control parameter latt ice start to 1 and creates the initial configuration by periodic replication of some crystalline unit cell In that case Moldy expects to find following the end which terminates the system specification an initial configuration specification of the following form a bc B y n n n species name Yi Zi dio qu 92 913 species name
141. reedom that is J 0 for some i the corresponding angular velocities etc are simply set to zero 2 10 Frameworks In addition to the bulk properties of solids and liquids much attention has been devoted in recent years to using MD simulation to model atoms or molecules interacting with surfaces or other structures such as zeolites The distinctive feature of such a situation is that the 2 or 3 dimensional surface or structure is larger than the interacting molecules by many orders of magnitude In fact the idealization of this system makes the structure infinite in extent atomistic model of this kind of a system necessitates choosing the periodic boundary conditions of the MD cell to be commensurate with the periodicity of the surface or structure This manual will refer to an infinitely repeating crystalline surface or structure as a framework 16 There two possible formulations of a system of this kind for use in a MD simulation Firstly the framework may be modelled as a set of independent atoms interacting via pair potentials This merely requires specifying the correct initial structure and MD cell plus a suitable set of pair potentials The latter model both the internal forces which determine the crystal structure of the framework and its interaction with the atoms or molecules of the fluid Conceptually there is no distinction between this situation and a regular solid or liquid system and the mechanics of initiating a simulat
142. rent accuracy is desired the cutoffs may be adjusted using equations p 23 The a parameter is specified by alpha in units of A and the direct and reciprocal space cutoff distances and ke by cutoff and k cutoff in units of and respectively The value of should only be changed after careful timing tests if the system size is large The power law given in equation gives a theoretical scaling of execution time with number of ions of T N15 In practise T 57 has been achieved over a range of from 72 to 7800 which is very close to optimal Important note The automatically determined value of is chosen to converge the Coulombic part of the potential only Due to the very general nature of the potentials it is not possible to choose automatically so as to guarantee convergence of the non electrostatic part Although in many cases the automatic value will be adequate it is the user s 29 responsibility to ensure that it is large enough If there are no charges the system specification file then re is not set and an error message is issued As an example of manual determination of the parameters for a simulation of 512 MCY water molecules the values a 0 3 re 9 and ke 1 9 give potential energies correct to approximately 1 part in 105 For a simulation including ions 1 1 Molal Magnesium Chloride solution the same accuracy is attained with a 0 45 re 9 and ke 3 7 The Ewald sum may
143. rewritten more usefully for an arbitrary multi component mixture by eliminating the molecular density and number N as gag r V 8 r ria rog 2 41 In the simulation this is evaluated by an expression very similar to equation 7 39 3Nnyis b _ 4npNaNg1 r 8r r 6 2 2 42 2 9 Initial Configuration One of the more trying aspects of initiating a molecular dynamics simulation is getting it started in the first place It is not hard to see the reason why An MD integration algorithm will only generate a good approximation to the correct equations of motion if the forces and velocities of the particles are less than some value The timestep is chosen so that this criterion is satisfied for near equilibrium configurations But if the configuration is far from equilibrium certain forces may be extremely large due to atoms approaching each other too closely And worse the breakdown of the integrator leads to breakdown of the conservation laws and the system evolves to a state even further from equilibrium One way around this difficulty is to start the system from a configuration known to be near equilibrium such as a known crystal structure For a solid state simulation this is the method of choice and Moldy allows any initial structure to be specified and replicated to fill the MD cell In the case of a liquid a crystalline initial state is less desirable and indeed none may be known Furthermore
144. rint normal output begin average integer 1001 When to start accumulating the thermodynamic averages average interval integer 5000 How frequently to calculate and print averages reset averages boolean false Discard accumulated averages in restart file begin rdf integer 1000000 When to start accumulating radial distribution function information rdf interval integer 20 How frequently binning calculation is performed rdf out integer 5000 How frequently to calculate and print RDFs rdf limit real 10 Calculate RDFs up to what distance nbins integer 100 Number of binning intervals between 0 and 21 rdf limit Table 3 1 Control Parameters continued name type default function xdr boolean true Write restart backup and dump files in portable binary format using Sun XDR dump file character null Template of file names used for data dumps begin dump integer 1 Timestep to begin dumping at dump interval integer 20 How frequently to perform dumps dump level integer 0 Amount of information to include in dump ndumps integer 250 Number of dump records in each dump file backup interval integer 500 Frequency to write backup file backup file character MDBACKUP Name of backup file temp file character MDTEMPX Name of temporary file used for writing restart configurations subcell real 0 Size of sub cell in to divide MD cell into for link cell force c
145. rmined after the input files have been partially read in but they must be allocated before the read is complete Thus reading of input files and array allocation must be interspersed of the general nature of the systems Moldy is able to simulate Many structures that might otherwise be hard coded must here be read from the input files and set up dynamically Moldy has three different modes of start up an initial run a restart and a backup restart The use of the stored values of the control parameters as defaults on a restart run is an added complication the ability to specify the input units the potential parameters are expressed in requires the appropriate conversion to be performed at start up Initialization is controlled by function start up whichis called directly from main It calls subsidiary functions to read the control and system specification restart or backup files It controls the allocation of dynamic memory for the system description and potential and dynamic variable arrays and the RDF and averages databases It oversees the 45 computation of quantities derived from the input parameters and system specification and generally transforms the input data from a form useful for human preparation to one useful for machine computation Figures and show a block diagram of start up and a brief description of the functions it calls Parameters read from the control file by read control which assigns their valu
146. rst column and the number of bonds angles within that interval in the second column The utility is invoked on the command line by bdist i input file b bond range a angle range o output file where the options are i name ofthe input file b bond interval limits in tenths of Angstroms entered in form minimum maximum increment Limits may be entered as integers only so smallest increment is 0 1A Defaults to 0 2 to 2A with a 0 5A increment a angle interval limits in degrees entered in form minimum maximum increment Limits may be entered as integers only so smallest increment is 1 Defaults to 0 to 180 with a 1 increment Note that if the mdbond output file to be analyzed contains multiple data sets from reading dump files the data for each time slice must be copied manually into separate files and the Time Slice headings removed before bdist can be used B 10 Ransub It is sometimes necessary to simulate a crystalline system in which one or more of the components have been randomly substituted onto lattice sites normally occupied by another species For large numbers of substitutions doing this by hand can be tedious so the utility ransub has been provided to speed up the construction of system specification files in which solute or dopant molecules are randomly distributed The command line options are ransub s system specification r restart file 1 dump file format t d
147. rwise when the new run terminates the new restart file may overwrite the old one Neither is it necessary to repeat the system specification since that too is stored in the restart file However there are occasions when it is desirable to do just that for example if the value of one of the potential parameters is to be modified In that case set the switch new sys spec to 1 true and provide a system specification as per a new simulation This is checked for consistency with the existing one and if correct replaces it The following checks are applied which only verify that it is sensible to assign the old dynamic variables to the new system 1 The number of species must be the same 2 Each species must have the same number of rotational degrees of freedom as its predecessor It is not possible to replace a polyatomic by a monatomic or linear molecule for example 3 The number of molecules of each species must not change This means that the order in the specification file must be identical too It is however possible to change the number of sites on a molecule subject to 2 3 3 1 Periodic Backup Closely related to restarting is the backup mechanism This is provided to guard against the complete loss of a simulation due to computer failure Periodically during a run Moldy writes its state to a backup file which is in fact just a restart file In the event of a crash the simulation can be restarted from the point the last backup was written r
148. s a single function argument avoiding long and cumbersome argument lists and the consequent risk of programming error All of the major struct types used in this way are defined in the header file structs h where each type is carefully annotated with the meaning of its members The struct control is shared between almost all modules in the program using external linkage It contains the values of the parameters from the control file the exact mapping is defined by the array of structs match declared in source file startup c The values are read from the control file by function read control in source file input c and adjusted by start up in startup c where the timestep related parameters are updated and the floating point values are transformed into program units These are the only functions which alter values in cont rol Most of the information needed to describe the system is stored in the struct system and the array of structs species Struct system contains counters to record the total numbers of species nspecies atomic sites nsites molecules nmols polyatomic molecules nmols_r and pointers to the arrays of dynamical variables Its counterpart species is a dynamically allocated array length nspecies of structs of type spec mt which contains individual data for each molecular or atomic species This includes the mass moments of inertia dipole moment and charge the number of sites belonging to this species and the number o
149. s quaternion normalization or constraints see section if the program runs out of memory or if a restart file can not be correctly opened or is of the wrong format Most of the messages are self explanatory However there are two fatal errors which occasionally arise and are somewhat cryptic F Quaternion n normalization error in beeman and F Quaternion n constraint error x Technically these refer to violations of the conditions that the quaternions representing the angular co ordinates be normalized to 1 and that equation E TT be satisfied Either may occur if the angular velocity of some molecule becomes too high for accurate integration of the equations of motion This may have a number of causes First the timestep may simply be too large Second the system may be in a state where atoms are so close as to generate large torques which accelerate the molecule to a high angular velocity This commonly arises if the starting configuration is very far from equilibrium particularly in the case of molecules with small moments of inertia such as methane In most cases the simulation may be restarted using strong rescaling or a Gaussian thermostat to limit velocities during the approach to equilibrium Occasionally a smaller timestep may help during equilibration The third cause of normalization or constraint errors is an error in the potentials or the units which allows for a state with high or infinite negative binding ener
150. s will vary on different processors and for different potentials its value is not critical since it enters the equations as a sixth root It must be emphasized that the is used as a cutoff for the short ranged potentials as well as for the electrostatic part The value chosen above does not take the nature of the non electrostatic part of the potential into account It is therefore the responsibility of the user to ensure that is adequate for this part too 2 4 2 Uniform Sheet Correction The Sth term in equation 19 is necessary only if the system has a non zero net electric charge and is useful in special cases such as framework systems In a periodic system the electrostatic energy is finite only if the total electric charge of the MD cell is zero The reciprocal space sum in equation for k 0 takes the form 2 4a i 1 which is zero the case of electroneutrality but infinite otherwise Its omission from the sum in equation 2 19 is physically equivalent to adding a uniform jelly of charge which exactly neutralizes the unbalanced point charges But though the form of the reciprocal space sum is unaffected by the uniform charge jelly the real space sum is not The real space part of the interaction of the jelly with each point charge as well as the self energy of the jelly itself must be included giving the fifth term in equation D T9 2 4 3 Surface Dipole Term The optional final term in equations 2 19 and
151. sary to perform the sum of contri butions within the inner loop of the force calculation c f equation 2 15 since it depends only on the site co ordinates and the total force acting on each site This gives a considerable reduction in the number of computations in the costliest part of the whole calculation 2 5 2 No minimum image convention One notable feature of the implementation of the link cell method in Moldy is that it does not follow the minimum image convention used in most MD codes Instead the list of neighbouring cells and hence interactions considered includes all periodic images of a particle which are within the cutoff This means that it is quite safe to use a cutoff of more than half of the MD cell side in any direction 2 5 3 Cellor Strict cutoff There are two options for the way in which site site interactions are selected for inclusion in the total energy and forces These are cell based cutoff and strict cutoff and are selected by the strict cutoff control parameter In cell based mode strict cutoff 0 the neighbour cell list is built to include only those cells whose centre is within the cutoff radius of the centre of the reference cell interactions between sites in the neighbouring cell list are computed This is a quick and dirty method as some interactions between sites closer than the cutoff will inevitably be excluded whereas some outside the cutoff range will be included In strict mode strict cutoff 1 al
152. st cases it should be possible to do configure make oo The configure script will determine the capabilities of the compiling system and create the files config h and Makefile tailored to the system It has built in information about the compiler names and optimization options it should recognise most major workstation types It should be able to determine the characteristics of an machine unknown to it and to create Makefile and config h to build a working version In that case you may need to hand tune compiler or optimization options which may be specified by setting the OPT and OPT2 environment variables Since these date rapidly as vendors improve their compilers and introduce new architectures feedback from users on appropriate options is welcome and will be included in future releases 32 You set or change the default compiler or options by calling configure with the environment variables CC CPPFLAGS CFLAGS OPT OPT2 LDFLAGS or LIBS defined Any values you supply for these will override the built in ones e g env 01 2 04 configure c shell 5 01 configure bourne shell korn shell bash OPT2 is used to compile only the most performance critical modules and will usually select a very high level of optimization It should be safe to select an optimization which means treat all function arguments as restricted pointers which are not aliased to any other o
153. such a starting configuration restricts the allowed number of molecules to a multiple of the number in the unit cell and worse may force the use of a non cubic MD cell Moldy incorporates a novel method of generating an initial configuration which in the main overcomes these prob lems 14 og oog te g oo g oo Figure 2 1 The Skew Start method N molecules are placed at evenly spaced intervals a a line joining a corner of the MD cell to its image cells away 5 in this illustration When mapped back into the original cell this guarantees a minimum separation of min d a 2 91 The Skew Start Method The essence of the Skew Start method is to generate a configuration which though not periodic in 3 dimensions nev ertheless is sufficiently regular to guarantee a minimum separation between molecular centres of mass Figure 2 demonstrates the principle in 2 dimensions The N molecules are placed at some suitable interval a on a line drawn between one corner of the MD cell of side L and one of its periodic images Clearly the index of the image cell corner should be chosen so that the molecule spacing a is close to the spacing of the line s images d For simplicity choose k 1 which leads to the condition LVh 1 L ax re N n 1 gt h m VN 1 2 43 Therefore the optimum value of h is the nearest integer to VN 1 The formalism may be exte
154. sure the virial mean square forces mean square torques and total dipole moments are calculated 3 51 Output units All of the various forms of output use the same units though for brevity they are not explicitly mentioned on the periodic output Lengths are measured in A energies are all expressed in kJ mol temperatures in Kelvin pressure and stress in MPa mean square forces and torques in N mol and Nm mol charge in electron charges and dipole moments in Debye Because energy is an extensive quantity the printed values refer to the whole system There is no practical way of expressing energies per mole of any particular constituent in a program capable of simulating arbitrary mixtures If these units do not suit they can be changed in the configuration file defs h where the conversion from internal to output units is parameterized Remember that the standard deviation is a measure of the fluctuations about the mean not the uncertainty in the mean For that the standard error in the mean is required which is more difficult to evaluate Theoretically it is the s d divided by VN where N is the number of independent observations But successive timesteps are highly correlated and do not count as independent See ref 0 section 6 4 page 191 onwards for a full discussion 26 Radial Distribution Functions Bin width 0 1 0 0 RDF 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000000 0 000
155. t re the interaction is computed This method suffers from several disadvantages Since for any given site the interaction with any other site is considered only once only the nearest periodic image of that site is included the minimum image convention However this restricts the cutoff radius to less than half the MD cell dimension r 2L More serious is the way the computational time scales with the number of sites If there are N sites there are O N separations to compute and the overall time therefore scales as O N The Verlet Neighbour List scheme 2 pp 146 149 makes use of the spatial locality of the interaction potential by maintaining a list for each site of all the neighbouring sites i e all those within the cutoff distance This can give considerable gains for moderate numbers of sites but it ultimately requires O N time to build the lists as well as O N storage for the lists Moldy uses an implementation of the link cell method of Quentrec al 45 described in Allen and Tildesley s book 2 pp 149 152 which is a true O N algorithm The fundamental idea is that the MD cell is partitioned into a number of smaller cells known as subcells Every timestep a linked list of all the particles contained in each subcell is constructed The selection of all pairs of particles within the cutoff is achieved by looping over all pairs of subcells within the cutoff and particles within the subcells Because of their regular arran
156. t sequentially except if the a flag is given whereupon they are averaged If no options or only a or m are given columnar data is written to standard output If either of the plot program commands pg or xg is given the selected plotting program is invoked to display the plot The columnar RDF data is left in a file named according to the input file with data appended If either of the ps or eps options is given the selected plot program is used to create a PostScript or encapsulated PostScript file containing the plot The np option modifies all of the above actions by stopping the invocation of the plot program and writing the script to standard output instead Columns may be selected by matching with a PERL regular expression without the usual delimiters given as the argument to the m option 5Gnuplot can be downloaded from http www gnuplot org and Xmgrace may be obtained from http plasma gate weizmann ac il 65 Bibliography 1 M Abramowitz and I A Stegun Handbook of mathematical functions Dover New York 1970 2 Allen and D J Tildesley Computer simulation of liquids Clarendon Press Oxford 1987 3 Hans C Andersen Molecular dynamics simulation at constant pressure and or temperature J Chem Phys 72 1980 4 2384 2393 4 D Beeman Some multistep methods for use in molecular dynamics calculations 1 Comp Phys 20 1976 130 139 5 H Bekker E J Dijkstra M K R Re
157. ted file names This also happens if an existing file does not appear to be a Moldy dump Existing dump files are also tested to ensure that there is no corruption due for example to a system crash and that they contain the correct number of records If the dump sequence can not be continued in these circumstances Moldy terminates with a fatal error rather than waste computer time This is actually the code sprintf C library function uses to signify converting an integer to a decimal character string This function is used to create the actual file name from the prototype and the integer dump sequence number See any C library manual for details timestamp is simply the number of seconds elapsed since midnight on January 1 1970 28 Two utility programs included in the distribution are dumpanal which identifies dump files by printing out the headers and dumpext which extracts atomic or molecular trajectories The latter should be useful as a prototype for writing programs to analyse dump data It is frequently convenient to perform analysis of dump data and perhaps graphical output on a different computer to that which generated the data In the past it has not usually been possible to sensibly transfer binary data between computers of different architectures However Moldy is able to write dump files in a portable format using XDR see section which may be read by dumpext on any machine The control parameter xdr enables
158. tems and frequently very important on machines with interleaved memory i e most vector machines 4 If the value of any array element is used more than once in a loop write the loop using temporary scalars to store results and assign them to the arrays at the end of the loop This allows the compiler to optimize memory references P 5 Minimize the floating point operation count in critical loops 6 Minimize integer arithmetic in critical code CRAY vector machines in particular have no integer multiplication hardware and integer operations are slow as a result The performance of Moldy has been carefully studied using profiling tools and all critical regions of code are written as efficiently vectorizable loops The most critical sections of code i e those which use the majority of the computer time are all to do with the site forces calculation Thus it is the inner loops in force c ewald c and kernel c to which most attention should be paid The pair distance loop of rd calc in rdf c should vectorize for efficient radial distribution function evaluation Others which are of minor importance are in beeman c matrix c quaterns c and algorith c Auxil c contains alternative versions of various sum dot product scatter and gather routines efc which are interfaces to machine specific libraries e g Cray scilib Convex veclib which usually have FORTRAN calling conventions There are also default versions coded in C which do vectorize
159. text is not the directive described in the manual but rather that directive after it has been passed through the preprocessor To determine what should be substituted on a new vector machine create a small test containing the documented directive and use the C preprocessor on that file The output will show the form that should be defined in defs h Unfortunately this method had been made obsolete by the ANSI C standard which makes it impossible to insert pragmas using the preprocessor Potential names in system specification files are converted to lowercase before being compared with potspec so it is impossible to match name containing uppercase letters 8By default the arrays are sized for up to seven parameters If this is not sufficient the limit set by the value of the constant NPOTP defined in defs h may be increased 37 be necessary to modify this file new port should declare a new preprocessor macro along the lines of MPI etc which should be used to conditionally compile its code only Any header files may be included in the appropriate place in parallel c Then the interface functions should be written to call the underlying message passing library These should again be conditionally compiled It should be obvious where to place them in the file and the existing versions will provide a model Their specifications are par sigintreset void Moldy sets a handler for SIGINT This function is called from the signal handler to r
160. than are assumed to be negligible and need not be evaluated In the case of a polyatomic molecular system it is usual to apply the cutoff according to the intermolecular separation R 2 3 Pressure and Stress The internal stress of a system of interacting molecules is given by Allen and Tildesley D pp 46 49 in terms of the molecular virial but we may rewrite it more conveniently in terms of the atomic site virial as Suv Yn 2 15 i Env EE ofa i joi p The quantity Pia rio is the co ordinate of each site relative to the molecular centre of mass The pressure is easily evaluated as one third of the trace of the stress tensor 2 3 2 Long Range Corrections The truncation of the interactions at the cut off radius does introduce some errors into the calculated potential energy and stress By neglecting density fluctuations on a scale longer than the cutoff radius we may approximate the errors and calculate correction terms 2 pp 64 65 2m i U V Li 2 2 16 a _ 2T 3 PV 2 3 Uc NaNprebap re 2 17 a p Pj 2 18 where the sums run over all the distinct kinds of sites is the total number of sites of type in the system and is the unit matrix 2 3 3 Potential Functions Moldy supports most common forms of potential function including Lennard Jones o r 2 9 2 2 gt 6
161. the main design criteria of the code The general nature of the systems to be accepted by Moldy namely the arbitrary mixtures of molecules with different numbers and kinds of atoms requires a number of large multidimensional arrays with system dependent bounds in more than one dimension It is impractical to declare these statically with dimensions fixed at some suitably large value because total memory use would then be infeasibly large The availability of standard portable dynamic memory allocation was one of the major reasons the author chose to write Moldy in ch Moldy uses C s capability of array emulation by pointers and heap based memory 28 rather than true C arrays which like FORTRAN s are restricted to bounds fixed at compile time Much of the dynamic memory used in Moldy is allocated on entry to a function and deallocated before exit to emulate local variably dimensioned arrays The main exceptions are the arrays of dynamical variables which are allocated once during the startup phase and not freed until program exit All dynamic memory is allocated using the function talloc which is a wrapper around the standard library function malloc Function talloc simply calls malloc tests the return value and calls an error exit function if it failed to allocate the memory requested Its interface takes advantage ofthe C preprocessor macros LINE and FILE to print out the location of a failing call and a wrapping macro aalloc
162. thod employed in Moldy This second order formulation was first used by Powles et al 3 and Sonnenschein showed 53 that it gives substantially better stability than if angular velocities and accelerations are used as dynamic variables Using equations D TO to describe the dynamics means that they are integrated as if all four components were in dependent Therefore the normalization condition qd 1 may not be exactly satisfied after performing an integration step Moldy adopts the usual practice of scaling all components of the quaternion after each timestep to satisfy the normalization condition IIA It is less widely realized that the second order equations 2 10 introduce a second unconstrained variable into the procedure Differentiating equation 2 6 gives a constraint on the quaternion derivatives gogo 4141 q2d2 9393 0 2 11 which is just the qo component of equation 2 9 Just as with the normalization condition the integration algorithm will not preserve this condition exactly unless explicit measures are taken After each timestep the constraint may be re established by subtracting the discrepancy from the quaternion derivatives If 6 gogo 4141 424 q343 then the corrected quaternion derivatives are given by 4 2 12 The definition of quaternions used here differs from that used in Evans paper I 3 equation 21 in the sign of 42 This error has been compounded by subsequent authors 53 52
163. tions only in the respect that it must correctly handle parallelepiped shaped MD cells and therefore a triclinic grid of k vectors A list of vectors satisfying k gt 0 and k lt ke is pre computed to simplify the control condition of the main k vector loop Again computational efficiency dictates much of the structure of ewald The values of cos k rj and sin k r must be computed for every site i and every vector k ha Ic with k lt Since the cosine and sine functions are relatively expensive to evaluate the values of cos ha cos kb etc are precomputed for each site i and h 0 1 max k 0 1 and stored in arrays of dimension hmax nsites etc Then within the loop over k vectors the trigonometric addition formulae are used to construct q cos k and q sin k r using only arithmetic operations That task is delegated to function qsincos since certain compilers can optimize it better that way The self energy and charged system terms of equations 19 and 2 24 are also evaluated in ewald 5 4 Parallelization It has been predicted since the early 1980s that parallel computers would offer the very highest performance for scientific computing That vision has been somewhat slow in coming about in part due to the difficulty of writing programs to explicitly exploit unique and idiosyncratic architectures and the lack of portability or re usability of the resulting code It
164. ulation and analysis of dump data see section B7 They are easily compiled on unix systems using the makefile the command is make progname for each program or make utilities to make the lot They are written with unix systems in mind using unix style option arguments but ought to compile and run under VMS if defined as foreign commands Several of them require you to specify lists of integer numbers for example selected molecules timeslices etc which share a common syntax for such options The numbers 1 5 17 to 20 inclusive and 34 44 54 94 are selected by the command line argument t 1 5 17 20 34 100 10 Selectors are separated by commas and each one may be a range separated by a hyphen with an optional increment following a colon B 1 Moldyext It is usual to plot the instantaneous values of energy temperature efc during the course of a simulation for example to monitor the approach to equilibrium Moldyext processes the periodic output produced by Moldy see section B 5 and extracts this information from each timestep recorded It is presented in tabular form for input into plotting programs The command is moldyext f fields output filel output file2 where fields is a list in the selector format above The numbering runs from left to right ignoring newlines and blank columns so that the translational KE s in figure B T are fields 1 and 14 and the o component of the stress tensor is 30 Moldyext can currently only
165. ump range m solvent species n no of substitutions u solute species w mass q charge z symbol c o output file where the arguments have the following meanings s read a system specification file r read a restart file Only one of s or r may be given d read configurational data from a dump file given as a prototype name containing a printf format string see section BZ t range of records in the dump file specified in selector format to average positions over m the name of the species to be replaced i e the solvent species u the name of the substituting or solute species n the number of molecules to be substituted w mass of the solute species Defaults to same mass as the first atom of the solvent species q charge of the solute species Defaults to same charge as the first atom of the solvent species z symbol of the solute species Defaults to same symbol as the first atom of the solvent species c ifreading system specification file skip any preceding control parameter info o name of optional output file System specification data created from multiple time slices are written consecu tively to the same file If the species data come from a restart file the user is prompted for the energy units to use when writing the potential parameters to the output N B These may vary slightly from the original inputted parameters depending on the exact value of the time unit used in the
166. urce version which might work under Mingw32 but only the pamm version ynder Cyg win E been tested The compilation procedure is exactly the same as under a unix environment configure followed by make N B The executables so created must be run under either a Cygwin bash shell or a MS DOS window Moldy has also been successfully compiled using Watcom C and the development environment Consult the Makefile to see which objects to link to build Moldy and the utilities It is also possible that the Borland C Makefile mak will work under Borland C for Windows 95 or 2 sure to delete or rename the unix make file Makefile since Turbo C make will attempt to execute this in preference Makefile mak The author s experience of Windows and Windows 95 platforms is somewhat limited I would very much welcome any reports from users how to build Moldy in these environments for inclusion in future versions of this manual 33 4 1 1 Parallel Version Message Passing The parallel version of Moldy relies on an interface with a suitable message passing library This is the recommended version and supersedes the shared memory parallel implementation described in section 1 J The current release contains interfaces to the MPI 7 the TCGMSG library and the Oxford BSP library MPI is the most recom mended interface since it is the new standard for message passing libraries and should become ubiquitous If none o
167. utput c if called with the FATAL or WARNING flags set in its argument list In the serial case this prints out a failure message and calls exit to terminate the run In a parallel run an error condition may occur on just one processor or on all simultaneously depending on precisely what caused the problem Furthermore if just one processor is involved it may not be the I O processor rank 0 The approach of printing the message only from processor 0 is inadequate since there would be no indication that something has gone wrong But neither is it appropriate to print on all since that would result in a multiplicity of messages in some cases The approach taken is that message always prints the error irrespective of which processor is executing it but that calls to message are made conditional on ithread 0 for those conditions known to occur on all processors simultaneously If the condition signalled was FATAL then message finally calls abort which terminates the run on all processors 5 44 Distributed forces calculation The two parts of the forces calculation which must be parallelized are the short ranged forces calculation and the reciprocal space part of the Ewald sum The outermost loop of the link cell forces calculation in function force inner runs over all subcells within the simulation cell see figure 5 In the serial case it has the form for icell 0 icell lt nx ny nz 11 Evaluate forces on
168. ver rot Get torque corr Nose Hoover None nhtherm Get N H G Thermostat euler Gaussian Calc angular accs f Iteratively update quatern derivatives using Beeman gtherm beeman_2 Get gaussian Step G step iv beeman 2 step quatern derivs hoover tr Calc N H acc corr e Iteratative loop to step quatern N C of M velocities using derivatives Beeman step iv converged beeman_2 step C of M vels step 2 Apply Beeman step iv to dynamic var derivatives velocities converged every dump interval steps Figure 5 4 b Flow diagram of function do step which performs a single timestep continued from Figure Qump write trajectories forces etc any polyatomics N 49 first call cell size shape changed Set up arrays of P B C relocation vectors Y neighbour list Create list of neighbouring cells site neighbour list force inner loop over all cells site neighbour list Build list of all sites within range of current cell Gather site co ords reloc vectors into neighbour lists Allocate and fill expanded array of pot params fill cells Build linked list of mols in cells Gather
169. version of a file These 35 functions make use of the file name syntax of the host operating system and are therefore system dependent Unix file systems do not have explicit version numbers but Moldy keeps a single previous version by appending a character to the name The pure ANSI versions just interface to rename and do nothing respectively 4 2 20 Optimization and Vectorization Moldy has been designed to run fast on a wide range of computers and in particular on those with vector multiprocessor and parallel architectures This is a difficult problem since the constructs which run fast on different architectures may be quite distinct and occasionally in conflict Nonetheless it has been found that following a few basic rules gives extremely good performance on a wide range of computer architectures In a rough order of importance these are 1 Minimize the number of memory references to floating point data in critical loops Memory access is the major bottleneck on almost every modern computer scalar vector or parallel 2 Minimize the number of memory references to floating point data in critical loops This cannot be emphasized enough 3 Ensure that memory is accessed contiguously within critical loops That is arrays should be accessed with a stride of 1 and with the last index varying most rapidly f This is absolutely critical on machines where memory is accessed via a cache i e all workstations and many parallel sys
170. w iterations into cache Otherwise a very slow cache load would be necessary on every loop iteration The link cell functions are all in force c figure 5 The outer level function is orce calc Function neighbour list constructs a list of subcells within the cutoff radius of a reference cell By adding an offset and applying periodic bound ary conditions this list yields the neighbour cells of subcell This should not be confused with the neighbour sites list above In fact the list contains only a hemisphere of cells since Newton s third law is exploited to halve the compu tational cost more rigorous list may be computed instead by strict neighbour list which is selected in strict cutoff mode see section 2 5 3 Next an array of potential parameters called potp is constructed This has an innermost dimension of the number of sites which maps one to one onto the co ordinates array It will be used as the argument of the gather operation using the site neighbour list Function 111 cells constructs the linked list of molecules in each subcell The loops over cells are all contained in force inner Within the loop over subcells the neighbour list of sites for this cell is constructed from the list of neighbour cells by site neighbour list This list is then used as the index array to gather the co ordinates charges potential parameters and any periodic boundary vectors The innermost loops compute the site pair distances
171. ystems it will be lost 5A backup file is also written if the run is terminated for exceeding its cpu limit This is not quite true Moldy does read the control file and any restart file but only to determine the name of the backup file Thus even if the backup has a non standard name it can still be found 24 from the backup name by appending lck A second run which attempts to use the same backup file will test for the presence of the lock file and abort the run if it finds it A restart or backup file is created by first writing the data to a temporary file which is then renamed to the final name This ensures that there is no possibility of a file being left incomplete or corrupt if the computer crashes part way through the write If the file already exists either it is replaced on systems which only keep one version of a file or new version is created on systems such as VMS which retain multiple versions In the unlikely event of it being necessary to change where the temporary file is it may be specified with the control parameter temp file 3 3 2 Portable Restart Files Moldy is ablef to read and write restart and dump files in a portable binary format which is transportable between computers of different architectures So a restart file written on for example a Sun may be used to initiate a new run on a Cray and the dump files generated on the Cray may be analyzed on the Sun This feature will also be of considerable us

Moldy User's Manual

Contents

Download Pdf Manuals

Related Search

Related Contents