GASP manual
From GridInfo
GASP manual. Stephen Wells, Royal Institution, London W1S 4BS. stephen@ri.ac.uk
GASP (GEOMETRIC ANALYSIS OF STRUCTURAL POLYHEDRA) was written by Stephen Wells to perform real-space rigid unit analyses on framework structures using Geometric Algebra. It is freely available to the scientific community. Your first publication using GASP should list Dr. Wells as a co-author; subsequent independent publications should reference [thesis].
1) Functions of GASP
GASP analyses mineral structures from the point of view of the Rigid Unit model. The primary aim of the program is to compare two polyhedral framework structures and analyse their differences as a combination of rigid unit motion (displacement and rotation of polyhedra) and distortion of polyhedra. The program will also perform a 'geometric relaxation' on a structure by attempting to make the polyhedra as regular as possible. This allows a purely geometrical investigation of the framework response to defects or to compression.
GASP is written entirely in Fortran 90 and does not make use of any libraries; it should therefore be easy to compile and maintain.
GASP has six modules (basefunctions, global, local, fitting, routines and gaspmain) which should be compiled in this order:
user@work: f90 -c basefunctions.f90 global.f90 local.f90 user@work: f90 -c fitting.f90 routines.f90 gaspmain.f90 user@work: f90 -o gasp *.o
to produce the gasp executable.
2) Files
All GASP files are text files. The program is controlled by a configuration file (foo.cfg) which contains geometrical information and commands. Structures can be read in native .dat format or in the common .xtl format; a parser is in development for CML structure files. Topological information in the form of lists of bonded neighbours is given in a .con file, while a representation of the structure as a polyhedral framework is kept in a .pol file. Two .pol files, when compared, generate an analysis of rigid unit and distortive motions, which is given in a .com file. The relaxer will output structures in .dat or .xtl format.
3) GASP commands in the .cfg file
On running the GASP executable, the program will prompt you for the name of a .cfg file. Enter foo and the program will read commands from foo.cfg.
The first line of the .cfg file is a title/comment line. The second must list all the elements (atom names) present in the structure. These names are restricted to two letters and need not bear any resemblance to real element names.
Geometrical information is given as polyhedral specifications, one to a line, following the element list. A specification includes the shape of the polyhedron (tetrahedron, octahedron or planar triangle), the names of the centre and vertex species, and the ideal centre-vertex bond length.
The remaining lines take the form "COMMAND arguments". Only the first three letters of a command are significant and case is not significant (SET, Set and set are all equivalent). Arguments are filenames or program parameters and should be given exactly. Some commands are switches without arguments. Blank lines are skipped. Execution is terminated either by one of the switches (END, EOF) or by reaching the end of the .cfg file. There is no specific comment syntax, but comments may be placed on the title line; after the complete argument list for a command; or after an END or EOF statement.
The four principal commands are PREcondition, FAStprecondition, COMpare and RELax. Principal switches are END, EOF, XTL, XML, CML. Parameters are adjusted using the SET command discussed in section [9].
-> Table of commands and arguments. <- COMMAND ARGUMENT1 ARGUMENT2 ARGUMENT3 Notes TET Centre species Vertex species Distance Tetrahedron, bond length=distance OCT Centre species Vertex species Distance Octahedron, bond length=distance TRI Centre species Vertex species Distance Triangle, bond length=distance PRE file1 Preconditioning. Reads structure from file1.dat Generates file1.con, file1.pol FAS file1 file2 Fast preconditioner. Reads file1.dat and file2.con Generates file1.pol COM file1 file2 Comparison. Reads file1.pol and file2.pol Generates file1_file2.com REL file1 file2 file3 Relaxation Reads file1.dat, file2.con, file3.for Generates file1R.dat END Ends execution EOF Ends execution SET variable value Resets parameters; see section 9 XML Not yet implemented XTL Use .xtl instead of .dat structures DAT Use .dat structures WHY ELV TAG
4) Structures in the .dat file
The structural information required is the dimensions of the unit cell and the names and positions of all atoms within the cell. Atom names are read as two characters only and need not adhere to the periodic table. The first line of a .dat file is a comment; the next line(s) should be either six cell parameters or a set of three Cartesian cell vectors. The remaining lines should each consist of a species and three coordinates.
If structures are given in .xtl format, they should be in P 1 symmetry (as is usual for a supercell). The command XTL should be placed in the .cfg file before any command reads a .xtl file.
In native .dat format, the unit cell may be given either by six crystallographic parameters or by three cartesian cell vectors. When cell vectors are used they should adhere to the following convention: the j cell vector should lie along the cartesian y axis; the k cell vector should lie in the cartesian y-z plane and extend in the positive z direction; the i cell vector is general but should of course extend in the positive x direction. Since large numbers of atoms are usually involved in a geometric analysis, the parameters for a crystal will usually be those of a supercell rather than a single unit cell.
Positions may be given in a .dat file either as fractional or cartesian coordinates. If fractional coordinates are given they should lie between 0 and 1, and the word 'fractional' should appear in the .dat file on a line between the cell vectors and the atomic positions.
5) Polyhedra and the .pol file
The .pol file contains a breakdown of a structure into polyhedra; for each polyhedron it gives about the positions of the atoms, and also the ideal geometric positions of the vertices; the difference between the real and ideal positions of the vertex atoms quantifies the distortion of the polyhedron. These mismatch scores (see section [7]) are summarised by polyhedral type at the end of the .pol file.
Example of a .pol file:
Polyhedra for ed_0_689
19.120399 -0.000001 -0.000001
0.000000 19.120399 0.000000
0.000000 -0.000001 13.023400
1 TET Si O 1.650000 1.647238 0.000339 0.000031
1 232 181 187 234
72 21 27 74
0.00000 0.00000 0.00000
-0.41854 -1.28374 -0.94355 -0.41759 -1.28088 -0.95261
-1.28374 0.41855 0.94342 -1.28088 0.41764 0.95259
1.28374 -0.41854 0.94355 1.28083 -0.41759 0.95267
0.41855 1.28374 -0.94355 0.41764 1.28084 -0.95265
Breakdown of .pol file:
Polyhedra for ed_0_689 <- comment line
=
19.120399 -0.000001 -0.000001 <- Cell vector i
0.000000 19.120399 0.000000 <- Cell vector j
0.000000 -0.000001 13.023400 <- Cell vector k
=
1 TET Si O 1.650000... <- First polyhedron, Si-O tetrahedron, 1.65 A bond length...
... 1.647238 0.000339 0.000031
... RMS bond length, angle-bending sum of squares, bond-stretching sum of squares
=
1 232 181 187 234 <- Member atoms: 1 (centre),232,182,187 and 234
72 21 27 74 <- connected to polyhedra 72,21,27 and 74
=
0.00000 0.00000 0.00000 <- position of central atom
=
-0.41854 -1.28374 -0.94355... <- position of vertex atom...
... -0.41759 -1.28088 -0.95261
... position of geometrically ideal vertex
=
-1.28374 0.41855 0.94342 -1.28088 0.41764 0.95259
1.28374 -0.41854 0.94355 1.28083 -0.41759 0.95267
0.41855 1.28374 -0.94355 0.41764 1.28084 -0.95265
Positions of vertex atoms and ideal vertices
=
Summary information at tail of .pol file:
Global results for ed_0_689R <- sum of results for all polyhedra Global mismatch value: 6.4016677E-02 <- mismatch between atoms and ideal vertices Bond-angle-bending term: 4.1426696E-02 Bond stretching term: 2.2589978E-02 Type 1 : Si O ; 80 polyhedra. <- breakdown by polyhedral type Ideal bondlength : 1.650000 RMS bond length : 1.646894 Type mismatch value : 6.4016677E-02 Type bond-angle-bending : 4.1426696E-02 Type bond stretching : 2.2589978E-02 RMS distortion distance : 1.4143978E-02 <- typical atom-to-vertex distance RMS bond-angle-bending distance : 1.1377980E-02 <- distortion due to angle bending RMS bond stretching distance : 8.4020048E-03 <- typical variation in bond length
6) Topology and the .con file
The connectivity of the structure in defined in the .con (connections) file. For each atom it lists its role as polyhedral centre, vertex or interstitial atom; vertices are labelled VT and central atoms are labelled TC, OC or 3C for tetrahedral, octahedral or triangular centres, while interstitials are labelled IN. This is followed by the numbers of the atoms which the current atom is connected to. Vertices may be non-bridging, in which case one of these numbers will be zero, and interstitials have no connected atoms.
The file lists the species of each atom but this information is not used by the program; the same .con file can therefore be used for structures with different elements but the same topology. This is useful for the study of aluminosilicates or AlPOs derived from silicate structures.
The preconditioning routine generates a .con file and a .pol file from the input .dat or .xtl file. The fast-preconditioning routine generates a .pol file from the input .dat or .xtl and an input .con file.
-> Example of a .con file <- Connections for ed_zeroP <- comment line Si TC 232 181 187 234 <- first atom, Tetrahedral Centre, bonds to atoms 232, 181, 187 and 234 Si TC 233 231 182 188 Si TC 185 230 183 236 ...
7) Comparison of structures; the .com file.
The comparison of two *.pol files generates the set of translations and rotations bringing the polyhedra from the first file as close as possible to their counterparts in the second file, and quantifies the residual distortion.
For each polyhedron, the comparison file lists the displacement (vector) and the rotation (bivector components) on one line, and the mismatch scores on the second line in the following order:
- Mismatch before fitting, the sum of the squares of the differences in position for each atom in the polyhedron.
- Mismatch with centres fitted, the sum of squares of differences in vertex atom positions when the centres of the polyhedra coincide. If rigid-unit translations are significant, this will lie below the first value.
- Mismatch after fitting, the residual sum of squares after finding the best-fit rotation.
- Angle-bending term, that part of the mismatch due to changes in internal angles.
- Stretching term, that part of the mismatch due to changes in bond length.
At the end of the file all of this information is summarised, both globally (for all polyhedra) and broken down by type of polyhedron, if more than one type is present.
7a) Comparison of results: RMS motion and mismatch scores.
The initial mismatch score is the sum of squares of atomic motion from one configuration to the other. This can be decomposed by atomic species and we can obtain the RMS motion for each species, which we may write as Delta-r.
The final mismatch score is the sum of squares of differences between vertex positions once all rigid-unit motion is accounted for- that is, the sum of squares of distortive motions. We can decompose this score by polyhedral type and obtain the RMS distortive motion of a vertex atom, which we may write as d.
Comparison of the values of Delta-r and d for the vertex atoms thus indicates the significance of rigid unit motion. For example, when analysing quartz structures obtained by Reverse Monte Carlo fitting to total neutron scattering data, we find that at room temperature Delta-r is approximately twice d, whereas at high temperatures in the beta phase Delta-r is approximately four times d.
Look in the summary information at the end of the .com file for the entries "RMS motion (distance)" for each atomic species and "RMS distortive motion (distance)" for each polyhedral type.
8) Geometric relaxation and the .for file
Geometric relaxation simply means relaxing the atoms so as to fit the ideal geometric forms of the polyhedra. Having fitted ideal polyhedra over the real atoms, we can apply a spring model tending to move vertex atoms towards the ideal geometric positions of the vertices, along with centre-centre springs which will optimise the bridging angles. After relaxing the atoms on these springs, we can re-fit the ideal polyhedra and iterate until the structure is relaxed. This offers a simple means of investigating the response of a framework to size effects, defects or compression.
The relaxation of the atoms on the springs is performed by an approximate single-stpe method, which eliminates any representation of mass from the model. The values of the spring constants are therefore artbitrary, onyl their ratios being significant.
Relaxation is controlled by a force file (*.for) which contains spring constants constraining bending and stretching motions for each polyhedral type, along with optional centre-centre spring constants and ideal distances to constrain bridging angles. Atoms can be held in place with the PIN command, and polyhedra can be fixed in place with the FIX command and rotated with the TURn command.
After relaxing for a number of iterations controlled by the variable niterations (see section [9]), the resulting structure is written to a file *R.dat or *R.xtl.
-> An example .for file <- Force model for silicates Si O 1.62 50000 200000 Bond length 1.62 A, k for bending motion 50000, k for bond stretching 200000 bridging Si O Si 3.09 100000 Si-Si distance 3.09 A, k for stretching 100000
9) Adjustable parameters and the SET command
"SET variablename variablevalue" in the .cfg file will set the value of certain internal variables.
Do not fiddle with these unless you really mean it. The only one that might regularly need to be changed is "niterations", which governs the number of steps in the relaxation routine.
-> Table of adjustable parameters <- NAME DEFAULT FUNCTION maxpolys 10000 Maximum number of polyhedra maxatoms 50000 Maximum number of atoms findwithin 1.3 Search for bonded atoms within (bond length * findwithin) bigshift 4.00 A difference in position greater than this will be considered to be out of cell, triggering periodic boundary conditions maxiter 100 Maximum iterations during rotor-fitting algorithm small 0.001 Fit is achieved when the gradient is less than this basestep 0.15 Step size in rotor-fitting is basestep * gradient scale 10.0 Basestep is halved over this many iterations. niterations 40 Number of iterations in relaxation routine.
10) Mathematical background and algorithms.
Fitting: the procedure for comparison of polyhedra is based on identifying a displacement and a rotation. The reported displacement of a polyhedron is the displacement of its central atom. Once polyhedra have been displaced so that their centres coincide, the rotation is obtained by a least-squares fit minimising the mismatch between one polyhedron and the other. The parameters which are varied in the fit are the components of a rotor operator defining the rotation. The process of fitting is illustrated in the diagram "Fitting".
Boundary conditions: the preconditioning function searches for atoms forming polyhedra by seeking the nearest neighbours of each polyhedral-centre atom. If it cannot find sufficient vertex atoms within the cell, it will apply periodic boundary conditions. This behaviour can be disabled by including the word "Cluster" on its own line in the .dat file between the cell vectors and the atomic positions.
Since the program was designed for the analysis of large supercells, there are cut-offs based on the position of the atom within the cell; for example, an atom in the upper half of the cell will not look for neighbours in the periodic image directly below the cell. This can cause confusion when the cell is small; as a guideline, cell vectors should be at least 8-10 Angstroms long. When analysing smaller cells, it is best to generate a supercell.
10a) Geometric Algebra.
Rotor operator: the operator has the form R=X-B/2 and acts on a vector via the two-sided operation r' = R r R~. X is a scalar equal to cos(theta/2) and B is a bivector defining the plane and magnitude of the rotation; the magnitude of B is 2 sin(theta/2). The bivector B=bx B1 + by B2 + bz B3 gives a rotation about an axis parallel to bx e1 + by e2 + bz e3. The product of the rotor with its reverse is R R~ = 1.
Bivector: from the Cartesian basis vectors e1,e2,e3 we obtain the bivectors using the geometric product; the three basis bivectors in three dimensions are B1 = e2 e3, B2 = e3 e1 and B3 = e1 e2. Bivectors commute with orthogonal vectors and rotate coplanar vectors by ninety degrees. The rotor operator has scalar and bivector components. It produces rotations of any vector lying in the plane of the bivector, and commutes with the perpendicular vector. The two-sided action of the rotor and its reverse (r' = R r R~) produces a rotation.
Geometric product: the geometric product is written by juxtaposition ab and is associative and distributive but not necessarily commutative. In three dimensions, a Cartesian basis set e1,e2,e3 generates an algebra of eight objects:
- the scalar, 1
- the vectors, e1 e2 e3
- the bivectors, e1e2 e2e3 e3e1, also written Bz Bx By
- the trivector (pseudoscalar), e1e2e3
Each object has a reverse, which is the same vectors multiplied in the opposite order. The reverse of the rotor operator R is R~ ; R R~ = 1.
11) Help and advice on using GASP
All queries should be addressed to the author and maintainer, Stephen Wells (stephen@ri.ac.uk, swel99@esc.cam.ac.uk). Bug reports are welcome; fixes even more so.
