~mok0/clipper/trunk

« back to all changes in this revision

Viewing changes to clipper/core/hkl_info.cpp

  • Committer: Morten Kjeldgaard
  • Date: 2011-02-09 21:25:56 UTC
  • Revision ID: mok@bioxray.dk-20110209212556-3zbon5ukh4zf3x00
Tags: 2.0.3
versionĀ 2.0.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* hkl_info.cpp: class file for hkl_info class + children */
 
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 "hkl_info.h"
 
44
 
 
45
 
 
46
namespace clipper {
 
47
 
 
48
 
 
49
Message_fatal message_recip_asu_error( "HKL_info: find_sym reciprocal space ASU error" );
 
50
Message_ctor message_ctor_hkl_info( " [HKL_info: constructed]" );
 
51
 
 
52
 
 
53
// methods for 'HKL_info'
 
54
// private methods:
 
55
 
 
56
/*! Update all the lookup tables to be consistent with the modified
 
57
  reflection list */
 
58
void HKL_info::update_hkl_list() {
 
59
  lookup.init( hkl );
 
60
  hkl_class_lookup.resize( num_reflections() );
 
61
  invresolsq_lookup.resize( num_reflections() );
 
62
  invresolsq_range_ = Range<ftype>();
 
63
  for ( int i = 0; i < num_reflections(); i++ ) {
 
64
    hkl_class_lookup[i] = spacegroup_.hkl_class( hkl_of( i ) );
 
65
    invresolsq_lookup[i] = hkl_of( i ).invresolsq( cell_ );
 
66
    invresolsq_range_.include( invresolsq_lookup[i] );
 
67
  }
 
68
}
 
69
 
 
70
 
 
71
// public methods:
 
72
// constructors/destructor
 
73
 
 
74
HKL_info::HKL_info()
 
75
{
 
76
  Message::message( message_ctor_hkl_info );
 
77
}
 
78
 
 
79
 
 
80
/*! Construct and initialise HKL_info object. This updates the
 
81
  spacegroup and cell and clears the reflection list. The resolution
 
82
  is used as a rejection criterion for reflections - no HKL will be
 
83
  stored beyond the given limit. Initially there are no reflections in
 
84
  the reflection list: see generate_hkl_list().
 
85
 
 
86
  If any of the parameters have null values, the existing values will
 
87
  be unchanged. The object will only be fully initialised once all
 
88
  parameters are available.
 
89
  \param spacegroup The spacegroup.
 
90
  \param cell The unit cell.
 
91
  \param resolution The resolution limit. */
 
92
HKL_info::HKL_info( const Spacegroup& spacegroup, const Cell& cell, const Resolution& resolution, const bool& generate )
 
93
{
 
94
  init( spacegroup, cell, resolution, generate );
 
95
  Message::message( message_ctor_hkl_info );
 
96
}
 
97
 
 
98
 
 
99
/*! Initialise the HKL_info object. This updates the spacegroup and
 
100
  cell and clears the reflection list. The resolution is used as a
 
101
  rejection criterion for reflections - no HKL will be stored beyond
 
102
  the given limit. Initially there are no reflections in the
 
103
  reflection list: see generate_hkl_list().
 
104
 
 
105
  If any of the parameters have null values, the existing values will
 
106
  be unchanged. The object will only be fully initialised once all
 
107
  parameters are available.
 
108
  \param spacegroup The spacegroup.
 
109
  \param cell The unit cell.
 
110
  \param resolution The resolution limit.
 
111
  \param generate If true, a reflection list will be generated for an ASU. */
 
112
void HKL_info::init( const Spacegroup& spacegroup, const Cell& cell, const Resolution& resolution, const bool& generate )
 
113
{
 
114
  // set spacegroup and cell
 
115
  spacegroup_ = spacegroup;
 
116
  cell_ = cell;
 
117
  resolution_ = resolution;
 
118
 
 
119
  // check parameters
 
120
  if ( is_null() ) return;
 
121
 
 
122
  // Create the intergised symops (grid irrelevent)
 
123
  Grid g( 24, 24, 24 );
 
124
  isymop.resize( spacegroup_.num_symops() );
 
125
  for ( int sym = 0; sym < spacegroup_.num_symops(); sym++ )
 
126
    isymop[sym] = Isymop( spacegroup_.symop(sym), g );
 
127
 
 
128
  // reflection lists
 
129
  hkl.clear();
 
130
  if ( generate ) generate_hkl_list();
 
131
  update_hkl_list();
 
132
}
 
133
 
 
134
 
 
135
/*! Initialise the HKL_info object. This updates the spacegroup and
 
136
  cell and clears the reflection list. The HKL_sampling determines
 
137
  the reflection list.
 
138
 
 
139
  If any of the parameters have null values, the existing values will
 
140
  be unchanged. The object will only be fully initialised once all
 
141
  parameters are available.
 
142
  \param spacegroup The spacegroup.
 
143
  \param cell The unit cell.
 
144
  \param hkl_sampling The resolution limit.
 
145
  \param generate If true, a reflection list will be generated for an ASU. */
 
146
void HKL_info::init( const Spacegroup& spacegroup, const Cell& cell, const HKL_sampling& hkl_sampling, const bool& generate )
 
147
{
 
148
  // set spacegroup and cell
 
149
  spacegroup_ = spacegroup;
 
150
  cell_ = cell;
 
151
  hkl_sampling_ = hkl_sampling;
 
152
 
 
153
  // check parameters
 
154
  if ( spacegroup_.is_null() || cell_.is_null() || hkl_sampling_.is_null() ) return;
 
155
 
 
156
  // estimate resolution
 
157
  resolution_ = hkl_sampling.resolution( cell );
 
158
 
 
159
  // Create the intergised symops (grid irrelevent)
 
160
  Grid g( 24, 24, 24 );
 
161
  isymop.resize( spacegroup_.num_symops() );
 
162
  for ( int sym = 0; sym < spacegroup_.num_symops(); sym++ )
 
163
    isymop[sym] = Isymop( spacegroup_.symop(sym), g );
 
164
 
 
165
  // reflection lists
 
166
  hkl.clear();
 
167
  if ( generate ) {
 
168
    std::vector<HKL> add;
 
169
    // make a reflection list for the given refln sampling
 
170
    HKL lim = hkl_sampling.hkl_limit();
 
171
    HKL rfl;
 
172
    // make a list of valid reflections
 
173
    for (rfl.h()=-lim.h(); rfl.h()<=lim.h(); rfl.h()++)
 
174
      for (rfl.k()=-lim.k(); rfl.k()<=lim.k(); rfl.k()++)
 
175
        for (rfl.l()=-lim.l(); rfl.l()<=lim.l(); rfl.l()++)
 
176
          if ( spacegroup_.recip_asu(rfl) &&
 
177
               hkl_sampling.in_resolution( rfl ) &&
 
178
               !( spacegroup_.hkl_class(rfl).sys_abs() ) ) add.push_back(rfl);
 
179
    // update the reflection data lists
 
180
    add_hkl_list( add );
 
181
  }
 
182
 
 
183
  update_hkl_list();
 
184
}
 
185
 
 
186
 
 
187
/*! \return true if the object has not been initalised. */
 
188
bool HKL_info::is_null() const
 
189
{ return ( spacegroup_.is_null() || cell_.is_null() || resolution_.is_null() ); }
 
190
 
 
191
 
 
192
/*! Using current cell, spacegroup, resolution. */
 
193
void HKL_info::generate_hkl_list()
 
194
{
 
195
  std::vector<HKL> add;
 
196
  // make a reflection list for the given cell, symm and resolution
 
197
  /* TO MAKE A BOX TO HOLD A SPHERE IN REAL OR RECIPROCAL SPACE
 
198
     In reciprocal space, use the real space cell dimension / resolution
 
199
     In real space, use the recip space dimension * radius */
 
200
  HKL rfl;
 
201
  int hmax = int(cell_.descr().a()/resolution_.limit());
 
202
  int kmax = int(cell_.descr().b()/resolution_.limit());
 
203
  int lmax = int(cell_.descr().c()/resolution_.limit());
 
204
  ftype s_lim = resolution_.invresolsq_limit();
 
205
  // make a list of valid reflections
 
206
  for (rfl.h()=-hmax; rfl.h()<=hmax; rfl.h()++)
 
207
    for (rfl.k()=-kmax; rfl.k()<=kmax; rfl.k()++)
 
208
      for (rfl.l()=-lmax; rfl.l()<=lmax; rfl.l()++)
 
209
        if ( spacegroup_.recip_asu(rfl) && rfl.invresolsq(cell_) < s_lim &&
 
210
             !( spacegroup_.hkl_class(rfl).sys_abs() ) ) add.push_back(rfl);
 
211
  // update the reflection data lists
 
212
  hkl.clear();
 
213
  add_hkl_list( add );
 
214
}
 
215
 
 
216
 
 
217
/*! The new HKLs are transformed to the default reciprocal ASU, and
 
218
  added to the reflection list. Duplicates and reflections outside the
 
219
  resoluution limit are ignored. Then the fast lookup tables for HKL,
 
220
  invresolsq, and reflection class are rebuilt.
 
221
  \param add The list of new reflections to add. */
 
222
void HKL_info::add_hkl_list( const std::vector<HKL>& add ) {
 
223
  HKL equiv; int sym; bool friedel;
 
224
  for ( int i = 0; i < add.size(); i++ ) {
 
225
    if ( add[i].invresolsq( cell_ ) <= resolution_.invresolsq_limit() ) {
 
226
      equiv = find_sym( add[i], sym, friedel );
 
227
      if ( lookup.index_of( equiv ) < 0 ) hkl.push_back( equiv );
 
228
    }
 
229
  }
 
230
  update_hkl_list();
 
231
}
 
232
 
 
233
 
 
234
/*! Returns the index of the reflection, the sym no. and Friedel flag.
 
235
 \internal */
 
236
HKL HKL_info::find_sym( const HKL& rfl, int& sym, bool& friedel ) const
 
237
{
 
238
  // find the symmetry operator mapping hkl into the reflection
 
239
  // list and determine the friedel flag
 
240
  HKL equiv;
 
241
 
 
242
  // now find the symop which gives a reflection from the list
 
243
  for (sym = 0; sym < spacegroup_.num_primops(); sym++) {
 
244
    equiv = rfl.transform(isymop[sym]);
 
245
    if ( spacegroup_.recip_asu(equiv) ) { friedel = false; return equiv; }
 
246
    equiv = -equiv;
 
247
    if ( spacegroup_.recip_asu(equiv) ) { friedel = true ; return equiv; }
 
248
  }
 
249
  Message::message( message_recip_asu_error );
 
250
  return equiv;
 
251
}
 
252
 
 
253
 
 
254
void HKL_info::debug() const
 
255
{
 
256
  std::cout << "Num reflns " << hkl.size() << "\n";
 
257
}
 
258
 
 
259
 
 
260
} // namespace clipper