1 /* symop.cpp: fundamental data types for the clipper libraries */
2 //C Copyright (C) 2000-2006 Kevin Cowtan and University of York
3 //L
4 //L This library is free software and is distributed under the terms
5 //L and conditions of version 2.1 of the GNU Lesser General Public
6 //L Licence (LGPL) with the following additional clause:
7 //L
8 //L `You may also combine or link a "work that uses the Library" to
9 //L produce a work containing portions of the Library, and distribute
10 //L that work under terms of your choice, provided that you give
11 //L prominent notice with each copy of the work that the specified
12 //L version of the Library is used in it, and that you include or
13 //L provide public access to the complete corresponding
14 //L machine-readable source code for the Library including whatever
15 //L changes were used in the work. (i.e. If you make changes to the
16 //L Library you must distribute those, but you do not need to
17 //L distribute source or object code to those portions of the work
18 //L not covered by this licence.)'
19 //L
20 //L Note that this clause grants an additional right and does not impose
21 //L any additional restriction, and so does not affect compatibility
22 //L with the GNU General Public Licence (GPL). If you wish to negotiate
23 //L other terms, please contact the maintainer.
24 //L
25 //L You can redistribute it and/or modify the library under the terms of
26 //L the GNU Lesser General Public License as published by the Free Software
27 //L Foundation; either version 2.1 of the License, or (at your option) any
28 //L later version.
29 //L
30 //L This library is distributed in the hope that it will be useful, but
31 //L WITHOUT ANY WARRANTY; without even the implied warranty of
32 //L MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
33 //L Lesser General Public License for more details.
34 //L
35 //L You should have received a copy of the CCP4 licence and/or GNU
36 //L Lesser General Public License along with this library; if not, write
37 //L to the CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK.
38 //L The GNU Lesser General Public can also be obtained by writing to the
39 //L Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
40 //L MA 02111-1307 USA
41
42
43 #include "symop.h"
44
45 #include "coords.h"
46
47
48 namespace clipper {
49
50
51 /*! Construct an RT operator from a string description,
52 e.g. 1/2x,z-y+2/3,x
53 '*' is optional for multiplication, commas are compulsory. */
RTop_frac(const String & strop)54 RTop_frac::RTop_frac( const String& strop )
55 {
56 std::vector<String> rows = strop.split(",");
57 if ( rows.size() != 3 )
58 Message::message(Message_fatal("RTop_frac: invalid string:"+strop));
59
60 rot() = Mat33<>(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
61 trn() = Vec3<> (0.0,0.0,0.0);
62
63 std::vector<String> cols;
64 String num;
65 ftype val;
66 int nrow, ncol, nchr, npart;
67 for ( nrow = 0; nrow < 3; nrow++ ) {
68 const String& row = rows[nrow];
69 cols.clear(); cols.push_back("");
70 for ( nchr = 0; nchr < row.length(); nchr++ ) {
71 const char& c = row[nchr];
72 if ( c == '+' || c == '-' ) cols.push_back( "" );
73 cols.back() += c;
74 }
75 for ( npart = 0; npart < cols.size(); npart++ ) {
76 const String& col = cols[npart];
77 ncol = 3;
78 num = "";
79 for ( nchr = 0; nchr < col.length(); nchr++ ) {
80 const char& c = col[nchr];
81 if ( c == 'x' || c == 'X' ) ncol = 0;
82 else if ( c == 'y' || c == 'Y' ) ncol = 1;
83 else if ( c == 'z' || c == 'Z' ) ncol = 2;
84 else if ( c == '-' || c == '/' || isdigit(c) ) num += c;
85 }
86 if ( ncol < 3 ) {
87 if ( num == "" || num == "+" || num == "-" ) num += "1";
88 val = num.rational();
89 rot()( nrow, ncol ) = val;
90 } else {
91 if ( num == "" || num == "+" || num == "-" ) num += "0";
92 val = num.rational();
93 trn()[ nrow ] = val;
94 }
95 }
96 }
97 }
98
99 /*! \param cell The cell concerned \return The transformed coordinate. */
rtop_orth(const Cell & cell) const100 RTop_orth RTop_frac::rtop_orth( const Cell& cell ) const
101 {
102 return RTop_orth( RTop<>(cell.matrix_orth()) * (*this) * RTop<>(cell.matrix_frac()) );
103 }
104
105 /*! \return The inverse of the operator. */
inverse() const106 RTop_frac RTop_frac::inverse() const
107 { return RTop_frac( RTop<>::inverse() ); }
108
109 /*! \return The identity operator. */
identity()110 RTop_frac RTop_frac::identity()
111 { return RTop_frac( RTop<>::identity() ); }
112
113 /*! \return The null (uninitialised) operator. */
null()114 RTop_frac RTop_frac::null()
115 { return RTop_frac( RTop<>::null() ); }
116
117
118 /*! Construct a symmetry operator and initialise it to the supplied RTop.
119 Translations are rounded to a basis of 48, and put on the range 0..1
120 \param mat The RTop to use. */
Symop(const RTop<> & rt)121 Symop::Symop( const RTop<>& rt )
122 {
123 // initialise to the supplied matrix
124 int i, j;
125 for ( i=0; i<3; i++) for ( j=0; j<3; j++)
126 rot()(i,j) = rint( rt.rot()(i,j) );
127 for ( i=0; i<3; i++)
128 trn()[i] = ftype( Util::mod( Util::intr(48.*rt.trn()[i]), 48 ) ) / 48.;
129 }
130
131 /*! Construct a symmetry operator and initialise it to the supplied matrix.
132 Translations are rounded to a basis of 48, and put on the range 0..1
133 \param mat The 4x4 matrix to use. The [i][3] elements contain the
134 translation. */
Symop(const ftype mat[4][4])135 Symop::Symop( const ftype mat[4][4] )
136 {
137 // initialise to the supplied matrix
138 int i, j;
139 for ( i=0; i<3; i++) for ( j=0; j<3; j++)
140 rot()(i,j) = mat[i][j];
141 for ( i=0; i<3; i++)
142 trn()[i] = ftype( Util::mod( Util::intr(48.*mat[i][3]), 48 ) ) / 48.;
143 }
144
145 /*! Return formatted representation of the symmetry operator.
146 \return The formatted text string, e.g. -x, -y+1/2, z. */
format() const147 String Symop::format() const
148 {
149 String s, t, xyz="xyz";
150 for ( int i = 0; i < 3; i++ ) {
151 t = "";
152 for ( int j = 0; j < 3; j++ )
153 if ( rot()(i,j) != 0.0 ) {
154 t += ( rot()(i,j) > 0.0 ) ? "+" : "-";
155 if ( Util::intr( fabs( rot()(i,j) ) ) != 1 )
156 t += String::rational( fabs( rot()(i,j) ), 24 );
157 t += xyz[j];
158 }
159 if ( trn()[i] != 0.0 )
160 t += String::rational( trn()[i], 24, true );
161 s += t.substr( ( t[0] == '+' ) ? 1 : 0 );
162 if ( i < 2 ) s+= ", ";
163 }
164 return s;
165 }
166
167
168 /*! Integerised symops are more efficient when handling integer
169 coordinate types, e.g. HKL, Coord_grid. The rotation parts of the
170 integerised symop are general and can be used for any recirpocal
171 space data. The translation part is specific to an individual grid.
172 \param symop The conventional symop.
173 \param grid The specific grid. */
Isymop(const Symop & symop,const Grid & grid)174 Isymop::Isymop( const Symop& symop, const Grid& grid )
175 {
176 for ( int i = 0; i < 3; i++ )
177 for ( int j = 0; j < 3; j++ )
178 rot()(i,j) = Util::intr( symop.rot()(i,j) );
179 trn()[0] = Util::intr( grid.nu() * symop.trn()[0] );
180 trn()[1] = Util::intr( grid.nv() * symop.trn()[1] );
181 trn()[2] = Util::intr( grid.nw() * symop.trn()[2] );
182 }
183
184
185 /* Construct the encoded form of the symop. */
Symop_code(const Symop & op)186 Symop_code::Symop_code( const Symop& op )
187 { init( Isymop( op, Grid(24,24,24) ) ); }
188
189 /* Construct the encoded form of the symop from integerised symop with
190 a grid (base) of (24,24,24).*/
Symop_code(const Isymop & op)191 Symop_code::Symop_code( const Isymop& op )
192 { init( op ); }
193
194 /* Initialise the encoded form of the symop from integerised symop with
195 a grid (base) of (24,24,24).*/
init(const Isymop & op)196 void Symop_code::init( const Isymop& op )
197 {
198 // initialise to the supplied code
199 int i, j, fac, code_r, code_t;
200 code_r = code_t = 0;
201 fac = 1;
202 for ( i = 0; i < 3; i++ ) {
203 code_t += ( Util::mod( op.trn()[i], 24 ) ) * fac;
204 fac *= 24;
205 }
206 fac = 1;
207 for ( i = 0; i < 3; i++ ) for ( j = 0; j < 3; j++ ) {
208 code_r += ( Util::mod( op.rot()(i,j) + 1, 3 ) ) * fac;
209 fac *= 3;
210 }
211 // xor to make identity zero
212 code_ = ( ( code_r ^ 0x4064 ) << 16 ) + code_t;
213 }
214
215 /*! Construct a symmetry operator and initialise it to the matrix
216 encoded in the given int.
217 \param code The integer code. */
symop() const218 Symop Symop_code::symop() const
219 {
220 Isymop iop = isymop();
221 Symop op;
222 for ( int i = 0; i < 3; i++ ) {
223 op.rot()(i,0) = ftype( iop.rot()(i,0) );
224 op.rot()(i,1) = ftype( iop.rot()(i,1) );
225 op.rot()(i,2) = ftype( iop.rot()(i,2) );
226 op.trn()[i] = ftype( iop.trn()[i] ) / 24.0;
227 }
228 return op;
229 }
230
231 /*! Construct an integerised symmetry operator and initialise it to
232 the matrix encoded in the given int, with a grid (base) of (24,24,24).
233 \param code The integer code. */
isymop() const234 Isymop Symop_code::isymop() const
235 {
236 Isymop op;
237 // initialise rotation and translation
238 int i, j, fac, code_r, code_t;
239 code_t = code_trn();
240 fac = 1;
241 for ( i = 0; i < 3; i++ ) {
242 op.trn()[i] = Util::mod( code_t/fac, 24 );
243 fac *= 24;
244 }
245 code_r = ( code_rot() >> 16 ) ^ 0x4064;
246 fac = 1;
247 for ( i = 0; i < 3; i++ ) for ( j = 0; j < 3; j++ ) {
248 op.rot()(i,j) = Util::mod( code_r/fac, 3 ) - 1;
249 fac *= 3;
250 }
251 return op;
252 }
253
code_rot() const254 Symop_code Symop_code::code_rot() const
255 { return Symop_code( code_ & 0xffff0000 ); }
256
code_trn() const257 Symop_code Symop_code::code_trn() const
258 { return Symop_code( code_ & 0x0000ffff ); }
259
260
261 } // namespace clipper
262