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