1
/* cns_hkl_io.cpp: class file for reflection data cns_hkl importer */
2
//C Copyright (C) 2000-2006 Kevin Cowtan and University of York
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:
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.)'
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.
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
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.
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,
42
#include "cns_hkl_io.h"
52
clipper::String cnstok( FILE* f ) {
55
while ( c > '\0' && c < ' ' ) c = fgetc( f );
56
while ( c > ' ' && c != '=') {
64
/*! Constructing an CNS_HKLfile does nothing except flag the object as not
65
attached to any file for either input or output */
66
CNS_HKLfile::CNS_HKLfile()
72
/*! Close any files which were left open. This is particularly
73
important since to access the CNS_HKL file efficiently, data reads and
74
writes are deferred until the file is closed. */
75
CNS_HKLfile::~CNS_HKLfile()
88
/*! The file is opened for reading. This CNS_HKLfile object will remain
89
attached to this file until it is closed. Until that occurs, no
90
other file may be opened with this object, however another CNS_HKLfile
91
object could be used to access another file.
92
\param filename_in The input filename or pathname. */
93
void CNS_HKLfile::open_read( const String filename_in )
96
Message::message( Message_fatal( "CNS_HKLfile: open_read - File already open" ) );
99
f_sigf_i = NULL; phi_wt_i = NULL; f_phi_i = NULL; abcd_i = NULL; flag_i = NULL;
100
filename = filename_in;
102
FILE* cns_hkl = fopen( filename.c_str(), "r" );
103
if ( cns_hkl == NULL )
104
Message::message( Message_fatal( "CNS_HKLfile: open_read - Could not read: "+filename ) );
111
/*! Close the file after reading. This command also actually fills in
112
the data in any HKL_data structures which have been marked for
114
void CNS_HKLfile::close_read()
117
Message::message( Message_fatal( "CNS_HKLfile: no file open for read" ) );
119
// make sure the data list is sized
120
if ( f_sigf_i != NULL ) f_sigf_i->update();
121
if ( phi_wt_i != NULL ) phi_wt_i->update();
122
if ( f_phi_i != NULL ) f_phi_i->update();
123
if ( abcd_i != NULL ) abcd_i->update();
124
if ( flag_i != NULL ) flag_i->update();
128
xtype fo[8], pw[8], fc[8], hl[8], fl[8];
132
FILE* cns_hkl = fopen( filename.c_str(), "r" );
133
if ( cns_hkl == NULL )
134
Message::message( Message_fatal( "CNS_HKLfile: import_hkl_data - Could not read: "+filename ) );
135
// read the data from the CNS_HKL
138
while ( (s=cnstok(cns_hkl)).length() > 0 ) {
139
t = (s+" ").substr(0,4);
142
hkl = HKL( h, k, l );
143
if ( f_sigf_i != NULL ) f_sigf_i->data_import( hkl, fo );
144
if ( phi_wt_i != NULL ) phi_wt_i->data_import( hkl, pw );
145
if ( f_phi_i != NULL ) f_phi_i->data_import( hkl, fc );
146
if ( abcd_i != NULL ) abcd_i->data_import( hkl, hl );
147
if ( flag_i != NULL ) flag_i->data_import( hkl, fl );
150
for ( int i = 0; i < 8; i++ ) fo[i] = pw[i] = fc[i] = hl[i] = fl[i] = 0.0;
152
if ( t == "INDE" ) code = 10;
153
else if ( t == "FOBS" ) code = 20;
154
else if ( t == "SIGM" ) code = 30;
155
else if ( t == "FOM " ) code = 40;
156
else if ( t == "FCAL" ) code = 50;
157
else if ( t == "ABCD" ) code = 60;
158
else if ( t == "TEST" ) code = 70;
159
else if ( code == 10 ) { h = s.i(); code++; }
160
else if ( code == 11 ) { k = s.i(); code++; }
161
else if ( code == 12 ) { l = s.i(); code++; }
162
else if ( code == 20 ) { fo[0] = s.f(); code++; }
163
else if ( code == 21 ) { pw[0] = s.f(); code++; }
164
else if ( code == 30 ) { fo[1] = s.f(); code++; }
165
else if ( code == 40 ) { pw[1] = s.f(); code++; }
166
else if ( code == 50 ) { fc[0] = s.f(); code++; }
167
else if ( code == 51 ) { fc[1] = s.f(); code++; }
168
else if ( code == 60 ) { hl[0] = s.f(); code++; }
169
else if ( code == 61 ) { hl[1] = s.f(); code++; }
170
else if ( code == 62 ) { hl[2] = s.f(); code++; }
171
else if ( code == 63 ) { hl[3] = s.f(); code++; }
172
else if ( code == 70 ) { fl[0] = s.f(); code++; }
174
if (h||k||l) { // get last rfl
175
hkl = HKL( h, k, l );
176
if ( f_sigf_i != NULL ) f_sigf_i->data_import( hkl, fo );
177
if ( phi_wt_i != NULL ) phi_wt_i->data_import( hkl, pw );
178
if ( f_phi_i != NULL ) f_phi_i->data_import( hkl, fc );
179
if ( abcd_i != NULL ) abcd_i->data_import( hkl, hl );
180
if ( flag_i != NULL ) flag_i->data_import( hkl, fl );
188
/*! The file is opened for writing. This will be a new file, created
189
entirely from data from within the program, rather than by extending
190
an existing file. Similar restrictions apply as for open_read().
192
\param filename_out The output filename or pathname. */
193
void CNS_HKLfile::open_write( const String filename_out )
196
Message::message( Message_fatal( "CNS_HKLfile: open_write - File already open" ) );
198
// open the output cns_hkl
200
f_sigf_o = NULL; phi_wt_o = NULL; f_phi_o = NULL; abcd_o = NULL; flag_o = NULL;
201
filename = filename_out;
203
FILE* cns_hkl = fopen( filename.c_str(), "w" );
204
if ( cns_hkl == NULL )
205
Message::message( Message_fatal( "CNS_HKLfile: open_write - Could not write: "+filename ) );
212
/*! Close the file after writing. This command also actually writes
213
the data reflection list from the HKL_info object and the data from
214
any HKL_data objects which have been marked for import. */
215
void CNS_HKLfile::close_write()
218
Message::message( Message_fatal( "CNS_HKLfile: close_write - no file open for write" ) );
220
// export the marked list data to an cns_hkl file
221
if ( hkl_ptr == NULL )
222
Message::message( Message_fatal( "CNS_HKLfile: close_write - no refln list exported" ) );
223
const HKL_info& hklinf = *hkl_ptr;
227
float f1, f2, f3, f4;
229
FILE* cns_hkl = fopen( filename.c_str(), "w" );
230
if ( cns_hkl == NULL )
231
Message::message( Message_fatal( "CNS_HKLfile: close_write - Could not write: "+filename ) );
232
fprintf( cns_hkl, "NREF=%i\n", hklinf.num_reflections() );
233
HKL_info::HKL_reference_index ih;
234
for ( ih = hklinf.first(); !ih.last(); ih.next() ) {
235
f1 = f2 = f3 = f4 = 0.0; // default sigf to 0 in case missing
237
fprintf( cns_hkl, "INDE %i %i %i", hkl.h(), hkl.k(), hkl.l() );
238
if ( f_sigf_o != NULL ) {
239
f_sigf_o->data_export( hkl, x );
240
for ( int i = 0; i < 4; i++ ) if ( Util::is_nan(x[i]) ) x[i] = 0.0;
244
if ( phi_wt_o != NULL ) {
245
phi_wt_o->data_export( hkl, x );
246
for ( int i = 0; i < 4; i++ ) if ( Util::is_nan(x[i]) ) x[i] = 0.0;
250
fprintf( cns_hkl, " FOBS=%.3f %.3f SIGM=%.3f FOM=%.3f",f1,f3,f2,f4 );
251
if ( f_phi_o != NULL ) {
252
f_phi_o->data_export( hkl, x );
253
for ( int i = 0; i < 4; i++ ) if ( Util::is_nan(x[i]) ) x[i] = 0.0;
254
fprintf( cns_hkl, " FCAL=%.3f %.3f",float(x[0]),float(x[1]) );
256
if ( abcd_o != NULL ) {
257
abcd_o->data_export( hkl, x );
258
for ( int i = 0; i < 4; i++ ) if ( Util::is_nan(x[i]) ) x[i] = 0.0;
259
fprintf( cns_hkl, " HLA=%.1f HLB=%.1f HLC=%.1f HLD=%.1f",float(x[0]),float(x[1]),float(x[2]),float(x[3]) );
261
if ( flag_o != NULL ) {
262
abcd_o->data_export( hkl, x );
263
for ( int i = 0; i < 4; i++ ) if ( Util::is_nan(x[i]) ) x[i] = 0.0;
264
fprintf( cns_hkl, " TEST=%i",Util::intr(x[0]) );
266
fprintf( cns_hkl, "\n" );
274
/*! Get the resolution limit from the CNS_HKL file.
275
Since a CNS_HKL file does not contain cell information, a Cell object
276
must be supplied, which will be used to determine the resultion.
277
The result is the resolution determined by the most extreme
278
reflection in the file.
279
\return The resolution. */
280
Resolution CNS_HKLfile::resolution( const Cell& cell ) const
283
Message::message( Message_fatal( "CNS_HKLfile: resolution - no file open for read" ) );
287
FILE* cns_hkl = fopen( filename.c_str(), "r" );
288
if ( cns_hkl == NULL )
289
Message::message( Message_fatal( "CNS_HKLfile: resolution - Could not read: "+filename ) );
290
// read the reflections from the cns_hkl
293
while ( (s=cnstok(cns_hkl)).length() > 0 ) {
294
t = (s+" ").substr(0,4);
296
h = cnstok(cns_hkl).i();
297
k = cnstok(cns_hkl).i();
298
l = cnstok(cns_hkl).i();
299
hkl = HKL( h, k, l );
300
slim = Util::max( slim, hkl.invresolsq(cell) );
305
return Resolution( 1.0/sqrt(slim) );
309
/*! Import the list of reflection HKLs from an CNS_HKL file into an
310
HKL_info object. If the resolution limit of the HKL_info object is
311
lower than the limit of the file, any excess reflections will be
312
rejected, as will any systematic absences or duplicates.
313
\param target The HKL_info object to be initialised. */
314
void CNS_HKLfile::import_hkl_info( HKL_info& target )
317
Message::message( Message_fatal( "CNS_HKLfile: import_hkl_info - no file open for read" ) );
319
std::vector<HKL> hkls;
322
FILE* cns_hkl = fopen( filename.c_str(), "r" );
323
if ( cns_hkl == NULL )
324
Message::message( Message_fatal( "CNS_HKLfile: import_hkl_info - Could not read: "+filename ) );
325
// read the reflections from the cns_hkl
326
ftype slim = target.resolution().invresolsq_limit();
328
while ( (s=cnstok(cns_hkl)).length() > 0 ) {
329
t = (s+" ").substr(0,4);
331
h = cnstok(cns_hkl).i();
332
k = cnstok(cns_hkl).i();
333
l = cnstok(cns_hkl).i();
334
hkl = HKL( h, k, l );
335
if ( hkl.invresolsq(target.cell()) < slim ) hkls.push_back( hkl );
340
target.add_hkl_list( hkls );
344
/*! Import data from an CNS_HKL file into an HKL_data object.
346
This routine does not actually read any data, but rather marks the
347
data to be read when the file is closed.
349
The data to be read (F_sigF or Phi_fom) will be selected based on
350
the type of the HKL_data object.
352
\param cdata The HKL_data object into which data is to be imported. */
353
void CNS_HKLfile::import_hkl_data( HKL_data_base& cdata )
356
Message::message( Message_fatal( "CNS_HKLfile: import_hkl_data - no file open for read" ) );
358
if ( cdata.type() == data32::F_sigF::type() ) f_sigf_i = &cdata;
359
else if ( cdata.type() == data32::Phi_fom::type() ) phi_wt_i = &cdata;
360
else if ( cdata.type() == data32::F_phi::type() ) f_phi_i = &cdata;
361
else if ( cdata.type() == data32::ABCD::type() ) abcd_i = &cdata;
362
else if ( cdata.type() == data32::Flag::type() ) flag_i = &cdata;
364
Message::message( Message_fatal( "CNS_HKLfile: import_hkl_data - data must be F_sigF/Phi_fom/F_phi/ABCD/Flag" ) );
368
/*! Export the list of reflection HKLs to an CNS_HKL file from an
370
void CNS_HKLfile::export_hkl_info( const HKL_info& target )
373
Message::message( Message_fatal( "CNS_HKLfile: export_hkl_info - no file open for write" ) );
378
/*! Export data from an HKL_data object into an CNS_HKL file.
380
This routine does not actually write any data, but rather marks the
381
data to be written when the file is closed.
383
The data to be read (F_sigF or Phi_fom) will be selected based on
384
the type of the HKL_data object.
386
\param cdata The HKL_data object from which data is to be exported. */
387
void CNS_HKLfile::export_hkl_data( const HKL_data_base& cdata )
390
Message::message( Message_fatal( "CNS_HKLfile: export_hkl_data - no file open for write" ) );
392
if ( cdata.type() == data32::F_sigF::type() ) f_sigf_o = &cdata;
393
else if ( cdata.type() == data32::Phi_fom::type() ) phi_wt_o = &cdata;
394
else if ( cdata.type() == data32::F_phi::type() ) f_phi_o = &cdata;
395
else if ( cdata.type() == data32::ABCD::type() ) abcd_o = &cdata;
396
else if ( cdata.type() == data32::Flag::type() ) flag_o = &cdata;
398
Message::message( Message_fatal( "CNS_HKLfile: export_hkl_data - data must be F_sigF/Phi_fom/F_phi/ABCD/Flag" ) );
402
} // namespace clipper
added
removed
clipper/cns/cns_hkl_io.cpp
