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_cachekey
00058 {
00059 public:
00060 Xmap_cachekey( const Spgr_descr& spgr_descr, const Grid_sampling& grid ) :
00061 spgr_descr_(spgr_descr), grid_sampling_(grid) {}
00062 const Spgr_descr& spgr_descr() const { return spgr_descr_; }
00063 const Grid_sampling& grid_sampling() const { return grid_sampling_; }
00064 private:
00065 Spgr_descr spgr_descr_;
00066 Grid_sampling grid_sampling_;
00067 };
00068
00069 class Xmap_cacheobj
00070 {
00071 public:
00072 Xmap_cacheobj( const Xmap_cachekey& xmap_cachekey );
00073 bool matches( const Xmap_cachekey& xmap_cachekey ) const;
00074 String format() const;
00075
00076 Xmap_cachekey 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 };
00087
00088
00090
00099 class Xmap_base
00100 {
00101 public:
00102 enum FFTtype { Default, Normal, Sparse };
00103
00105 bool is_null() const;
00106
00108 const Cell& cell() const { return cell_; }
00110 const Spacegroup& spacegroup() const { return spacegroup_; }
00112 const Grid_sampling& grid_sampling() const { return grid_sam_; }
00114 const Grid_range& grid_asu() const { return cacheref.data().asu_grid; }
00116
00117 inline Coord_grid coord_of( const int& index ) const
00118 { return cacheref.data().map_grid.deindex( index ); }
00120
00122 inline int index_of( const Coord_grid& coord ) const {
00123 if ( cacheref.data().asu_grid.in_grid( coord ) ) {
00124 const int i = cacheref.data().map_grid.index( coord );
00125 if ( asu[ i ] == 0 ) return i;
00126 }
00127 return -1;
00128 }
00130 Coord_grid to_map_unit( const Coord_grid& pos ) const
00131 { return pos.unit( grid_sam_ ); }
00132
00134 const RTop<>& operator_orth_grid() const { return rt_orth_grid; }
00136 const RTop<>& operator_grid_orth() const { return rt_grid_orth; }
00138
00140 inline Coord_orth coord_orth( const Coord_map& cm ) const
00141 { return Coord_orth( rt_grid_orth.rot()*cm ); }
00143
00145 inline Coord_map coord_map( const Coord_orth& co ) const
00146 { return Coord_map ( rt_orth_grid.rot()*co ); }
00147
00149 bool in_map( const Coord_grid& pos ) const { return true; }
00151 template<class I> bool in_map( const Coord_map& cm ) const { return true; }
00152
00154 int multiplicity( const Coord_grid& pos ) const;
00155
00157
00161 class Map_reference_base
00162 {
00163 public:
00165 inline const Xmap_base& base_xmap() const { return *map_; }
00167 inline const int& index() const { return index_; }
00169 bool last() const { return ( index_ >= map_->map_grid.size() ); }
00170 protected:
00172 const Xmap_base* map_;
00174 int index_;
00175 };
00176
00178
00189 class Map_reference_index : public Map_reference_base
00190 {
00191 public:
00193 Map_reference_index() {}
00195 explicit Map_reference_index( const Xmap_base& map )
00196 { map_ = ↦ index_=0; next(); }
00198 Map_reference_index( const Xmap_base& map, const Coord_grid& pos ) { map_ = ↦ int sym; map_->find_sym( pos, index_, sym ); }
00200 inline Coord_grid coord() const
00201 { return map_->map_grid.deindex(index_); }
00203 inline const Coord_orth coord_orth() const
00204 { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
00206 inline Map_reference_index& set_coord( const Coord_grid& pos )
00207 { int sym; map_->find_sym( pos, index_, sym ); return *this; }
00209 inline Map_reference_index& next() {
00210 do {
00211 index_++; if ( last() ) break;
00212 } while ( map_->asu[index_] != 0 );
00213 return *this;
00214 }
00216
00217
00218 inline int index_offset(const int& du,const int& dv,const int& dw) const {
00219 int i = index_ + du*map_->du[0] + dv*map_->dv[0] + dw*map_->dw[0];
00220 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] ) ) ); }
00221 return i;
00222 }
00223
00224
00225
00226
00227 };
00228
00230
00241 class Map_reference_coord : public Map_reference_base
00242 {
00243 public:
00245 Map_reference_coord() {}
00247 explicit Map_reference_coord( const Xmap_base& map )
00248 { map_ = ↦ index_ = 0; next(); }
00250 Map_reference_coord( const Xmap_base& map, const Coord_grid& pos ) {
00251 map_ = ↦
00252 pos_ = pos;
00253 map_->find_sym( pos_, index_, sym_ );
00254 }
00256 inline const Coord_grid& coord() const { return pos_; }
00258 inline const Coord_orth coord_orth() const
00259 { return Coord_orth( map_->rt_grid_orth.rot() * coord().coord_map() ); }
00261 inline const int& sym() const { return sym_; }
00263 Map_reference_coord& set_coord( const Coord_grid& pos );
00265
00266 inline Map_reference_coord& next() {
00267 sym_ = 0;
00268 do {
00269 index_++; if ( last() ) break;
00270 } while ( map_->asu[index_] != 0 );
00271 pos_ = map_->map_grid.deindex(index_);
00272 return *this;
00273 }
00274
00275 inline Map_reference_coord& next_u() { pos_.u()++; index_ += map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00276 inline Map_reference_coord& next_v() { pos_.v()++; index_ += map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00277 inline Map_reference_coord& next_w() { pos_.w()++; index_ += map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00278 inline Map_reference_coord& prev_u() { pos_.u()--; index_ -= map_->du[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00279 inline Map_reference_coord& prev_v() { pos_.v()--; index_ -= map_->dv[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00280 inline Map_reference_coord& prev_w() { pos_.w()--; index_ -= map_->dw[sym_]; if (map_->asu[index_] != 0) edge(); return *this; }
00281
00282 inline Map_reference_coord& operator =( const Coord_grid& pos )
00283 { return set_coord( pos ); }
00284
00285
00286
00287
00288
00289 protected:
00291 int sym_;
00293 Coord_grid pos_;
00294
00296 void edge();
00297 };
00298
00300 Map_reference_index first() const { return Map_reference_index( *this ); }
00302 Map_reference_coord first_coord() const { return Map_reference_coord( *this ); }
00304 static FFTtype& default_type() { return default_type_; }
00305 protected:
00306 ObjectCache<Xmap_cacheobj>::Reference cacheref;
00307 const unsigned char* asu;
00308 const Isymop* isymop;
00309 const int* du;
00310 const int* dv;
00311 const int* dw;
00312 Grid_range asu_grid;
00313 Grid_range map_grid;
00314 int nsym;
00315
00316 Cell cell_;
00317 Spacegroup spacegroup_;
00318 Grid_sampling grid_sam_;
00319
00320 RTop<> rt_orth_grid;
00321 RTop<> rt_grid_orth;
00322
00324 Xmap_base();
00326 void init( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam );
00327 inline void find_sym( const Coord_grid& base, int& index, int& sym ) const;
00328 void asu_error() const;
00329
00330 static FFTtype default_type_;
00331
00332 friend class Xmap_base::Map_reference_base;
00333 friend class Xmap_base::Map_reference_index;
00334 friend class Xmap_base::Map_reference_coord;
00335 };
00336
00337
00338
00339
00341
00355 template<class T> class Xmap : public Xmap_base
00356 {
00357 public:
00359 Xmap() {}
00361 Xmap( const Spacegroup& spacegroup, const Cell& cell, const Grid_sampling& grid_sam ) { init( spacegroup, cell, grid_sam ); }
00363 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() ); }
00364
00366 inline const T& operator[] (const Xmap_base::Map_reference_index& ix) const
00367 { return list[ix.index()]; }
00369 inline T& operator[] (const Xmap_base::Map_reference_index& ix)
00370 { return list[ix.index()]; }
00371
00373 inline const T& operator[] (const Xmap_base::Map_reference_coord& ix) const
00374 { return list[ix.index()]; }
00376 inline T& operator[] (const Xmap_base::Map_reference_coord& ix)
00377 { return list[ix.index()]; }
00378
00380 const T& get_data( const Coord_grid& pos ) const;
00382 void set_data( const Coord_grid& pos, const T& val );
00384 inline const T& get_data( const int& index ) const;
00386 bool set_data( const int& index, const T& val );
00387
00389 template<class I> T interp( const Coord_frac& pos ) const;
00391 template<class I> void interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const;
00393 template<class I> void interp_curv( const Coord_frac& pos, T& val, Grad_frac<T>& grad, Curv_frac<T>& curv ) const;
00395 template<class I> T interp( const Coord_map& pos ) const;
00397 template<class I> void interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const;
00399 template<class I> void interp_curv( const Coord_map& pos, T& val, Grad_map<T>& grad, Curv_map<T>& curv ) const;
00400
00402 template<class H> void fft_from( const H& fphidata, const FFTtype type = Default );
00404 template<class H> void fft_to ( H& fphidata, const FFTtype type = Default ) const;
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00419 const T& operator =( const T& value );
00421 const Xmap<T>& operator +=( const Xmap<T>& other );
00423 const Xmap<T>& operator -=( const Xmap<T>& other );
00424
00425 private:
00426 std::vector<T> list;
00427 };
00428
00429
00430
00431
00432
00433 void Xmap_base::find_sym( const Coord_grid& base, int& index, int& sym ) const
00434 {
00435
00436 Coord_grid rot = to_map_unit( base );
00437 if ( asu_grid.in_grid( rot ) ) {
00438 index = map_grid.index( rot );
00439 if ( asu[ index ] == 0 ) {
00440 sym = 0;
00441 } else {
00442 sym = asu[ index ] - 1;
00443 index = map_grid.index( to_map_unit( base.transform(isymop[sym]) ) );
00444 }
00445 } else {
00446
00447 for ( sym = 1; sym < nsym; sym++ ) {
00448 rot = to_map_unit( base.transform( isymop[sym] ) );
00449 if ( asu_grid.in_grid( rot ) ) {
00450 index = map_grid.index( rot );
00451 if ( asu[ index ] == 0 ) return;
00452 }
00453 }
00454 index = 0;
00455 asu_error();
00456 }
00457 return;
00458 }
00459
00460
00465 template<class T> const T& Xmap<T>::get_data( const Coord_grid& pos ) const
00466 {
00467 int index, sym;
00468 find_sym( pos, index, sym );
00469 return list[ index ];
00470 }
00471
00476 template<class T> void Xmap<T>::set_data( const Coord_grid& pos, const T& val )
00477 {
00478 int index, sym;
00479 find_sym( pos, index, sym );
00480 list[ index ] = val;
00481 }
00482
00489 template<class T> const T& Xmap<T>::get_data( const int& index ) const
00490 { return list[index]; }
00491
00499 template<class T> bool Xmap<T>::set_data( const int& index, const T& val )
00500 {
00501 if ( index >= 0 && index < list.size() )
00502 if ( asu[index] == 0 ) {
00503 list[index] = val;
00504 return true;
00505 }
00506 return false;
00507 }
00508
00517 template<class T> template<class I> T Xmap<T>::interp( const Coord_frac& pos ) const
00518 {
00519 T val;
00520 I::interp( *this, pos.coord_map( grid_sam_ ), val );
00521 return val;
00522 }
00523
00524
00532 template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_frac& pos, T& val, Grad_frac<T>& grad ) const
00533 {
00534 Grad_map<T> g;
00535 I::interp_grad( *this, pos.coord_map( grid_sam_ ), val, g );
00536 grad = g.grad_frac( grid_sam_ );
00537 }
00538
00539
00549 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
00550 {
00551 Grad_map<T> g;
00552 Curv_map<T> c;
00553 I::interp_curv( *this, pos.coord_map( grid_sam_ ), val, g, c );
00554 grad = g.grad_frac( grid_sam_ );
00555 curv = c.curv_frac( grid_sam_ );
00556 }
00557
00558
00567 template<class T> template<class I> T Xmap<T>::interp( const Coord_map& pos ) const
00568 {
00569 T val;
00570 I::interp( *this, pos, val );
00571 return val;
00572 }
00573
00574
00582 template<class T> template<class I> void Xmap<T>::interp_grad( const Coord_map& pos, T& val, Grad_map<T>& grad ) const
00583 { I::interp_grad( *this, pos, val, grad ); }
00584
00585
00595 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
00596 { I::interp_curv( *this, pos, val, grad, curv ); }
00597
00598
00603 template<class T> template<class H> void Xmap<T>::fft_from( const H& fphidata, const FFTtype type )
00604 {
00605 if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
00606
00607 FFTmap_sparse_p1_hx fftmap( grid_sampling() );
00608
00609 typename H::HKL_reference_index ih;
00610 ffttype f, phi0, phi1;
00611 int sym;
00612 for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
00613 f = fphidata[ih].f();
00614 if ( f != 0.0 ) {
00615 phi0 = fphidata[ih].phi();
00616 const HKL& hkl = ih.hkl();
00617 fftmap.set_hkl( hkl,
00618 std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
00619 for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
00620 phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
00621 fftmap.set_hkl( hkl.transform( isymop[sym] ),
00622 std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
00623 }
00624 }
00625 }
00626
00627 for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
00628 fftmap.require_real_data( ix.coord() );
00629
00630 fftmap.fft_h_to_x(1.0/cell().volume());
00631
00632 for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
00633 (*this)[ix] = fftmap.real_data( ix.coord() );
00634 } else {
00635
00636 FFTmap_p1 fftmap( grid_sampling() );
00637
00638 typename H::HKL_reference_index ih;
00639 ffttype f, phi0, phi1;
00640 int sym;
00641 for ( ih = fphidata.first_data(); !ih.last(); fphidata.next_data( ih ) ) {
00642 f = fphidata[ih].f();
00643 if ( f != 0.0 ) {
00644 phi0 = fphidata[ih].phi();
00645 const HKL& hkl = ih.hkl();
00646 fftmap.set_hkl( hkl,
00647 std::complex<ffttype>( f*cos(phi0), f*sin(phi0) ) );
00648 for ( sym = 1; sym < spacegroup_.num_primops(); sym++ ) {
00649 phi1 = phi0 + hkl.sym_phase_shift( spacegroup_.symop(sym) );
00650 fftmap.set_hkl( hkl.transform( isymop[sym] ),
00651 std::complex<ffttype>( f*cos(phi1), f*sin(phi1) ) );
00652 }
00653 }
00654 }
00655
00656 fftmap.fft_h_to_x(1.0/cell().volume());
00657
00658 for ( Map_reference_index ix = first(); !ix.last(); ix.next() )
00659 (*this)[ix] = fftmap.real_data( ix.coord() );
00660 }
00661 }
00662
00663
00672 template<class T> template<class H> void Xmap<T>::fft_to ( H& fphidata, const FFTtype type ) const
00673 {
00674 if ( type == Sparse || ( type == Default && default_type() == Sparse ) ) {
00675
00676 FFTmap_sparse_p1_xh fftmap( grid_sampling() );
00677
00678 ffttype f;
00679 int sym;
00680 for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
00681 f = (*this)[ix];
00682 if ( f != 0.0 ) {
00683 fftmap.real_data( ix.coord() ) = f;
00684 for ( sym = 1; sym < cacheref.data().nsym; sym++ )
00685 fftmap.real_data(
00686 ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
00687 }
00688 }
00689
00690 typename H::HKL_reference_index ih;
00691 for ( ih = fphidata.first(); !ih.last(); ih.next() )
00692 fftmap.require_hkl( ih.hkl() );
00693
00694 fftmap.fft_x_to_h(cell().volume());
00695
00696 for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
00697 std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
00698 fphidata[ih].f() = std::abs(c);
00699 fphidata[ih].phi() = std::arg(c);
00700 }
00701 } else {
00702
00703 FFTmap_p1 fftmap( grid_sampling() );
00704
00705 ffttype f;
00706 int sym;
00707 for ( Map_reference_index ix = first(); !ix.last(); ix.next() ) {
00708 f = (*this)[ix];
00709 if ( f != 0.0 ) {
00710 fftmap.real_data( ix.coord() ) = f;
00711 for ( sym = 1; sym < cacheref.data().nsym; sym++ )
00712 fftmap.real_data(
00713 ix.coord().transform( isymop[sym] ).unit( grid_sam_ ) ) = f;
00714 }
00715 }
00716
00717 fftmap.fft_x_to_h(cell().volume());
00718
00719 typename H::HKL_reference_index ih;
00720 for ( ih = fphidata.first(); !ih.last(); ih.next() ) {
00721 std::complex<ffttype> c = fftmap.get_hkl( ih.hkl() );
00722 fphidata[ih].f() = std::abs(c);
00723 fphidata[ih].phi() = std::arg(c);
00724 }
00725 }
00726 }
00727
00728
00731 template<class T> const T& Xmap<T>::operator =( const T& value )
00732 {
00733
00734 for ( Map_reference_index im = first(); !im.last(); im.next() )
00735 list[im.index()] = value;
00736 return value;
00737 }
00738
00739
00741 template<class T> const Xmap<T>& Xmap<T>::operator +=( const Xmap<T>& other )
00742 {
00743 if ( spacegroup().hash() != other.spacegroup().hash() ||
00744 grid_sampling() != other.grid_sampling() )
00745 Message::message( Message_fatal( "Xmap: map mismatch in +=" ) );
00746 for ( Map_reference_index im = first(); !im.last(); im.next() )
00747 list[im.index()] += other[im];
00748 return (*this);
00749 }
00750
00752 template<class T> const Xmap<T>& Xmap<T>::operator -=( const Xmap<T>& other )
00753 {
00754 if ( spacegroup().hash() != other.spacegroup().hash() ||
00755 grid_sampling() != other.grid_sampling() )
00756 Message::message( Message_fatal( "Xmap: map mismatch in -=" ) );
00757 for ( Map_reference_index im = first(); !im.last(); im.next() )
00758 list[im.index()] -= other[im];
00759 return (*this);
00760 }
00761
00762
00763 }
00764
00765 #endif