00001
00002 #ifndef CCTBX_MAPTBX_AVERAGE_DENSITIES_H
00003 #define CCTBX_MAPTBX_AVERAGE_DENSITIES_H
00004
00005 #include <cctbx/uctbx.h>
00006 #include <scitbx/array_family/accessors/c_grid.h>
00007 #include <scitbx/math/utils.h>
00008
00009 namespace cctbx { namespace maptbx {
00010
00012 inline
00013 int
00014 nint(double x)
00015 {
00016 return int(std::ceil(x+0.5)-(std::fmod(x*0.5+0.25,1.0)!=0));
00017 }
00018
00019 af::versa<double, af::c_grid<3> > box_map_averaging(
00020 cctbx::uctbx::unit_cell const& unit_cell,
00021 af::const_ref<double, af::c_grid<3> > const& data,
00022 double const& rad)
00023 {
00024 int lx,ly,lz,kx,ky,kz,mx,my,mz,x1box,x2box,y1box,y2box,z1box,z2box;
00025 int NX = data.accessor()[0];
00026 int NY = data.accessor()[1];
00027 int NZ = data.accessor()[2];
00028 double xrad = rad*unit_cell.reciprocal_parameters()[0]*NX;
00029 double yrad = rad*unit_cell.reciprocal_parameters()[1]*NY;
00030 double zrad = rad*unit_cell.reciprocal_parameters()[2]*NZ;
00031 af::versa<double, af::c_grid<3> > new_data;
00032 new_data.resize(af::c_grid<3>(NX,NY,NZ), 0.0);
00033 af::ref<double, af::c_grid<3> > new_data_ref = new_data.ref();
00034
00035 for (lx = 0; lx < NX; lx++) {
00036 for (ly = 0; ly < NY; ly++) {
00037 for (lz = 0; lz < NZ; lz++) {
00038 double r_ave_xyz = 0.0;
00039 int counter = 0;
00040 x1box=nint(static_cast<double>(lx)-xrad) - 1;
00041 x2box=nint(static_cast<double>(lx)+xrad) + 1;
00042 y1box=nint(static_cast<double>(ly)-yrad) - 1;
00043 y2box=nint(static_cast<double>(ly)+yrad) + 1;
00044 z1box=nint(static_cast<double>(lz)-zrad) - 1;
00045 z2box=nint(static_cast<double>(lz)+zrad) + 1;
00046 for (kx = x1box; kx <= x2box; kx++) {
00047 for (ky = y1box; ky <= y2box; ky++) {
00048 for (kz = z1box; kz <= z2box; kz++) {
00049 mx = scitbx::math::mod_positive(kx, NX);
00050 my = scitbx::math::mod_positive(ky, NY);
00051 mz = scitbx::math::mod_positive(kz, NZ);
00052 r_ave_xyz += data(mx,my,mz);
00053 counter += 1;
00054 }}}
00055 new_data_ref(lx,ly,lz) = r_ave_xyz / counter;
00056 }}}
00057 return new_data;
00058 }
00059
00060
00061 class grid_points_in_sphere_around_atom_and_distances {
00062 public:
00063 grid_points_in_sphere_around_atom_and_distances(
00064 cctbx::uctbx::unit_cell const& uc,
00065 af::const_ref<double, af::c_grid<3> > const& data,
00066 double const& rad,
00067 double const& shell,
00068 scitbx::vec3<double> const& site_frac)
00069 {
00070 int nx = data.accessor()[0];
00071 int ny = data.accessor()[1];
00072 int nz = data.accessor()[2];
00073 double distsq;
00074 double grid_step = uc.parameters()[3] / nx;
00075 int mx,my,mz,kx,ky,kz;
00076 double abc = uc.parameters()[0]*uc.parameters()[1]*uc.parameters()[2];
00077 double uc_volume_over_abc = uc.volume() / abc;
00078 double sin_alpha = std::sin(scitbx::deg_as_rad(uc.parameters()[3]));
00079 double sin_beta = std::sin(scitbx::deg_as_rad(uc.parameters()[4]));
00080 double sin_gamma = std::sin(scitbx::deg_as_rad(uc.parameters()[5]));
00081 double ascale = uc_volume_over_abc / sin_alpha;
00082 double bscale = uc_volume_over_abc / sin_beta;
00083 double cscale = uc_volume_over_abc / sin_gamma;
00084 double as = uc.parameters()[0]/ascale;
00085 double bs = uc.parameters()[1]/bscale;
00086 double cs = uc.parameters()[2]/cscale;
00087 double xshell = shell/uc.parameters()[0]/ascale;
00088 double yshell = shell/uc.parameters()[1]/bscale;
00089 double zshell = shell/uc.parameters()[2]/cscale;
00090 double mr1= uc.metrical_matrix()[0];
00091 double mr2= uc.metrical_matrix()[3];
00092 double mr3= uc.metrical_matrix()[4];
00093 double mr5= uc.metrical_matrix()[1];
00094 double mr6= uc.metrical_matrix()[5];
00095 double mr9= uc.metrical_matrix()[2];
00096 double tmr2 = mr2*2.0;
00097 double tmr3 = mr3*2.0;
00098 double tmr6 = mr6*2.0;
00099 double xL = -xshell;
00100 double xR = 1.0+xshell;
00101 double yL = -yshell;
00102 double yR = 1.0+yshell;
00103 double zL = -zshell;
00104 double zR = 1.0+zshell;
00105 double xf = site_frac[0];
00106 double yf = site_frac[1];
00107 double zf = site_frac[2];
00108 double radsq = rad * rad;
00109 bool close_to_au=(xf>=xL||xf<=xR)&&(yf>=yL||yf<=yR)&&(zf>=zL||zf<=zR);
00110 if(close_to_au) {
00111 double coas = rad/as;
00112 double cobs = rad/bs;
00113 double cocs = rad/cs;
00114 int x1box = scitbx::math::nearest_integer( nx*(xf-coas) ) - 1;
00115 int x2box = scitbx::math::nearest_integer( nx*(xf+coas) ) + 1;
00116 int y1box = scitbx::math::nearest_integer( ny*(yf-cobs) ) - 1;
00117 int y2box = scitbx::math::nearest_integer( ny*(yf+cobs) ) + 1;
00118 int z1box = scitbx::math::nearest_integer( nz*(zf-cocs) ) - 1;
00119 int z2box = scitbx::math::nearest_integer( nz*(zf+cocs) ) + 1;
00120 for(kx = x1box; kx <= x2box; kx++) {
00121 double xn=xf-double(kx)/nx;
00122 for(ky = y1box; ky <= y2box; ky++) {
00123 double yn=yf-double(ky)/ny;
00124 for(kz = z1box; kz <= z2box; kz++) {
00125 double zn=zf-double(kz)/nz;
00126 distsq=mr1*xn*xn+mr5*yn*yn+mr9*zn*zn+
00127 tmr2*xn*yn+tmr3*xn*zn+tmr6*yn*zn;
00128 if(distsq <= radsq) {
00129 mx = scitbx::math::mod_positive(kx, nx);
00130 my = scitbx::math::mod_positive(ky, ny);
00131 mz = scitbx::math::mod_positive(kz, nz);
00132
00133 data_at_grid_points_.push_back(data(mx,my,mz));
00134
00135 distances_.push_back(std::sqrt(distsq));
00136 }
00137 }
00138 }
00139 }
00140 }
00141
00142 double tolerance = grid_step/25. ;
00143 for(std::size_t i = 0; i < data_at_grid_points_.size(); i++) {
00144 double dist = distances_[i];
00145 double data_ave = data_at_grid_points_[i];
00146 int counter = 1;
00147 for(std::size_t j = 0; j < data_at_grid_points_.size(); j++) {
00148 if(distances_[j]<dist+tolerance && distances_[j]>dist-tolerance && i!=j && std::abs(distances_[i]-distances_[j]) > 1.e-6) {
00149 counter++;
00150 data_ave = data_ave + data_at_grid_points_[j];
00151
00152
00153 }
00154 }
00155 data_at_grid_points_averaged_.push_back(data_ave / counter);
00156 }
00157
00158 }
00159 af::shared<double> data_at_grid_points() { return data_at_grid_points_; }
00160 af::shared<double> data_at_grid_points_averaged() { return data_at_grid_points_averaged_; }
00161 af::shared<double> distances() { return distances_; }
00162 protected:
00163 af::shared<double> data_at_grid_points_, data_at_grid_points_averaged_;
00164 af::shared<double> distances_;
00165 };
00166
00167 class find_gaussian_parameters {
00168 public:
00169 find_gaussian_parameters() {}
00170 find_gaussian_parameters(af::const_ref<double> const& data_at_grid_points,
00171 af::const_ref<double> const& distances,
00172 double const& cutoff_radius,
00173 double const& allowed_region_radius,
00174 double weight_power)
00175 {
00176 CCTBX_ASSERT(cutoff_radius <= allowed_region_radius);
00177 double max_distances = af::max(distances);
00178 CCTBX_ASSERT(cutoff_radius <= max_distances &&
00179 allowed_region_radius <= max_distances);
00180 double p=0.0, q=0.0, r=0.0, s=0.0, t=0.0;
00181 int n = 0;
00182 int number_of_points = data_at_grid_points.size();
00183 for(int i = 0; i < number_of_points; i++) {
00184 if(data_at_grid_points[i] > 0.0 && distances[i] <= cutoff_radius) {
00185 n = n + 1;
00186 double data_i = data_at_grid_points[i];
00187 double distance_i = distances[i];
00188 double distance_i_sq = distance_i*distance_i;
00189 double pow_distance_i = std::pow(distance_i, weight_power);
00190 double log_of_data_i = std::log(data_i);
00191 if(pow_distance_i < 1.e-9) pow_distance_i = 1.0;
00192 p=p+log_of_data_i/pow_distance_i;
00193 q=q+distance_i_sq/pow_distance_i;
00194 r=r+(distance_i_sq*distance_i_sq)/pow_distance_i;
00195 s=s+distance_i_sq*log_of_data_i/pow_distance_i;
00196 t=t+1.0/pow_distance_i;
00197 }
00198 }
00199 CCTBX_ASSERT(r != 0.0);
00200 double den = t-q*q/r;
00201 CCTBX_ASSERT(den != 0.0);
00202 double alpha_opt = (p-s*q/r) / den;
00203 a_real_space_ = std::exp( alpha_opt );
00204 b_real_space_ = 1./r*(alpha_opt*q - s);
00205 double tmp = b_real_space_/scitbx::constants::pi;
00206 double pi_sq = scitbx::constants::pi*scitbx::constants::pi;
00207 CCTBX_ASSERT(tmp != 0.0);
00208 a_reciprocal_space_ = a_real_space_/std::sqrt(tmp*tmp*tmp);
00209 CCTBX_ASSERT(b_real_space_ != 0.0);
00210 b_reciprocal_space_ = pi_sq/b_real_space_*4;
00211 double num = 0.0;
00212 double denum = 0.0;
00213 for(int i = 0; i < number_of_points; i++) {
00214 if(data_at_grid_points[i]>0.0 && distances[i]<=allowed_region_radius) {
00215 double data_i = data_at_grid_points[i];
00216 double distance_i = distances[i];
00217 double distance_i_sq = distance_i*distance_i;
00218 double data_i_approx =
00219 a_real_space_ * std::exp(-b_real_space_*distance_i_sq);
00220 num=num+std::abs(data_i-data_i_approx);
00221 denum=denum+data_i;
00222 }
00223 }
00224 CCTBX_ASSERT(denum != 0.0);
00225 gof_ = num/denum*100.;
00226 }
00227 double a_real_space() { return a_real_space_; }
00228 double b_real_space() { return b_real_space_; }
00229 double a_reciprocal_space() { return a_reciprocal_space_; }
00230 double b_reciprocal_space() { return b_reciprocal_space_; }
00231 double gof() { return gof_; }
00232 protected:
00233 double a_real_space_,b_real_space_,a_reciprocal_space_,b_reciprocal_space_;
00234 double gof_;
00235 };
00236
00237 class one_gaussian_peak_approximation {
00238 public:
00239 one_gaussian_peak_approximation(
00240 af::const_ref<double> const& data_at_grid_points,
00241 af::const_ref<double> const& distances,
00242 bool const& use_weights,
00243 bool const& optimize_cutoff_radius)
00244 {
00245 first_zero_radius_ =find_first_zero_radius(data_at_grid_points, distances);
00246 double radius_increment = 0.01, weight_power = 0.0, weight_increment =0.05;
00247 weight_power_ = -1;
00248 cutoff_radius_ = -1;
00249 gof_ = 999.;
00250 if(use_weights && optimize_cutoff_radius) {
00251 for(double c_radius = 0.1; c_radius <= first_zero_radius_;
00252 c_radius = c_radius+radius_increment) {
00253 for(double w_power=0.0;w_power<=10.0;w_power=w_power+weight_increment) {
00254 find_gaussian_parameters fgp_obj(data_at_grid_points,
00255 distances,
00256 c_radius,
00257 first_zero_radius_,
00258 w_power);
00259 if(fgp_obj.gof() < gof_) {
00260 gof_ = fgp_obj.gof();
00261 weight_power_ = w_power;
00262 cutoff_radius_ = c_radius;
00263 fgp_obj_ = fgp_obj;
00264 }
00265 }}
00266
00267 }
00268 else if(use_weights && !optimize_cutoff_radius) {
00269 for(double w_power=0.0;w_power<=20.0;w_power=w_power+weight_increment) {
00270 find_gaussian_parameters fgp_obj(data_at_grid_points,
00271 distances,
00272 first_zero_radius_,
00273 first_zero_radius_,
00274 w_power);
00275 if(fgp_obj.gof() < gof_) {
00276 gof_ = fgp_obj.gof();
00277 weight_power_ = w_power;
00278 cutoff_radius_ = first_zero_radius_;
00279 fgp_obj_ = fgp_obj;
00280 }
00281 }
00282 }
00283 else if(optimize_cutoff_radius && !use_weights) {
00284 for(double c_radius = 0.1; c_radius <= first_zero_radius_;
00285 c_radius = c_radius+radius_increment) {
00286 find_gaussian_parameters fgp_obj(data_at_grid_points,
00287 distances,
00288 c_radius,
00289 first_zero_radius_,
00290 0.0);
00291 if(fgp_obj.gof() < gof_) {
00292 gof_ = fgp_obj.gof();
00293 cutoff_radius_ = c_radius;
00294 fgp_obj_ = fgp_obj;
00295 }
00296 }
00297 }
00298 else {
00299 find_gaussian_parameters fgp_obj(data_at_grid_points,
00300 distances,
00301 first_zero_radius_,
00302 first_zero_radius_,
00303 0.0);
00304 gof_ = fgp_obj.gof();
00305 cutoff_radius_ = first_zero_radius_;
00306 fgp_obj_ = fgp_obj;
00307 }
00308 }
00309 double find_first_zero_radius(
00310 af::const_ref<double> const& data_at_grid_points,
00311 af::const_ref<double> const& distances)
00312 {
00313 double shift = 1.e-3;
00314 af::shared<double> distances_for_negative_data_at_grid_points;
00315 for(int i = 0; i < data_at_grid_points.size(); i++) {
00316 if(data_at_grid_points[i] < 0.0) {
00317 distances_for_negative_data_at_grid_points.push_back(distances[i]);
00318 }
00319 }
00320 double result = af::max(distances);
00321 if(distances_for_negative_data_at_grid_points.size() > 0) {
00322 result =
00323 af::min(distances_for_negative_data_at_grid_points.const_ref())-shift;
00324 }
00325 return result;
00326 }
00327 double a_real_space() { return fgp_obj_.a_real_space(); }
00328 double b_real_space() { return fgp_obj_.b_real_space(); }
00329 double a_reciprocal_space() { return fgp_obj_.a_reciprocal_space(); }
00330 double b_reciprocal_space() { return fgp_obj_.b_reciprocal_space(); }
00331 double gof()
00332 {
00333 CCTBX_ASSERT(gof_ == fgp_obj_.gof());
00334 return gof_;
00335 }
00336 double cutoff_radius() { return cutoff_radius_; }
00337 double weight_power() { return weight_power_; }
00338 double first_zero_radius() { return first_zero_radius_; }
00339 protected:
00340 double a_real_space_,b_real_space_,a_reciprocal_space_,b_reciprocal_space_;
00341 double gof_, cutoff_radius_, weight_power_, first_zero_radius_;
00342 find_gaussian_parameters fgp_obj_;
00343 };
00344
00345 template <typename DataType>
00346 af::shared<DataType>
00347 average_densities(
00348 uctbx::unit_cell const& unit_cell,
00349 af::const_ref<DataType, af::c_grid<3> > const& data,
00350 af::const_ref<scitbx::vec3<double> > const& sites_frac,
00351 float radius)
00352 {
00353
00354 af::shared<DataType> result((af::reserve(sites_frac.size())));
00355 typedef float f_t;
00356 int nx = static_cast<int>(data.accessor()[0]);
00357 int ny = static_cast<int>(data.accessor()[1]);
00358 int nz = static_cast<int>(data.accessor()[2]);
00359 f_t mr1= static_cast<f_t>(unit_cell.metrical_matrix()[0]);
00360 f_t mr5= static_cast<f_t>(unit_cell.metrical_matrix()[1]);
00361 f_t mr9= static_cast<f_t>(unit_cell.metrical_matrix()[2]);
00362
00363 f_t mr2= static_cast<f_t>(unit_cell.metrical_matrix()[3]);
00364
00365 f_t mr3= static_cast<f_t>(unit_cell.metrical_matrix()[4]);
00366
00367 f_t mr6= static_cast<f_t>(unit_cell.metrical_matrix()[5]);
00368 f_t tmr2 = mr2*2;
00369 f_t tmr3 = mr3*2;
00370 f_t tmr6 = mr6*2;
00371 f_t sx = 1/static_cast<f_t>(nx); f_t tsx= sx*2; f_t sxsq=mr1*sx*sx;
00372 f_t sy = 1/static_cast<f_t>(ny); f_t tsy= sy*2; f_t sysq=mr5*sy*sy;
00373 f_t sz = 1/static_cast<f_t>(nz); f_t tsz= sz*2; f_t szsq=mr9*sz*sz;
00374 f_t w1=mr1*sx*tsx; f_t w4=mr5*sy*tsy;
00375 f_t w2=mr2*sx*tsy; f_t w5=mr6*sy*tsz;
00376 f_t w3=mr3*sx*tsz; f_t w6=mr9*sz*tsz;
00377 f_t tsxg1=tsx*mr1; f_t tsyg4=tsy*mr2; f_t tszg3=tsz*mr3;
00378 f_t tsxg4=tsx*mr2; f_t tsyg5=tsy*mr5; f_t tszg8=tsz*mr6;
00379 f_t tsxg7=tsx*mr3; f_t tsyg8=tsy*mr6; f_t tszg9=tsz*mr9;
00380 f_t rp[3];
00381 for(unsigned i=0;i<3;i++) {
00382 rp[i] = static_cast<f_t>(unit_cell.reciprocal_parameters()[i]);
00383 }
00384 std::vector<int> mys;
00385 std::vector<int> mzs;
00386 f_t cutoff = radius;
00387 typedef scitbx::math::float_int_conversions<f_t, int> fic;
00388 for(std::size_t i_site=0;i_site<sites_frac.size();i_site++) {
00389 fractional<> const& site = sites_frac[i_site];
00390 f_t xfi=static_cast<f_t>(site[0]);
00391 f_t yfi=static_cast<f_t>(site[1]);
00392 f_t zfi=static_cast<f_t>(site[2]);
00393 f_t cutoffsq=cutoff*cutoff;
00394 f_t coas = cutoff*rp[0];
00395 int x1box=fic::ifloor(nx*(xfi-coas));
00396 int x2box=fic::iceil(nx*(xfi+coas));
00397 f_t cobs = cutoff*rp[1];
00398 int y1box=fic::ifloor(ny*(yfi-cobs));
00399 int y2box=fic::iceil(ny*(yfi+cobs));
00400 f_t cocs = cutoff*rp[2];
00401 int z1box=fic::ifloor(nz*(zfi-cocs));
00402 int z2box=fic::iceil(nz*(zfi+cocs));
00403 f_t sxbcen=xfi-x1box*sx;
00404 f_t sybcen=yfi-y1box*sy;
00405 f_t szbcen=zfi-z1box*sz;
00406 f_t distsm=mr1*sxbcen*sxbcen+mr5*sybcen*sybcen+mr9*szbcen*szbcen
00407 +tmr2*sxbcen*sybcen+tmr3*sxbcen*szbcen+tmr6*sybcen*szbcen;
00408 f_t w7=tsxg1*sxbcen+tsxg4*sybcen+tsxg7*szbcen;
00409 f_t w8=tsyg4*sxbcen+tsyg5*sybcen+tsyg8*szbcen;
00410 f_t w9=tszg3*sxbcen+tszg8*sybcen+tszg9*szbcen;
00411 mys.clear();
00412 mys.reserve(y2box-y1box+1);
00413 for (int ky = y1box; ky <= y2box; ky++) {
00414 int my = ky % ny;
00415 if (my < 0) my += ny;
00416 mys.push_back(my);
00417 }
00418 mzs.clear();
00419 mzs.reserve(z2box-z1box+1);
00420 for (int kz = z1box; kz <= z2box; kz++) {
00421 int mz = kz % nz;
00422 if (mz < 0) mz += nz;
00423 mzs.push_back(mz);
00424 }
00425 DataType density_sum = 0;
00426 std::size_t n_density_contributions = 0;
00427 f_t distsx = distsm;
00428 f_t s1xx = sxsq - w7;
00429 f_t s1xy = sysq - w8;
00430 f_t s1xz = szsq - w9;
00431 for (int kx = x1box; kx <= x2box; kx++) {
00432 int mx = kx % nx;
00433 if (mx < 0) mx += nx;
00434 int mxny = mx * ny;
00435 f_t s2yz = s1xz;
00436 f_t s2_incr = s1xy;
00437 f_t s2 = distsx;
00438 std::vector<int>::const_iterator mye = mys.end();
00439 for (std::vector<int>::const_iterator
00440 myi=mys.begin();
00441 myi!=mye;
00442 myi++) {
00443 f_t s3_incr = s2yz;
00444 f_t dist = s2;
00445 const DataType* data_mxnymynz = &data[(mxny + (*myi)) * nz];
00446 std::vector<int>::const_iterator mze = mzs.end();
00447 for (std::vector<int>::const_iterator
00448 mzi=mzs.begin();
00449 mzi!=mze;
00450 mzi++) {
00451 if (dist < cutoffsq) {
00452 density_sum += data_mxnymynz[*mzi];
00453 n_density_contributions++;
00454 }
00455 dist += s3_incr;
00456 s3_incr += w6;
00457 }
00458 s2 += s2_incr;
00459 s2_incr += w4;
00460 s2yz += w5;
00461 }
00462 distsx += s1xx;
00463 s1xx += w1;
00464 s1xy += w2;
00465 s1xz += w3;
00466 }
00467 if (n_density_contributions > 0) {
00468 result.push_back( density_sum
00469 / static_cast<DataType>(n_density_contributions));
00470 }
00471 else {
00472 result.push_back(0);
00473 }
00474 }
00475 return result;
00476 }
00477
00478 }}
00479
00480 #endif // CCTBX_MAPTBX_AVERAGE_DENSITIES_H
