1
#include <clipper/clipper.h>
2
#include <clipper/clipper-contrib.h>
3
#include <clipper/clipper-cctbx.h>
4
//#include <clipper/clipper-mmdbold.h>
9
using namespace clipper;
14
const char* hallsymbols[] = {
15
// the 530 tabulated settings
16
"P 1", "-P 1", "P 2y", "P 2", "P 2x",
17
"P 2yb", "P 2c", "P 2xa", "C 2y", "A 2y",
18
"I 2y", "A 2", "B 2", "I 2", "B 2x",
19
"C 2x", "I 2x", "P -2y", "P -2", "P -2x",
20
"P -2yc", "P -2yac", "P -2ya", "P -2a", "P -2ab",
21
"P -2b", "P -2xb", "P -2xbc", "P -2xc", "C -2y",
22
"A -2y", "I -2y", "A -2", "B -2", "I -2",
23
"B -2x", "C -2x", "I -2x", "C -2yc", "A -2yac",
24
"I -2ya", "A -2ya", "C -2ybc", "I -2yc", "A -2a",
25
"B -2bc", "I -2b", "B -2b", "A -2ac", "I -2a",
26
"B -2xb", "C -2xbc", "I -2xc", "C -2xc", "B -2xbc",
27
"I -2xb", "-P 2y", "-P 2", "-P 2x", "-P 2yb",
28
"-P 2c", "-P 2xa", "-C 2y", "-A 2y", "-I 2y",
29
"-A 2", "-B 2", "-I 2", "-B 2x", "-C 2x",
30
"-I 2x", "-P 2yc", "-P 2yac", "-P 2ya", "-P 2a",
31
"-P 2ab", "-P 2b", "-P 2xb", "-P 2xbc", "-P 2xc",
32
"-P 2ybc", "-P 2yn", "-P 2yab", "-P 2ac", "-P 2n",
33
"-P 2bc", "-P 2xab", "-P 2xn", "-P 2xac", "-C 2yc",
34
"-A 2yac", "-I 2ya", "-A 2ya", "-C 2ybc", "-I 2yc",
35
"-A 2a", "-B 2bc", "-I 2b", "-B 2b", "-A 2ac",
36
"-I 2a", "-B 2xb", "-C 2xbc", "-I 2xc", "-C 2xc",
37
"-B 2xbc", "-I 2xb", "P 2 2", "P 2c 2", "P 2a 2a",
38
"P 2 2b", "P 2 2ab", "P 2bc 2", "P 2ac 2ac", "P 2ac 2ab",
39
"C 2c 2", "A 2a 2a", "B 2 2b", "C 2 2", "A 2 2",
40
"B 2 2", "F 2 2", "I 2 2", "I 2b 2c", "P 2 -2",
41
"P -2 2", "P -2 -2", "P 2c -2", "P 2c -2c", "P -2a 2a",
42
"P -2 2a", "P -2 -2b", "P -2b -2", "P 2 -2c", "P -2a 2",
43
"P -2b -2b", "P 2 -2a", "P 2 -2b", "P -2b 2", "P -2c 2",
44
"P -2c -2c", "P -2a -2a", "P 2c -2ac", "P 2c -2b", "P -2b 2a",
45
"P -2ac 2a", "P -2bc -2c", "P -2a -2ab", "P 2 -2bc", "P 2 -2ac",
46
"P -2ac 2", "P -2ab 2", "P -2ab -2ab", "P -2bc -2bc", "P 2ac -2",
47
"P 2bc -2bc", "P -2ab 2ab", "P -2 2ac", "P -2 -2bc", "P -2ab -2",
48
"P 2 -2ab", "P -2bc 2", "P -2ac -2ac", "P 2c -2n", "P 2c -2ab",
49
"P -2bc 2a", "P -2n 2a", "P -2n -2ac", "P -2ac -2n", "P 2 -2n",
50
"P -2n 2", "P -2n -2n", "C 2 -2", "A -2 2", "B -2 -2",
51
"C 2c -2", "C 2c -2c", "A -2a 2a", "A -2 2a", "B -2 -2b",
52
"B -2b -2", "C 2 -2c", "A -2a 2", "B -2b -2b", "A 2 -2",
53
"B 2 -2", "B -2 2", "C -2 2", "C -2 -2", "A -2 -2",
54
"A 2 -2c", "B 2 -2c", "B -2c 2", "C -2b 2", "C -2b -2b",
55
"A -2c -2c", "A 2 -2a", "B 2 -2b", "B -2b 2", "C -2c 2",
56
"C -2c -2c", "A -2a -2a", "A 2 -2ac", "B 2 -2bc", "B -2bc 2",
57
"C -2bc 2", "C -2bc -2bc", "A -2ac -2ac", "F 2 -2", "F -2 2",
58
"F -2 -2", "F 2 -2d", "F -2d 2", "F -2d -2d", "I 2 -2",
59
"I -2 2", "I -2 -2", "I 2 -2c", "I -2a 2", "I -2b -2b",
60
"I 2 -2a", "I 2 -2b", "I -2b 2", "I -2c 2", "I -2c -2c",
61
"I -2a -2a", "-P 2 2", "P 2 2 -1n", "-P 2ab 2bc", "-P 2 2c",
62
"-P 2a 2", "-P 2b 2b", "P 2 2 -1ab", "-P 2ab 2b", "P 2 2 -1bc",
63
"-P 2b 2bc", "P 2 2 -1ac", "-P 2a 2c", "-P 2a 2a", "-P 2b 2",
64
"-P 2 2b", "-P 2c 2c", "-P 2c 2", "-P 2 2a", "-P 2a 2bc",
65
"-P 2b 2n", "-P 2n 2b", "-P 2ab 2c", "-P 2ab 2n", "-P 2n 2bc",
66
"-P 2ac 2", "-P 2bc 2bc", "-P 2ab 2ab", "-P 2 2ac", "-P 2 2bc",
67
"-P 2ab 2", "-P 2a 2ac", "-P 2b 2c", "-P 2a 2b", "-P 2ac 2c",
68
"-P 2bc 2b", "-P 2b 2ab", "-P 2 2ab", "-P 2bc 2", "-P 2ac 2ac",
69
"-P 2ab 2ac", "-P 2ac 2bc", "-P 2bc 2ab", "-P 2c 2b", "-P 2c 2ac",
70
"-P 2ac 2a", "-P 2b 2a", "-P 2a 2ab", "-P 2bc 2c", "-P 2 2n",
71
"-P 2n 2", "-P 2n 2n", "P 2 2ab -1ab", "-P 2ab 2a", "P 2bc 2 -1bc",
72
"-P 2c 2bc", "P 2ac 2ac -1ac", "-P 2c 2a", "-P 2n 2ab", "-P 2n 2c",
73
"-P 2a 2n", "-P 2bc 2n", "-P 2ac 2b", "-P 2b 2ac", "-P 2ac 2ab",
74
"-P 2bc 2ac", "-P 2ac 2n", "-P 2bc 2a", "-P 2c 2ab", "-P 2n 2ac",
75
"-P 2n 2a", "-P 2c 2n", "-C 2c 2", "-C 2c 2c", "-A 2a 2a",
76
"-A 2 2a", "-B 2 2b", "-B 2b 2", "-C 2bc 2", "-C 2bc 2bc",
77
"-A 2ac 2ac", "-A 2 2ac", "-B 2 2bc", "-B 2bc 2", "-C 2 2",
78
"-A 2 2", "-B 2 2", "-C 2 2c", "-A 2a 2", "-B 2b 2b",
79
"-C 2b 2", "-C 2b 2b", "-A 2c 2c", "-A 2 2c", "-B 2 2c",
80
"-B 2c 2", "C 2 2 -1bc", "-C 2b 2bc", "C 2 2 -1bc", "-C 2b 2c",
81
"A 2 2 -1ac", "-A 2a 2c", "A 2 2 -1ac", "-A 2ac 2c", "B 2 2 -1bc",
82
"-B 2bc 2b", "B 2 2 -1bc", "-B 2b 2bc", "-F 2 2", "F 2 2 -1d",
83
"-F 2uv 2vw", "-I 2 2", "-I 2 2c", "-I 2a 2", "-I 2b 2b",
84
"-I 2b 2c", "-I 2a 2b", "-I 2b 2", "-I 2a 2a", "-I 2c 2c",
85
"-I 2 2b", "-I 2 2a", "-I 2c 2", "P 4", "P 4w",
86
"P 4c", "P 4cw", "I 4", "I 4bw", "P -4",
87
"I -4", "-P 4", "-P 4c", "P 4ab -1ab", "-P 4a",
88
"P 4n -1n", "-P 4bc", "-I 4", "I 4bw -1bw", "-I 4ad",
89
"P 4 2", "P 4ab 2ab", "P 4w 2c", "P 4abw 2nw", "P 4c 2",
90
"P 4n 2n", "P 4cw 2c", "P 4nw 2abw", "I 4 2", "I 4bw 2bw",
91
"P 4 -2", "P 4 -2ab", "P 4c -2c", "P 4n -2n", "P 4 -2c",
92
"P 4 -2n", "P 4c -2", "P 4c -2ab", "I 4 -2", "I 4 -2c",
93
"I 4bw -2", "I 4bw -2c", "P -4 2", "P -4 2c", "P -4 2ab",
94
"P -4 2n", "P -4 -2", "P -4 -2c", "P -4 -2ab", "P -4 -2n",
95
"I -4 -2", "I -4 -2c", "I -4 2", "I -4 2bw", "-P 4 2",
96
"-P 4 2c", "P 4 2 -1ab", "-P 4a 2b", "P 4 2 -1n", "-P 4a 2bc",
97
"-P 4 2ab", "-P 4 2n", "P 4ab 2ab -1ab", "-P 4a 2a", "P 4ab 2n -1ab",
98
"-P 4a 2ac", "-P 4c 2", "-P 4c 2c", "P 4n 2c -1n", "-P 4ac 2b",
99
"P 4n 2 -1n", "-P 4ac 2bc", "-P 4c 2ab", "-P 4n 2n", "P 4n 2n -1n",
100
"-P 4ac 2a", "P 4n 2ab -1n", "-P 4ac 2ac", "-I 4 2", "-I 4 2c",
101
"I 4bw 2bw -1bw", "-I 4bd 2", "I 4bw 2aw -1bw", "-I 4bd 2c", "P 3",
102
"P 31", "P 32", "R 3", "P 3*", "-P 3",
103
"-R 3", "-P 3*", "P 3 2", "P 3 2\"", "P 31 2c (0 0 1)",
104
"P 31 2\"", "P 32 2c (0 0 -1)", "P 32 2\"", "R 3 2\"", "P 3* 2",
105
"P 3 -2\"", "P 3 -2", "P 3 -2\"c", "P 3 -2c", "R 3 -2\"",
106
"P 3* -2", "R 3 -2\"c", "P 3* -2n", "-P 3 2", "-P 3 2c",
107
"-P 3 2\"", "-P 3 2\"c", "-R 3 2\"", "-P 3* 2", "-R 3 2\"c",
108
"-P 3* 2n", "P 6", "P 61", "P 65", "P 62",
109
"P 64", "P 6c", "P -6", "-P 6", "-P 6c",
110
"P 6 2", "P 61 2 (0 0 -1)", "P 65 2 (0 0 1)", "P 62 2c (0 0 1)", "P 64 2c (0 0 -1)",
111
"P 6c 2c", "P 6 -2", "P 6 -2c", "P 6c -2", "P 6c -2c",
112
"P -6 2", "P -6c 2", "P -6 -2", "P -6c -2c", "-P 6 2",
113
"-P 6 2c", "-P 6c 2", "-P 6c 2c", "P 2 2 3", "F 2 2 3",
114
"I 2 2 3", "P 2ac 2ab 3", "I 2b 2c 3", "-P 2 2 3", "P 2 2 3 -1n",
115
"-P 2ab 2bc 3", "-F 2 2 3", "F 2 2 3 -1d", "-F 2uv 2vw 3", "-I 2 2 3",
116
"-P 2ac 2ab 3", "-I 2b 2c 3", "P 4 2 3", "P 4n 2 3", "F 4 2 3",
117
"F 4d 2 3", "I 4 2 3", "P 4acd 2ab 3", "P 4bd 2ab 3", "I 4bd 2c 3",
118
"P -4 2 3", "F -4 2 3", "I -4 2 3", "P -4n 2 3", "F -4c 2 3",
119
"I -4bd 2c 3", "-P 4 2 3", "P 4 2 3 -1n", "-P 4a 2bc 3", "-P 4n 2 3",
120
"P 4n 2 3 -1n", "-P 4bc 2bc 3", "-F 4 2 3", "-F 4c 2 3", "F 4d 2 3 -1d",
121
"-F 4vw 2vw 3", "F 4d 2 3 -1cd", "-F 4cvw 2vw 3", "-I 4 2 3", "-I 4bd 2c 3",
122
// the 51 pathological point groups
123
"-P 1", "-P 2", "-P 2y", "-P 2x", "-P 2\"", "-P 2y\"", "-P 2x\"", "-P 2'", "-P 2y'", "-P 2x'", "-P 2 2", "-P 2 2\"", "-P 2 2\"(y,z,x)", "-P 2 2\"(z,x,y)", "-P 3", "-P 3 (y,z,x)", "-P 3 (z,x,y)", "-P 3 (-x,y,z)", "-P 3 (y,z,-x)", "-P 3 (z,-x,y)", "-P 3*", "-P 3* (-x,y,z)", "-P 3* (x,-y,z)", "-P 3* (x,y,-z)", "-P 3 2", "-P 3 2 (y,z,x)", "-P 3 2 (z,x,y)", "-P 3* 2", "-P 3* 2 (-x,y,z)", "-P 3* 2 (x,-y,z)", "-P 3* 2 (-x,-y,z)", "-P 3 2\"", "-P 3 2\"(z,x,y)", "-P 3 2\"(y,z,x)", "-P 3 2\"(-x,y,z)", "-P 3 2\"(z,-x,y)", "-P 3 2\"(y,z,-x)", "-P 4", "-P 4 (y,z,x)", "-P 4 (z,x,y)", "-P 4 2", "-P 4 2 (y,z,x)", "-P 4 2 (z,x,y)", "-P 6", "-P 6 (y,z,x)", "-P 6 (z,x,y)", "-P 6 2", "-P 6 2 (y,z,x)", "-P 6 2 (z,x,y)", "-P 2 2 3", "-P 4 2 3",
128
Cell c( Cell_descr( 11, 12, 13, 70, 80, 90 ) );
129
std::cout << c.format() << "\n" << CCTBX::cell(CCTBX::cell(c)).format() << "\n";
133
std::cout << hkl.format() << "\t" << CCTBX::hkl(CCTBX::Hkl(hkl)).format() << "\n";
134
std::cout << hkl.format() << "\t" << CCTBX::Hkl(CCTBX::hkl(hkl)).format() << "\n";
136
// test fft in each spacegroup
137
Cell cellc( Cell_descr( 37, 37, 37 ) );
138
Cell cellha( Cell_descr( 37, 37, 37, 120, 90, 90 ) );
139
Cell cellhb( Cell_descr( 37, 37, 37, 90, 120, 90 ) );
140
Cell cellhc( Cell_descr( 37, 37, 37, 90, 90, 120 ) );
141
Cell cellha1( Cell_descr( 37, 37, 37, 60, 90, 90 ) );
142
Cell cellhb1( Cell_descr( 37, 37, 37, 90, 60, 90 ) );
143
Cell cellhc1( Cell_descr( 37, 37, 37, 90, 90, 60 ) );
145
for ( int s = 0; s < sizeof(hallsymbols)/sizeof(hallsymbols[0]); s++ ) {
147
String symbol = hallsymbols[s];
148
cctbx::sgtbx::SpaceGroup sgops(hallsymbols[s]);
149
Spacegroup s1( Spgr_descr( symbol, Spgr_descr::Hall ) );
150
Spacegroup s2 = CCTBX::spacegroup( sgops );
151
Spacegroup s3 = CCTBX::spacegroup( CCTBX::spacegroup( s2 ) );
153
// identify trigonal/hexagonal groups
155
for ( int sym = 1; sym < s1.num_symops(); sym++ ) {
156
if ( ( s1.symop(sym).rot()(1,1) * s1.symop(sym).rot()(1,2) == -1 ) ||
157
( s1.symop(sym).rot()(2,1) * s1.symop(sym).rot()(2,2) == -1 ) )
159
if ( ( s1.symop(sym).rot()(0,0) * s1.symop(sym).rot()(0,2) == -1 ) ||
160
( s1.symop(sym).rot()(2,0) * s1.symop(sym).rot()(2,2) == -1 ) )
162
if ( ( s1.symop(sym).rot()(0,0) * s1.symop(sym).rot()(0,1) == -1 ) ||
163
( s1.symop(sym).rot()(1,0) * s1.symop(sym).rot()(1,1) == -1 ) )
165
if ( ( s1.symop(sym).rot()(1,1) * s1.symop(sym).rot()(1,2) == 1 ) ||
166
( s1.symop(sym).rot()(2,1) * s1.symop(sym).rot()(2,2) == 1 ) )
168
if ( ( s1.symop(sym).rot()(0,0) * s1.symop(sym).rot()(0,2) == 1 ) ||
169
( s1.symop(sym).rot()(2,0) * s1.symop(sym).rot()(2,2) == 1 ) )
171
if ( ( s1.symop(sym).rot()(0,0) * s1.symop(sym).rot()(0,1) == 1 ) ||
172
( s1.symop(sym).rot()(1,0) * s1.symop(sym).rot()(1,1) == 1 ) )
178
Atom atom = Atom::null();
179
atom.set_occupancy( 1.0 );
180
atom.set_u_iso( 0.5 );
181
atom.set_element( "C" );
182
atom.set_coord_orth( Coord_orth( 12, 8, 5 ) );
183
atoms.push_back( atom );
184
atom.set_element( "N" );
185
atom.set_coord_orth( Coord_orth( 11, 6, 4 ) );
186
atoms.push_back( atom );
187
atom.set_element( "O" );
188
atom.set_coord_orth( Coord_orth( 13, 5, 4 ) );
189
atoms.push_back( atom );
193
MMDB mmdb( s1, cell );
194
clipper::NDBModel model;
195
clipper::NDBChain chain;
196
clipper::NDBResidue residue;
197
clipper::NDBAtom atom;
200
residue.set_type("GLY");
201
atom = clipper::NDBAtom::null();
202
atom.set_occupancy(1.0);
204
clipper::DBModel m1 = mmdb.add_model( model );
205
clipper::DBChain c1 = m1.add_chain( chain );
206
clipper::DBResidue r1 = c1.add_residue( residue );
207
atom.set_element( "C" );
208
atom.set_coord_orth( Coord_orth( 12, 8, 5 ) );
210
atom.set_element( "N" );
211
atom.set_coord_orth( Coord_orth( 11, 6, 4 ) );
213
atom.set_element( "O" );
214
atom.set_coord_orth( Coord_orth( 13, 5, 4 ) );
216
mmdb.finalise_edit();
217
clipper::DBAtom_selection atoms = mmdb.select_atoms_serial( 0, 0 );
220
// calc structure factors by slow fft
221
HKL_info hkl_info( s1, cell, Resolution( 3.0 ), true );
222
HKL_data<data32::F_phi> f_phi1( hkl_info );
223
HKL_data<data32::F_phi> f_phi2( hkl_info );
224
SFcalc_iso_sum<float>( f_phi1, atoms );
225
SFcalc_iso_fft<float>( f_phi2, atoms );
226
Grid_sampling grid( s1, cell, Resolution( 3.0 ) );
227
Xmap<float> xmap( s1, cell, grid );
228
FFTmap fftmap( s1, cell, grid );
229
xmap.fft_from( f_phi2, Xmap_base::Normal );
230
xmap.fft_to( f_phi2, Xmap_base::Sparse );
231
fftmap.fft_rfl_to_map( f_phi2, xmap );
232
xmap.fft_to( f_phi2, Xmap_base::Normal );
233
xmap.fft_from( f_phi2, Xmap_base::Sparse );
234
fftmap.fft_map_to_rfl( xmap, f_phi2 );
237
for (int h=-2; h<=2; h++)
238
for (int k=-2; k<=2; k++)
239
for (int l=-2; l<=2; l++) {
241
if ( !HKL_class( s1, rfl ).sys_abs() )
242
if ( std::abs( std::complex<float>(fftmap.get_recip_data(rfl) ) -
243
std::complex<float>(f_phi2[rfl]) ) > 1.0e-3 ) {
245
std::cout << rfl.format() <<
246
"\t: " << fftmap.get_recip_data(rfl).f() <<
247
"\t" << fftmap.get_recip_data(rfl).phi() <<
248
"\t: " << f_phi2[rfl].f() <<
249
"\t" << f_phi2[rfl].phi() <<
250
"\t" << f_phi2[rfl].missing() << ":" << HKL_class( s1, rfl ).sys_abs() << "\n";
252
if ( HKL_class( s1, rfl ).sys_abs() != f_phi2[rfl].missing() ) {
254
int sym; bool friedel;
255
int ih = hkl_info.index_of( hkl_info.find_sym( rfl, sym, friedel) );
256
std::cout << rfl.format() << ih << " " << sym << " " << friedel <<
257
"\t" << hkl_info.hkl_of( ih ).format() <<
258
"\t" << f_phi2[rfl].missing() <<
259
":" << HKL_class( s1, rfl ).sys_abs() << "\n";
263
HKL_info::HKL_reference_index ih;
264
float tol = 0.002 * f_phi1[ HKL( 0, 0, 0 ) ].f();
265
for ( ih = hkl_info.first(); !ih.last(); ih.next() )
266
if ( std::abs( std::complex<float>(f_phi1[ih]) -
267
std::complex<float>(f_phi2[ih]) ) > tol ) {
269
std::cout << ih.hkl().format() <<
270
"\t: " << f_phi1[ih].f() << "\t" << f_phi1[ih].phi() <<
271
"\t: " << f_phi2[ih].f() << "\t" << f_phi2[ih].phi() <<
276
if ( s1.hash() != s2.hash() || s1.hash() != s3.hash() || nerr > 1 )
277
std::cout << "Fail ";
280
String name = symbol + " ";
281
std::cout << name.substr(0,17) << cell.alpha_deg() << " " << cell.beta_deg() << " " << cell.gamma_deg() << "\t" << s1.asu_max().format() << " " << nerr << "\n";
added
removed
examples/cctbxtest.cpp
