1 // Clipper structure factor calculation demo
2 /* (C) 2002 Kevin Cowtan */
3 // This is more of a demo application than a serious version
4 
5 #include <clipper/clipper.h>
6 #include <clipper/clipper-ccp4.h>
7 #include <clipper/clipper-contrib.h>
8 #include <clipper/clipper-mmdb.h>
9 
10 #include <iostream>
11 #include <algorithm>
12 
13 
14 using namespace clipper;
15 
sfcalc(HKL_data<datatypes::F_phi<float>> & fphidata,const Atom_list & atoms)16 void sfcalc( HKL_data<datatypes::F_phi<float> >& fphidata, const Atom_list& atoms )
17 {
18   // prepare target map
19   const HKL_info& hkls = fphidata.base_hkl_info();
20   const Grid_sampling grid( hkls.spacegroup(), hkls.cell(), hkls.resolution() );
21   Xmap<float> xmap( hkls.spacegroup(), hkls.cell(), grid );
22 
23   // work out how big a box we need to calc density over for each atom
24   Grid_range gd( hkls.cell(), grid, 3.0 );
25   Xmap<float>::Map_reference_coord i0, iu, iv, iw;
26   // loop over atoms
27   for ( int i = 0; i < atoms.size(); i++ ) if ( !atoms[i].is_null() ) {
28     AtomShapeFn sf( atoms[i] );  // get atom shape fn
29     Coord_frac uvw = atoms[i].coord_orth().coord_frac( hkls.cell() );
30     Coord_grid g0 = uvw.coord_grid( grid ) + gd.min();
31     Coord_grid g1 = uvw.coord_grid( grid ) + gd.max();
32     i0 = Xmap<float>::Map_reference_coord( xmap, g0 );
33     // sum all map contributions from this atoms
34     for ( iu = i0; iu.coord().u() <= g1.u(); iu.next_u() )
35       for ( iv = iu; iv.coord().v() <= g1.v(); iv.next_v() )
36         for ( iw = iv; iw.coord().w() <= g1.w(); iw.next_w() )
37           xmap[iw] += sf.rho( iw.coord_orth() );
38   }
39 
40   // now correct for multiplicity of points on special positions
41   for ( Xmap<float>::Map_reference_index ix = xmap.first();
42 	!ix.last(); ix.next() )
43     xmap[ix] *= xmap.multiplicity( ix.coord() );
44 
45   // calc structure factors from map by fft
46   xmap.fft_to( fphidata );
47 }
48 
49 
main(int argc,char ** argv)50 int main( int argc, char** argv )
51 {
52   // make all the data objects
53   clipper::CSpacegroup spgr( "base spgr" );
54   clipper::CCell cell( spgr, "base cell" );
55   clipper::CResolution reso( cell, "base reso", clipper::Resolution(3.0) );
56   clipper::CHKL_info hkls( reso, "hkls", true );
57   clipper::CHKL_data<clipper::data32::F_phi> fphi1( hkls );
58   clipper::CHKL_data<clipper::data32::F_phi> fphi2( hkls );
59   clipper::CHKL_data<clipper::data32::F_phi> fphi3( hkls );
60   clipper::CHKL_data<clipper::data32::F_phi> fphi4( hkls );
61   clipper::CHKL_data<clipper::data32::F_phi> fphi5( hkls );
62   clipper::CHKL_data<clipper::data32::F_phi> fphi6( hkls );
63   clipper::CGrid_sampling grid( reso );
64   clipper::CXmap<float> xmap( grid );
65 
66   // atomic model
67   clipper::MMDBManager mmdb;
68   mmdb.ReadPDBASCII( argv[1] );
69   //mmdb.write_file( "out.cif", 1 );
70 
71   // get info from mmdb
72   spgr.init( mmdb.spacegroup() );
73   cell.init( mmdb.cell() );
74 
75   // make a second mmdb as a test
76   clipper::MMDBManager mmdb2;
77   mmdb2.set_spacegroup( spgr ); mmdb2.set_cell( cell );
78   mmdb2.cell().debug();         mmdb2.spacegroup().debug();
79 
80   // calc Z's
81   {
82     clipper::AtomSF sf( "C", 0.25 );
83     clipper::HKL_info::HKL_reference_index ih;
84     for ( ih = hkls.first(); !ih.last(); ih.next() )
85       if (ih.index()%100000 == 0) {
86 	std::cout << ih.hkl().format() << " " <<
87 	  sf.f_iso(ih.hkl().coord_reci_orth(cell).invresolsq()) << " " <<
88 	  sf.f_aniso(ih.hkl().coord_reci_orth(cell)) << " " <<
89 	  ih.hkl().coord_reci_orth(cell).invresolsq() << " " <<
90 	  sf.f_iso(ih.hkl().coord_reci_orth(cell).invresolsq()) /
91 	  sf.f_aniso(ih.hkl().coord_reci_orth(cell)) << "\n";
92       }
93   }
94 
95 
96   // debug info
97   spgr.clipper::Container::debug();
98   spgr.clipper::Spacegroup::debug();
99   cell.clipper::Cell::debug();
100   hkls.clipper::HKL_info::debug();
101 
102   // get a list of all the atoms
103   //clipper::DBAtom_selection atoms = mmdb.select_atoms_serial( 0, 0 );
104   clipper::mmdb::PPCAtom psel;
105   int hndl, nsel;
106   hndl = mmdb.NewSelection();
107   mmdb.SelectAtoms( hndl, 0, 0, ::mmdb::SKEY_NEW );
108   mmdb.GetSelIndex( hndl, psel, nsel );
109   clipper::MMDBAtom_list atoms( psel, nsel );
110   mmdb.DeleteSelection( hndl );
111 
112   clipper::SFcalc_iso_sum<float>( fphi1, atoms );
113   clipper::SFcalc_iso_fft<float>( fphi2, atoms, 2.5, 2.0, 0.25 );
114   clipper::SFcalc_aniso_sum<float>( fphi3, atoms );
115   clipper::SFcalc_aniso_fft<float>( fphi4, atoms, 2.5, 2.0, 0.25 );
116 
117   /*
118   sfcalc( fphi5, atoms );
119   for ( HKL_info::HKL_reference_index ih = hkls.first(); !ih.last(); ih.next() ) {
120     std::cout << ih.hkl().format() << " " << fphi4[ih].f() << " " << fphi3[ih].f() << "\n";
121     std::cout << "                        " << clipper::Util::rad2d(fphi4[ih].phi()) << " " << clipper::Util::rad2d(fphi3[ih].phi()) << "\n";
122   }
123   */
124 
125   /*
126   clipper::CCP4MAPfile mapout;
127   mapout.open_write( "out.map" );
128   mapout.export_xmap( xmap );
129   mapout.close_write();
130   clipper::CCP4MTZfile mtzout;
131   mtzout.open_write( "out.mtz" );
132   mtzout.export_hkl_info( hkls );
133   mtzout.export_hkl_data( fphi1, clipper::MTZdataset("dset",1.0), clipper::MTZcrystal("xtal","proj",cell), "[FC1 PHIC1]" );
134   mtzout.export_hkl_data( fphi2, clipper::MTZdataset("dset",1.0), clipper::MTZcrystal("xtal","proj",cell), "[FC2 PHIC2]" );
135   mtzout.export_hkl_data( fphi3, clipper::MTZdataset("dset",1.0), clipper::MTZcrystal("xtal","proj",cell), "[FC3 PHIC3]" );
136   mtzout.export_hkl_data( fphi4, clipper::MTZdataset("dset",1.0), clipper::MTZcrystal("xtal","proj",cell), "[FC4 PHIC4]" );
137   mtzout.close_write();
138   */
139 
140   clipper::HKL_info::HKL_reference_index ih;
141   clipper::HKL_data<clipper::data32::E_sigE> esigdata( hkls );
142   clipper::AtomSF sf( atoms[0].element(), 0.0 );
143   for ( ih = hkls.first(); !ih.last(); ih.next() ) {
144     esigdata[ih].E() = fphi3[ih].f() / sf.f_iso( ih.invresolsq() );
145     esigdata[ih].sigE() = 1.0;
146   }
147 
148   std::vector<clipper::ftype> p( 7, 0.0 );
149   clipper::BasisFn_aniso_gaussian basisfn;
150   clipper::TargetFn_meanEnth<clipper::data32::E_sigE> targetfn( esigdata, 2.0 );
151   clipper::ResolutionFn_nonlinear rfn( hkls, basisfn, targetfn, p );
152   std::cout << basisfn.scale( rfn.params() ) << "\n" << basisfn.u_aniso_orth( rfn.params() ).format() << "\n";
153 
154   // check the answers
155 
156   for ( ih = hkls.first(); !ih.last(); ih.next() ) {
157     std::cout << ih.hkl().format() << " " << fphi1[ih].f() << " " << fphi2[ih].f() << " " << fphi3[ih].f() << " " << fphi4[ih].f() << "\n";
158     std::cout << "                        " << clipper::Util::rad2d(fphi1[ih].phi()) << " " << clipper::Util::rad2d(fphi2[ih].phi()) << " " << clipper::Util::rad2d(fphi3[ih].phi()) << " " << clipper::Util::rad2d(fphi4[ih].phi()) << "\n";
159   }
160 
161   // now test ed calc
162 
163   clipper::Grid_sampling gsam( spgr, cell, reso, 2.0 );
164   clipper::Grid_range gmap( clipper::Coord_grid( 1, 2, 3 ),
165 			  clipper::Coord_grid( 9, 8, 7 ) );
166   clipper::Xmap<float> exmap( spgr, cell, gsam );
167   clipper::NXmap<float> enxmap( cell, gsam, gmap );
168 
169   clipper::EDcalc_iso<float> ediso;
170   clipper::EDcalc_aniso<float> edaniso;
171 
172   ediso( enxmap, atoms );
173 
174   ediso( exmap, atoms );
175   exmap.fft_to( fphi5 );
176   edaniso( exmap, atoms );
177   exmap.fft_to( fphi6 );
178 
179   for ( ih = hkls.first(); !ih.last(); ih.next() )
180     if ( std::abs(std::complex<float>(fphi5[ih]) - std::complex<float>(fphi2[ih]) ) > 5.0 ) std::cout << "err" << ih.hkl().format() << " " << fphi5[ih].f() << " " << fphi2[ih].f() << " " << fphi6[ih].f() << " " << fphi4[ih].f() << "\n";
181 
182   clipper::CCP4MAPfile iomap;
183   iomap.open_write( "map.map" );
184   iomap.export_xmap( exmap );
185   iomap.close_write();
186   clipper::EDcalc_mask<float> edmask;
187   edmask( exmap, atoms );
188   iomap.open_write( "mask.map" );
189   iomap.export_xmap( exmap );
190   iomap.close_write();
191 
192   /*
193   for ( int i = 0; i < atoms.size(); i++ )
194     if ( atoms[i].is_atom() )
195       std::cout << atoms[i].u_aniso_orth().format() << "\n" << atoms[i].u_aniso_orth().sqrt().format() << "\n" << (atoms[i].u_aniso_orth().sqrt()*atoms[i].u_aniso_orth().sqrt()).format() << "\n\n" ;
196   */
197 
198   // now test some selections
199   /*
200   clipper::DBAtom_selection s1 = mmdb.select_atoms( "16-17" );
201   clipper::DBAtom_selection s2 = mmdb.select_atoms( "15-20/CA[C]" );
202   clipper::DBAtom_selection s3 = s1 & s2;
203   clipper::DBAtom_selection s4 = s1 | s2;
204   clipper::DBAtom_selection s5 = s1 ^ s2;
205   for ( int i = 0; i < s3.size(); i++ )
206     std::cout << i << " " << s3[i].residue().seqnum() << "\t" << s3[i].type() << "\n";
207   for ( int i = 0; i < s4.size(); i++ )
208     std::cout << i << " " << s4[i].residue().seqnum() << "\t" << s4[i].type() << "\n";
209   for ( int i = 0; i < s5.size(); i++ )
210     std::cout << i << " " << s5[i].residue().seqnum() << "\t" << s5[i].type() << "\n";
211   */
212 
213   // test AtomSF objects and gradients
214   clipper::Coord_orth co( 1.0, 2.0, 3.0 );
215   clipper::AtomSF sf1( "N", 0.25, 1.0 );
216   clipper::AtomShapeFn sf2( co, "N", 0.25, 1.0 );
217   clipper::AtomShapeFn sfx( co+clipper::Coord_orth(0.001,0.0,0.0), "N", 0.25, 1.0 );
218   clipper::AtomShapeFn sfy( co+clipper::Coord_orth(0.0,0.001,0.0), "N", 0.25, 1.0 );
219   clipper::AtomShapeFn sfz( co+clipper::Coord_orth(0.0,0.0,0.001), "N", 0.25, 1.0 );
220   clipper::AtomShapeFn sfo( co, "N", 0.25, 1.001 );
221   clipper::AtomShapeFn sfu( co, "N", 0.251, 1.0 );
222   clipper::AtomShapeFn sfx2( co+clipper::Coord_orth(0.002,0.0,0.0), "N", 0.25, 1.0 );
223   clipper::AtomShapeFn sfy2( co+clipper::Coord_orth(0.0,0.002,0.0), "N", 0.25, 1.0 );
224   clipper::AtomShapeFn sfz2( co+clipper::Coord_orth(0.0,0.0,0.002), "N", 0.25, 1.0 );
225   clipper::AtomShapeFn sfo2( co, "N", 0.25, 1.002 );
226   clipper::AtomShapeFn sfu2( co, "N", 0.252, 1.0 );
227   clipper::AtomShapeFn sfxy( co+clipper::Coord_orth(0.001,0.001,0.0), "N", 0.25, 1.0 );
228   clipper::AtomShapeFn sfyz( co+clipper::Coord_orth(0.0,0.001,0.001), "N", 0.25, 1.0 );
229   clipper::AtomShapeFn sfzx( co+clipper::Coord_orth(0.001,0.0,0.001), "N", 0.25, 1.0 );
230   std::vector<clipper::AtomShapeFn::TYPE> params;
231   params.push_back( clipper::AtomShapeFn::X );
232   params.push_back( clipper::AtomShapeFn::Y );
233   params.push_back( clipper::AtomShapeFn::Z );
234   params.push_back( clipper::AtomShapeFn::Occ );
235   params.push_back( clipper::AtomShapeFn::Uiso );
236   sf2.agarwal_params() = params;
237   for ( int i = 0; i < 50; i++ ) {
238     clipper::Coord_orth c1 = co + clipper::Coord_orth( 0.1*double(i), 0.0, 0.0 );
239     clipper::Coord_orth c2 = co + clipper::Coord_orth( 0.0, 0.07*double(i), 0.1*double(i) );
240     std::cout << sf1.rho_iso( (c1-co).lengthsq() ) << "\t" << sf2.rho( c1 ) << "\t-\t" << sf1.rho_iso( (c2-co).lengthsq() ) << "\t" << sf2.rho( c2 ) << "\n";
241     std::cout << sf1.f_iso( 0.01*double(i) ) << "\t" << sf2.f( 0.01*double(i) ) << "\n";
242     double f;
243     std::vector<double> g(6);
244     Matrix<double> c(6,6);
245     sf2.rho_curv( c2, f, g, c );
246     std::cout << "G> " << (sfx.rho(c2)-sf2.rho(c2))/0.001 << ":" << g[0] << "\t" << (sfy.rho(c2)-sf2.rho(c2))/0.001 << ":" << g[1] << "\t" << (sfz.rho(c2)-sf2.rho(c2))/0.001 << ":" << g[2] << "\t" << (sfo.rho(c2)-sf2.rho(c2))/0.001 << ":" << g[3] << "\t" << (sfu.rho(c2)-sf2.rho(c2))/0.001 << ":" << g[4] << "\n";
247     std::cout << "C> " << (sfx2.rho(c2)-2*sfx.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(0,0) << "\t" << (sfy2.rho(c2)-2*sfy.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(1,1) << "\t" << (sfz2.rho(c2)-2*sfz.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(2,2) << "\t" << (sfo2.rho(c2)-2*sfo.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(3,3) << "\t" << (sfu2.rho(c2)-2*sfu.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(4,4) << "\n";
248     std::cout << "c> " << (sfxy.rho(c2)-sfx.rho(c2)-sfy.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(0,1) << "\t" << (sfyz.rho(c2)-sfy.rho(c2)-sfz.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(1,2) << "\t" << (sfzx.rho(c2)-sfz.rho(c2)-sfx.rho(c2)+sf2.rho(c2))/0.000001 << ":" << c(2,0) << "\n";
249   }
250 
251   std::cout << "----------------------------------------\n";
252 
253   clipper::U_aniso_orth uani( 0.25, 0.3, 0.35, 0.1, 0.05, 0.0 );
254   clipper::AtomSF sf3( "N", uani, 1.0 );
255   clipper::AtomShapeFn sf4( co, "N", uani, 1.0 );
256   for ( double x = -1.0; x <= 1.0; x+= 0.333333 )
257     for ( double y = -1.0; y <= 1.0; y+= 0.333333 )
258       for ( double z = -1.0; z <= 1.0; z+= 0.333333 ) {
259 	clipper::Coord_orth c(x,y,z);
260 	clipper::Coord_reci_orth r(x,y,z);
261 	std::cout << c.format() << sf3.rho_aniso(c) << "\t" << sf4.rho(co+c) << "\t" << sf3.f_aniso(r) << "\t" << sf4.f(r) << "\n";
262       }
263 
264   std::cout << "----------------------------------------\n";
265 
266   // now test ordinal functions
267   std::vector<clipper::ftype> resols;
268   for ( ih = hkls.first(); !ih.last(); ih.next() )
269     resols.push_back( ih.invresolsq() );
270   clipper::Generic_ordinal ordinal;
271   ordinal.init(resols);
272   clipper::Generic_ordinal ordinv = ordinal;
273   ordinv.invert();
274   clipper::Resolution_ordinal resord;
275   resord.init(hkls,1.0);
276   std::sort( resols.begin(), resols.end() );
277   for ( int i = 0; i < resols.size(); i++ )
278     std::cout << i << " " << resols[i] << " " << ordinal.ordinal(resols[i]) << " " << resord.ordinal(resols[i]) << " " << ordinv.ordinal(ordinal.ordinal(resols[i])) << "\n";
279 
280   std::cout << "----------------------------------------\n";
281 
282   // now test the Ramachandran plot class
283   clipper::Ramachandran rama;
284   std::cout << "\nNonGlyPro\n";
285   rama.init( clipper::Ramachandran::NonGlyPro );
286   for ( int psi = 180; psi >= -180; psi -= 10 ) {
287     std::cout << clipper::String(psi,4) << " " ;
288     for ( int phi = -180; phi <= 180; phi += 10 ) {
289       if      ( rama.favored( clipper::Util::d2rad(phi),
290 			      clipper::Util::d2rad(psi) ) ) std::cout << "#";
291       else if ( rama.allowed( clipper::Util::d2rad(phi),
292 			      clipper::Util::d2rad(psi) ) ) std::cout << "+";
293       else                                                  std::cout << "-";
294     }
295     std::cout << "\n";
296   }
297   std::cout << "\nPro\n";
298   rama.init( clipper::Ramachandran::Pro );
299   for ( int psi = 180; psi >= -180; psi -= 10 ) {
300     std::cout << clipper::String(psi,4) << " " ;
301     for ( int phi = -180; phi <= 180; phi += 10 ) {
302       if      ( rama.favored( clipper::Util::d2rad(phi),
303 			      clipper::Util::d2rad(psi) ) ) std::cout << "#";
304       else if ( rama.allowed( clipper::Util::d2rad(phi),
305 			      clipper::Util::d2rad(psi) ) ) std::cout << "+";
306       else                                                  std::cout << "-";
307     }
308     std::cout << "\n";
309   }
310   std::cout << "\nGly\n";
311   rama.init( clipper::Ramachandran::Gly );
312   for ( int psi = 180; psi >= -180; psi -= 10 ) {
313     std::cout << clipper::String(psi,4) << " " ;
314     for ( int phi = -180; phi <= 180; phi += 10 ) {
315       if      ( rama.favored( clipper::Util::d2rad(phi),
316 			      clipper::Util::d2rad(psi) ) ) std::cout << "#";
317       else if ( rama.allowed( clipper::Util::d2rad(phi),
318 			      clipper::Util::d2rad(psi) ) ) std::cout << "+";
319       else                                                  std::cout << "-";
320     }
321     std::cout << "\n";
322   }
323   std::cout << "\nAll\n";
324   rama.init( clipper::Ramachandran::All );
325   for ( int psi = 180; psi >= -180; psi -= 10 ) {
326     std::cout << clipper::String(psi,4) << " " ;
327     for ( int phi = -180; phi <= 180; phi += 10 ) {
328       if      ( rama.favored( clipper::Util::d2rad(phi),
329 			      clipper::Util::d2rad(psi) ) ) std::cout << "#";
330       else if ( rama.allowed( clipper::Util::d2rad(phi),
331 			      clipper::Util::d2rad(psi) ) ) std::cout << "+";
332       else                                                  std::cout << "-";
333     }
334     std::cout << "\n";
335   }
336 }
337