00001 #ifndef REFERENCE_DENSITY_KERNELS_H_
00002 #define REFERENCE_DENSITY_KERNELS_H_
00003
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 #include <ost/mol/mm/density_kernels.hh>
00036 #include <openmm/Platform.h>
00037 #include <vector>
00038 #include <set>
00039 #include <limits>
00040
00041 namespace ost{ namespace mol{ namespace mm{
00042
00046 class ReferenceCalcDensityForceKernel : public CalcDensityForceKernel {
00047 public:
00048 ReferenceCalcDensityForceKernel(std::string name, const OpenMM::Platform& platform) : CalcDensityForceKernel(name, platform), data_(NULL) {
00049 }
00050
00051 ~ReferenceCalcDensityForceKernel();
00052
00059 void initialize(const OpenMM::System& system, const DensityForce& force);
00068 double execute(OpenMM::ContextImpl& context, bool includeForces, bool includeEnergy);
00075 void copyParametersToContext(OpenMM::ContextImpl& context, const DensityForce& force);
00076 private:
00077
00078 inline Real GetValue(int a, int b, int c){
00079 return data_[a * idx_helper_one_ + b * idx_helper_two_ + c];
00080 }
00081
00082 inline Real GetIntpolValue(Real x, Real y, Real z){
00083
00084
00085 Real dx = x - x_origin_;
00086 Real dy = y - y_origin_;
00087 Real dz = z - z_origin_;
00088
00089 int x_bin = std::floor(dx * one_over_x_sampling_);
00090 int y_bin = std::floor(dy * one_over_y_sampling_);
00091 int z_bin = std::floor(dz * one_over_z_sampling_);
00092
00093 if(x_bin < 0 || x_bin >= (x_extent_-1) ||
00094 y_bin < 0 || y_bin >= (y_extent_-1) ||
00095 z_bin < 0 || z_bin >= (z_extent_-1)) return 0.0;
00096
00097
00098 dx = (dx - x_bin * x_sampling_) * one_over_x_sampling_;
00099 dy = (dy - y_bin * y_sampling_) * one_over_y_sampling_;
00100 dz = (dz - z_bin * z_sampling_) * one_over_z_sampling_;
00101
00102 intpol_values_[0] = this->GetValue(x_bin,y_bin,z_bin);
00103 intpol_values_[1] = this->GetValue(x_bin+1,y_bin,z_bin);
00104 intpol_values_[2] = this->GetValue(x_bin,y_bin+1,z_bin);
00105 intpol_values_[3] = this->GetValue(x_bin+1,y_bin+1,z_bin);
00106 intpol_values_[4] = this->GetValue(x_bin,y_bin,z_bin+1);
00107 intpol_values_[5] = this->GetValue(x_bin+1,y_bin,z_bin+1);
00108 intpol_values_[6] = this->GetValue(x_bin,y_bin+1,z_bin+1);
00109 intpol_values_[7] = this->GetValue(x_bin+1,y_bin+1,z_bin+1);
00110
00111
00112 Real f11 = (1.0 - dx) * intpol_values_[0] + dx * intpol_values_[1];
00113 Real f12 = (1.0 - dx) * intpol_values_[2] + dx * intpol_values_[3];
00114 Real f21 = (1.0 - dx) * intpol_values_[4] + dx * intpol_values_[5];
00115 Real f22 = (1.0 - dx) * intpol_values_[6] + dx * intpol_values_[7];
00116
00117 Real f1 = (1.0 - dy) * f11 + dy * f12;
00118 Real f2 = (1.0 - dy) * f21 + dy * f22;
00119
00120 return (1.0 - dz) * f1 + dz * f2;
00121 }
00122
00123
00124 Real* data_;
00125 Real intpol_values_[8];
00126 int idx_helper_one_;
00127 int idx_helper_two_;
00128 Real x_sampling_;
00129 Real y_sampling_;
00130 Real z_sampling_;
00131 Real one_over_x_sampling_;
00132 Real one_over_y_sampling_;
00133 Real one_over_z_sampling_;
00134 Real half_x_sampling_;
00135 Real half_y_sampling_;
00136 Real half_z_sampling_;
00137 Real x_origin_;
00138 Real y_origin_;
00139 Real z_origin_;
00140 int x_extent_;
00141 int y_extent_;
00142 int z_extent_;
00143
00144
00145 Real resolution_;
00146 Real s_;
00147 Real one_s_;
00148 Real square_one_s_;
00149 Real one_sqrt_2_pi_s_;
00150 Real padding_dist_;
00151 Real scaling_;
00152 std::vector<Real> density_values_;
00153 std::vector<Real> body_values_;
00154
00155
00156 int num_bodies_;
00157 int num_particles_;
00158 std::vector<float> particle_masses_;
00159 std::vector<std::vector<int> > bodies_;
00160 };
00161
00162 }}}
00163
00164 #endif