00001
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #ifndef CLIPPER_XMAP
00046 #define CLIPPER_XMAP
00047
00048
00049 #include "fftmap.h"
00050 #include "fftmap_sparse.h"
00051 #include "derivs.h"
00052
00053
00054 namespace clipper
00055 {
00056
00057 class Xmap_cacheobj
00058 {
00059 public:
00060 class Key
00061 {
00062 public:
00063 Key( const Spgr_descr& spgr_descr, const Grid_sampling& grid ) :
00064 spgr_descr_(spgr_descr), grid_sampling_(grid) {}
00065 const Spgr_descr& spgr_descr() const { return spgr_descr_; }
00066 const Grid_sampling& grid_sampling() const { return grid_sampling_; }
00067 private:
00068 Spgr_descr spgr_descr_;
00069 Grid_sampling grid_sampling_;
00070 };
00071
00072 Xmap_cacheobj( const Key& xmap_cachekey );
00073 bool matches( const Key& xmap_cachekey ) const;
00074 String format() const;
00075
00076 Key key;
00077 Grid_sampling xtl_grid;
00078 Grid_range asu_grid;
00079 Grid_range map_grid;
00080 int nsym;
00081 std::vector<unsigned char> asu;
00082 std::vector<Isymop> isymop;
00083 std::vector<int> du, dv, dw;
00084 Array2d<unsigned char> symperm;
00085 Mat33<> mat_grid_orth;
00086 static Mutex mutex;
00087 };
00088
00089
00091
00100 class Xmap_base
00101 {
00102 public:
00103 enum FFTtype { Default, Normal, Sparse };
00104
00106 bool is_null() const;
00107
00109 const Cell& cell() const { return cell_; }
00111 const Spacegroup& spacegroup() const { return spacegroup_; }
00113 const Grid_sampling& grid_sampling() const { return grid_sam_; }
00115 const Grid_range& grid_asu() const { return cacheref.data().asu_grid; }
00117
00118 inline Coord_grid coord_of( const int& index ) const
00119 { return cacheref.data().map_grid.deindex( index ); }
00121
00123 inline int index_of( const Coord_grid& coord ) const {
00124 if ( cacheref.data().asu_grid.in_grid( coord ) ) {
00125 const int i = cacheref.data().map_grid.index( coord );
00126 if ( asu[ i ] == 0 ) return i;
00127 }
00128 return -1;
00129 }
00131 Coord_grid to_map_unit( const Coord_grid& pos ) const
00132 { return pos.unit( grid_sam_ ); }
00133
00135 const RTop<>& operator_orth_grid() const { return rt_orth_grid; }
00137 const RTop<>& operator_grid_orth() const { return rt_grid_orth; }
00139
00141 inline Coord_orth coord_orth( const Coord_map& cm ) const
00142 { return Coord_orth( rt_grid_orth.rot()*cm ); }
00144
00146 inline Coord_map coord_map( const Coord_orth& co ) const
00147 { return Coord_map ( rt_orth_grid.rot()*co ); }
00148
00150 bool in_map( const Coord_grid& pos ) const { return true; }
00152 template<class I> bool in_map( const Coord_map& cm ) const { return true; }
00153
00155 int multiplicity( const Coord_grid& pos ) const;
00156
00158
00162 class Map_reference_base
00163 {
00164 public:
00166 inline const Xmap_base& base_xmap() const { return *map_; }
00168 inline const int& index() const { return index_; }
00170 bool last() const { return ( index_ >= map_->map_grid.size() ); }
00171 protected:
00173 const Xmap_base* map_;
00175 int index_;
00176 };
00177
00179
00190 class Map_reference_index : public Map_reference_base
00191 {
00192 public:
00194 Map_reference_index() {}
00196 explicit Map_reference_index( const Xmap_base& map )
00197 { map_ = ↦ index_=0; next(); }
00199 Map_reference_index( const Xmap_base& map, const Coord_grid& pos ) { map_ = ↦ int sym; map_->find_sym( pos, index_, sym ); }
00201 inline Coord_grid coord() const
00202 { return map_->map_grid.deindex(index_); }
00204 inline const Coord_orth coord_orth() const
00205 { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
00207 inline Map_reference_index& set_coord( const Coord_grid& pos )
00208 { int sym; map_->find_sym( pos, index_, sym ); return *this; }
00210 inline Map_reference_index& next() {
00211 do {
00212 index_++; if ( last() ) break;
00213 } while ( map_->asu[index_] != 0 );
00214 return *this;
00215 }
00217
00218
00219 inline int index_offset(const int& du,const int& dv,const int& dw) const {
00220 int i = index_ + du*map_->du[0] + dv*map_->dv[0] + dw*map_->dw[0];
00221 if ( map_->asu[i] != 0 ) { i = map_->map_grid.index( map_->to_map_unit( map_->map_grid.deindex(i).transform( map_->isymop[map_->asu[i]-1] ) ) ); }
00222 return i;
00223 }
00224
00225
00226
00227
00228 };
00229
00231
00242 class Map_reference_coord : public Map_reference_base
00243 {
00244 public:
00246 Map_reference_coord() {}
00248 explicit Map_reference_coord( const Xmap_base& map )
00249 { map_ = ↦ index_ = 0; next(); }
00251 Map_reference_coord( const Xmap_base& map, const Coord_grid& pos ) {
00252 map_ = ↦
00253 pos_ = pos;
00254 map_->find_sym( pos_, index_, sym_ );
00255 }
00257 inline const Coord_grid& coord() const { return pos_; }
00259 inline const Coord_orth coord_orth() const
00260 { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
00262 inline const int& sym() const { return sym_; }
00264 Map_reference_coord& set_coord( const Coord_grid& pos );
00266
00267 inline Map_reference_coord& next() {
00268 sym_ = 0;
00269 do {
00270 index_++; if ( last() ) break;
00271 } while ( map_->asu[index_] != 0 );
00272 pos_ = map_->map_grid.deindex(index_);
00273 return *this;
00274 }
00275
00276 inline Map_reference_coord& next_u() { pos_.u()++; index_ += map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00277 inline Map_reference_coord& next_v() { pos_.v()++; index_ += map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00278 inline Map_reference_coord& next_w() { pos_.w()++; index_ += map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00279 inline Map_reference_coord& prev_u() { pos_.u()--; index_ -= map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00280 inline Map_reference_coord& prev_v() { pos_.v()--; index_ -= map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00281 inline Map_reference_coord& prev_w() { pos_.w()--; index_ -= map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00282
00283 inline Map_reference_coord& operator =( const Coord_grid& pos )
00284 { return set_coord( pos ); }
00285
00286
00287
00288
00289
00290 protected:
00292 int sym_;
00294 Coord_grid pos_;
00295
00297 void edge();
00298 };
00299
00301 Map_reference_index first() const { return Map_reference_index( *this ); }
00303 Map_reference_coord first_coord() const { return Map_reference_coord( *this ); }
00305 static FFTtype& default_type() { return default_type_; }
00306 protected:
00307 ObjectCache<Xmap_cacheobj>::Reference cacheref;
00308 const unsigned char* asu;
00309 const Isymop* isymop;
00310 const int* du;
00311 const int* dv;
00312 const int* dw;
00313 Grid_range asu_grid;
00314 Grid_range map_grid;
00315 int nsym;
00316
00317 Cell cell_;
00318 Spacegroup spacegroup_;
00319 Grid_sampling grid_sam_;
00320
00321 RTop<> rt_orth_grid;
00322 RTop<> rt_grid_orth;
00323
00325 Xmap_base();
00327 void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam );
00328 inline void find_sym( const Coord_grid& base, int& index, int& sym ) const;
00329 void asu_error( const Coord_grid& pos ) const;
00330
00331 static FFTtype default_type_;
00332
00333 friend class Xmap_base::Map_reference_base;
00334 friend class Xmap_base::Map_reference_index;
00335 friend class Xmap_base::Map_reference_coord;
00336 };
00337
00338
00339
00340
00342
00356 template<class T> class Xmap : public Xmap_base
00357 {
00358 public:
00360 Xmap() {}
00362 Xmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam ) { init( spacegroup, cell, grid_sam ); }
00364 void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam ) { Xmap_base::init( spacegroup, cell, grid_sam ); list.resize( cacheref.data().asu.size() ); }
00365
00367 inline const T& operator[] (const Xmap_base::Map_reference_index& ix) const
00368 { return list[ix.index()]; }
00370 inline T& operator[] (const Xmap_base::Map_reference_index& ix)
00371 { return list[ix.index()]; }
00372
00374 inline const T& operator[] (const Xmap_base::Map_reference_coord& ix) const
00375 { return list[ix.index()]; }
00377 inline T& operator[] (const Xmap_base::Map_reference_coord& ix)
00378 { return list[ix.index()]; }
00379
00381 const T& get_data( const Coord_grid& pos ) const;
00383 void set_data( const Coord_grid& pos, const T& val );
00385 inline const T& get_data( const int& index ) const;
00387 bool set_data( const int& index, const T& val );
00388
00390 template<class I> T interp( const Coord_frac& pos ) const;
00392 template<class I> void interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const;
00394 template<class I> void interp_curv( const Coord_frac& pos, T& val, Grad_frac<T>& grad, Curv_frac<T>& curv ) const;
00396 template<class I> T interp( const Coord_map& pos ) const;
00398 template<class I> void interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const;
00400 template<class I> void interp_curv( const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv ) const;
00401
00403 template<class H> void fft_from( const H& fphidata, const FFTtype type = Default );
00405 template<class H> void fft_to ( H& fphidata, const FFTtype type = Default ) const;
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00420 const T& operator =( const T& value );
00422 const Xmap<T>& operator +=( const Xmap<T>& other );
00424 const Xmap<T>& operator -=( const Xmap<T>& other );
00425
00426 private:
00427 std::vector<T> list;
00428 };
00429
00430
00431
00432
00433
00434 void Xmap_base::find_sym( const Coord_grid& base, int& index, int& sym ) const
00435 {
00436
00437 Coord_grid rot = to_map_unit( base );
00438 if ( asu_grid.in_grid( rot ) ) {
00439 index = map_grid.index( rot );
00440 if ( asu[ index ] == 0 ) {
00441 sym = 0;
00442 } else {
00443 sym = asu[ index ] - 1;
00444 index = map_grid.index( to_map_unit( base.transform(isymop[sym]) ) );
00445 }
00446 } else {
00447
00448 for ( sym = 1; sym < nsym; sym++ ) {
00449 rot = to_map_unit( base.transform( isymop[sym] ) );
00450 if ( asu_grid.in_grid( rot ) ) {
00451 index = map_grid.index( rot );
00452 if ( asu[ index ] == 0 ) return;
00453 }
00454 }
00455 index = 0;
00456 asu_error( base );
00457 }
00458 return;
00459 }
00460
00461
00466 template<class T> const T& Xmap<T>::get_data( const Coord_grid& pos ) const
00467 {
00468 int index, sym;
00469 find_sym( pos, index, sym );
00470 return list[ index ];
00471 }
00472
00477 template<class T> void Xmap<T>::set_data( const Coord_grid& pos, const T& val )
00478 {
00479 int index, sym;
00480 find_sym( pos, index, sym );
00481 list[ index ] = val;
00482 }
00483
00490 template<class T> const T& Xmap<T>::get_data( const int& index ) const
00491 { return list[index]; }
00492
00500 template<class T> bool Xmap<T>::set_data( const int& index, const T& val )
00501 {
00502 if ( index >= 0 && index < list.size() )
00503 if ( asu[index] == 0 ) {
00504 list[index] = val;
00505 return true;
00506 }
00507 return false;
00508 }
00509
00518 template<class T> template<class I> T Xmap<T>::interp( const Coord_frac& pos ) const
00519 {
00520 T val;
00521 I::interp( *this, pos.coord_map( grid_sam_ ), val );
00522 return val;
00523 }
00524
00525
00533 template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const
00534 {
00535 Grad_map<T> g;
00536 I::interp_grad( *this, pos.coord_map( grid_sam_ ), val, g );
00537 grad = g.grad_frac( grid_sam_ );
00538 }
00539
00540
00550 template<class T> template<class I> void Xmap<T>::interp_curv( const Coord_frac& pos, T& val, Grad_frac<T>& grad, Curv_frac<T>& curv ) const
00551 {
00552 Grad_map<T> g;
00553 Curv_map<T> c;
00554 I::interp_curv( *this, pos.coord_map( grid_sam_ ), val, g, c );
00555 grad = g.grad_frac( grid_sam_ );
00556 curv = c.curv_frac( grid_sam_ );
00557 }
00558
00559
00568 template<class T> template<class I> T Xmap<T>::interp( const Coord_map& pos ) const
00569 {
00570 T val;
00571 I::interp( *this, pos, val );
00572 return val;
00573 }
00574
00575
00583 template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const
00584 { I::interp_grad( *this, pos, val, grad ); }
00585
00586
00596 template<class T> template<class I> void Xmap<T>::interp_curv( const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv ) const
00597 { I::interp_curv( *this, pos, val, grad, curv ); }
00598
00599
00604 template<class T> template<class H> void Xmap<T>::fft_from( const H& fphidata, const FFTtype type )
00605 {
00606 if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
00607
00608 FFTmap_sparse_p1_hx fftmap( grid_sampling() );
00609
00610 typename H::HKL_reference_index ih;
00611 ffttype f, phi0, phi1;
00612 int sym;
00613 for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
00614 f = fphidata[ih].f();
00615 if ( f != 0.0 ) {
00616 phi0 = fphidata[ih].phi();
00617 const HKL& hkl = ih.hkl();
00618 fftmap.set_hkl( hkl,
00619 std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
00620 for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
00621 phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
00622 fftmap.set_hkl( hkl.transform( isymop[sym] ),
00623 std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
00624 }
00625 }
00626 }
00627
00628 for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
00629 fftmap.require_real_data( ix.coord() );
00630
00631 fftmap.fft_h_to_x(1.0/cell().volume());
00632
00633 for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
00634 (*this)[ix] = fftmap.real_data( ix.coord() );
00635 } else {
00636
00637 FFTmap_p1 fftmap( grid_sampling() );
00638
00639 typename H::HKL_reference_index ih;
00640 ffttype f, phi0, phi1;
00641 int sym;
00642 for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
00643 f = fphidata[ih].f();
00644 if ( f != 0.0 ) {
00645 phi0 = fphidata[ih].phi();
00646 const HKL& hkl = ih.hkl();
00647 fftmap.set_hkl( hkl,
00648 std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
00649 for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
00650 phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
00651 fftmap.set_hkl( hkl.transform( isymop[sym] ),
00652 std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
00653 }
00654 }
00655 }
00656
00657 fftmap.fft_h_to_x(1.0/cell().volume());
00658
00659 for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
00660 (*this)[ix] = fftmap.real_data( ix.coord() );
00661 }
00662 }
00663
00664
00673 template<class T> template<class H> void Xmap<T>::fft_to ( H& fphidata, const FFTtype type ) const
00674 {
00675 if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
00676
00677 FFTmap_sparse_p1_xh fftmap( grid_sampling() );
00678
00679 ffttype f;
00680 int sym;
00681 for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
00682 f = (*this)[ix];
00683 if ( f != 0.0 ) {
00684 fftmap.real_data( ix.coord() ) = f;
00685 for ( sym = 1; sym < cacheref.data().nsym; sym++ )
00686 fftmap.real_data(
00687 ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
00688 }
00689 }
00690
00691 typename H::HKL_reference_index ih;
00692 for ( ih = fphidata.first(); !ih.last(); ih.next() )
00693 fftmap.require_hkl( ih.hkl() );
00694
00695 fftmap.fft_x_to_h(cell().volume());
00696
00697 for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
00698 std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
00699 fphidata[ih].f() = std::abs(c);
00700 fphidata[ih].phi() = std::arg(c);
00701 }
00702 } else {
00703
00704 FFTmap_p1 fftmap( grid_sampling() );
00705
00706 ffttype f;
00707 int sym;
00708 for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
00709 f = (*this)[ix];
00710 if ( f != 0.0 ) {
00711 fftmap.real_data( ix.coord() ) = f;
00712 for ( sym = 1; sym < cacheref.data().nsym; sym++ )
00713 fftmap.real_data(
00714 ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
00715 }
00716 }
00717
00718 fftmap.fft_x_to_h(cell().volume());
00719
00720 typename H::HKL_reference_index ih;
00721 for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
00722 std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
00723 fphidata[ih].f() = std::abs(c);
00724 fphidata[ih].phi() = std::arg(c);
00725 }
00726 }
00727 }
00728
00729
00732 template<class T> const T& Xmap<T>::operator =( const T& value )
00733 {
00734
00735 for ( Map_reference_index im = first(); !im.last(); im.next() )
00736 list[im.index()] = value;
00737 return value;
00738 }
00739
00740
00742 template<class T> const Xmap<T>& Xmap<T>::operator +=( const Xmap<T>& other )
00743 {
00744 if ( spacegroup().hash() != other.spacegroup().hash() ||
00745 grid_sampling() != other.grid_sampling() )
00746 Message::message( Message_fatal( "Xmap: map mismatch in +=" ) );
00747 for ( Map_reference_index im = first(); !im.last(); im.next() )
00748 list[im.index()] += other[im];
00749 return (*this);
00750 }
00751
00753 template<class T> const Xmap<T>& Xmap<T>::operator -=( const Xmap<T>& other )
00754 {
00755 if ( spacegroup().hash() != other.spacegroup().hash() ||
00756 grid_sampling() != other.grid_sampling() )
00757 Message::message( Message_fatal( "Xmap: map mismatch in -=" ) );
00758 for ( Map_reference_index im = first(); !im.last(); im.next() )
00759 list[im.index()] -= other[im];
00760 return (*this);
00761 }
00762
00763
00764 }
00765
00766 #endif