00001 #ifndef CCTBX_MAPTBX_REAL_SPACE_TARGET_AND_GRADIENTS_H
00002 #define CCTBX_MAPTBX_REAL_SPACE_TARGET_AND_GRADIENTS_H
00003
00004 #include <cctbx/maptbx/eight_point_interpolation.h>
00005
00006 namespace cctbx { namespace maptbx {
00007
00008 class target_and_gradients {
00009 public:
00010 target_and_gradients(
00011 uctbx::unit_cell const& unit_cell,
00012 af::const_ref<double, af::c_grid_padded<3> > const& map_target,
00013 af::const_ref<double, af::c_grid_padded<3> > const& map_current,
00014 double const& step,
00015 af::const_ref<scitbx::vec3<double> > const& sites_frac)
00016 {
00017 int nx = static_cast<int>(map_target.accessor().focus()[0]);
00018 int ny = static_cast<int>(map_target.accessor().focus()[1]);
00019 int nz = static_cast<int>(map_target.accessor().focus()[2]);
00020 af::versa<double, af::c_grid_padded<3> > diff_density_array;
00021 diff_density_array.resize(af::c_grid_padded<3>(nx,ny,nz), 0);
00022 int x1box=0 ;
00023 int x2box=nx-1;
00024 int y1box=0 ;
00025 int y2box=ny-1;
00026 int z1box=0 ;
00027 int z2box=nz-1;
00028 target_ = 0;
00029 double* diff_density_array_ = diff_density_array.begin();
00030 for(int kx = x1box; kx <= x2box; kx++) {
00031 int mx = scitbx::math::mod_positive(kx,nx);
00032 int mxny = mx*ny;
00033 for(int ky = y1box; ky <= y2box; ky++) {
00034 int my = scitbx::math::mod_positive(ky,ny);
00035 int mxnypmynz = (mxny+my)*nz;
00036 for(int kz = z1box; kz <= z2box; kz++) {
00037 int mz = scitbx::math::mod_positive(kz,nz);
00038 double diff = map_target[mxnypmynz+mz]-map_current[mxnypmynz+mz];
00039 target_ += (diff * diff);
00040 diff_density_array_[mxnypmynz+mz] = -2 * diff;
00041 }
00042 }
00043 }
00044
00045 af::const_ref<double, af::c_grid_padded<3> > diffd =
00046 diff_density_array.const_ref();
00047 cctbx::cartesian<> step_x = cctbx::cartesian<>(step,0,0);
00048 cctbx::cartesian<> step_y = cctbx::cartesian<>(0,step,0);
00049 cctbx::cartesian<> step_z = cctbx::cartesian<>(0,0,step);
00050 double two_step = 2*step;
00051 gradients_.resize(sites_frac.size(), scitbx::vec3<double>(0,0,0));
00052 for(std::size_t i_site=0;i_site<sites_frac.size();i_site++) {
00053 cctbx::fractional<> const& site_frac = sites_frac[i_site];
00054 cctbx::cartesian<> site_cart = unit_cell.orthogonalize(site_frac);
00055 cctbx::fractional<> sxp = unit_cell.fractionalize(site_cart+step_x);
00056 cctbx::fractional<> sxm = unit_cell.fractionalize(site_cart-step_x);
00057 cctbx::fractional<> syp = unit_cell.fractionalize(site_cart+step_y);
00058 cctbx::fractional<> sym = unit_cell.fractionalize(site_cart-step_y);
00059 cctbx::fractional<> szp = unit_cell.fractionalize(site_cart+step_z);
00060 cctbx::fractional<> szm = unit_cell.fractionalize(site_cart-step_z);
00061 double gx = (eight_point_interpolation(diffd, sxp) -
00062 eight_point_interpolation(diffd, sxm)) / two_step;
00063 double gy = (eight_point_interpolation(diffd, syp) -
00064 eight_point_interpolation(diffd, sym)) / two_step;
00065 double gz = (eight_point_interpolation(diffd, szp) -
00066 eight_point_interpolation(diffd, szm)) / two_step;
00067 gradients_[i_site] = scitbx::vec3<double>(gx,gy,gz);
00068 }
00069 }
00070
00071 double target() { return target_; }
00072 af::shared<scitbx::vec3<double> > gradients() { return gradients_; }
00073
00074 protected:
00075 double target_;
00076 af::shared<scitbx::vec3<double> > gradients_;
00077 };
00078
00079
00080 }}
00081
00082 #endif // GUARD
