Merge pull request #64 from dptech-corp/feature/Hbond

Feature/Hbond :add H bond constraint

.gitignore

@@ -11,6 +11,10 @@ unidock/example/screening_test/result
 unidock/example/screening_test/def_ligands.txt
 unidock/example/screening_test/rundock.sh
 unidock/example/screening_test/unidock
+unidock/example/screening_test/fail_result
+unidock/example/screening_test/receptor_grids
+unidock/example/screening_test/unidock_root
+unidock/example/screening_test/test_dock
 unidock_tools/dist
 unidock_tools/dist/*
 unidock_tools/unidock_tools.egg-info

unidock/src/cuda/atom_constants_gpu.cuh

@@ -49,7 +49,11 @@ __device__ __constant__ sz AD_TYPE_CG1 = 27;
 __device__ __constant__ sz AD_TYPE_CG2 = 28;
 __device__ __constant__ sz AD_TYPE_CG3 = 29;
 __device__ __constant__ sz AD_TYPE_W = 30;  // hydrated ligand
-__device__ __constant__ sz AD_TYPE_SIZE = 31;
+__device__ __constant__ sz AD_TYPE_OXA  = 31; // biased protein hydrogen bond acceptor
+__device__ __constant__ sz AD_TYPE_NXA  = 32; // biased protein hydrogen bond acceptor
+__device__ __constant__ sz AD_TYPE_OXD  = 33; // biased protein hydrogen bond donor
+__device__ __constant__ sz AD_TYPE_NXD  = 34; // biased protein hydrogen bond donor
+__device__ __constant__ sz AD_TYPE_SIZE = 35;
 
 // X-Score
 __device__ __constant__ sz XS_TYPE_C_H = 0;
@@ -84,7 +88,11 @@ __device__ __constant__ sz XS_TYPE_C_H_CG3 = 28;
 __device__ __constant__ sz XS_TYPE_C_P_CG3 = 29;
 __device__ __constant__ sz XS_TYPE_G3 = 30;
 __device__ __constant__ sz XS_TYPE_W = 31;  // hydrated ligand
-__device__ __constant__ sz XS_TYPE_SIZE = 32;
+__device__ __constant__  sz XS_TYPE_O_XA    = 32;
+__device__ __constant__  sz XS_TYPE_N_XA    = 33;
+__device__ __constant__  sz XS_TYPE_O_XD	 = 34;
+__device__ __constant__  sz XS_TYPE_N_XD    = 35;
+__device__ __constant__  sz XS_TYPE_SIZE = 36;
 
 // DrugScore-CSD
 __device__ __constant__ sz SY_TYPE_C_3 = 0;
@@ -155,7 +163,11 @@ __device__ atom_kind_gpu atom_kind_data_gpu[] = {
     {"CG1", 2.00000, 0.15000, 0.0, 0.0, -0.00143, 33.51030, 0.77},  // 27
     {"CG2", 2.00000, 0.15000, 0.0, 0.0, -0.00143, 33.51030, 0.77},  // 28
     {"CG3", 2.00000, 0.15000, 0.0, 0.0, -0.00143, 33.51030, 0.77},  // 29
-    {"W", 0.00000, 0.00000, 0.0, 0.0, 0.00000, 0.00000, 0.00}       // 30
+    {"W", 0.00000, 0.00000, 0.0, 0.0, 0.00000, 0.00000, 0.00},       // 30
+    {"OXA", 1.60000, 0.20000, 0.0, 0.0,   -0.00251,   17.15730,   0.73}, //  31
+	{"NXA", 1.75000, 0.16000, 0.0, 0.0,   -0.00162,   22.44930,   0.75}, //  32
+	{"OXD", 1.60000, 0.20000, 0.0, 0.0,   -0.00251,   17.15730,   0.73}, //  33
+	{"NXD", 1.75000, 0.16000, 0.0, 0.0,   -0.00162,   22.44930,   0.75}  //  34
 };
 
 __device__ __constant__ fl metal_solvation_parameter_gpu = -0.00110;
@@ -258,6 +270,14 @@ __device__ __forceinline__ sz ad_type_to_el_type_gpu(sz t) {
             return EL_TYPE_Dummy;
         case AD_TYPE_W:
             return EL_TYPE_Dummy;
+        case AD_TYPE_OXA: 
+            return EL_TYPE_O;
+		case AD_TYPE_NXA: 
+            return EL_TYPE_N;
+		case AD_TYPE_OXD: 
+            return EL_TYPE_O;
+		case AD_TYPE_NXD: 
+            return EL_TYPE_N;
         case AD_TYPE_SIZE:
             return EL_TYPE_SIZE;
         default:
@@ -299,7 +319,11 @@ __device__ __constant__ fl xs_vdw_radii_gpu[] = {
     0.0,  // G1
     0.0,  // G2
     0.0,  // G3
-    0.0   // W
+    0.0,   // W
+    1.7, // O_XA
+	1.8, // N_XA
+	1.7, // O_XD
+	1.8  // N_XD
 };
 
 __device__ __constant__ fl xs_vinardo_vdw_radii_gpu[] = {
@@ -334,7 +358,11 @@ __device__ __constant__ fl xs_vinardo_vdw_radii_gpu[] = {
     0.0,  // G1
     0.0,  // G2
     0.0,  // G3
-    0.0   // W
+    0.0,   // W
+    1.6, // O_XA
+	1.7, // N_XA
+	1.6, // O_XD
+	1.7
 };
 
 __device__ __forceinline__ fl xs_radius_gpu(sz t) {
@@ -366,12 +394,12 @@ __device__ __forceinline__ bool xs_is_hydrophobic_gpu(sz xs) {
 }
 
 __device__ __forceinline__ bool xs_is_acceptor_gpu(sz xs) {
-    return xs == XS_TYPE_N_A || xs == XS_TYPE_N_DA || xs == XS_TYPE_O_A || xs == XS_TYPE_O_DA;
+    return xs == XS_TYPE_N_A || xs == XS_TYPE_N_DA || xs == XS_TYPE_O_A || xs == XS_TYPE_O_DA||xs == XS_TYPE_O_XA ||xs == XS_TYPE_N_XA;
 }
 
 __device__ __forceinline__ bool xs_is_donor_gpu(sz xs) {
     return xs == XS_TYPE_N_D || xs == XS_TYPE_N_DA || xs == XS_TYPE_O_D || xs == XS_TYPE_O_DA
-           || xs == XS_TYPE_Met_D;
+           || xs == XS_TYPE_Met_D||xs == XS_TYPE_O_XD ||xs == XS_TYPE_N_XD;
 }
 
 __device__ __forceinline__ bool xs_donor_acceptor_gpu(sz t1, sz t2) {

unidock/src/cuda/kernel.h

@@ -36,7 +36,7 @@ void check(T result, char const *const func, const char *const file, int const l
 #define MAX_NUM_OF_LIG_PAIRS 4096
 #define MAX_NUM_OF_BFGS_STEPS 64
 #define MAX_NUM_OF_RANDOM_MAP 1000  // not too large (stack overflow!)
-#define GRIDS_SIZE 34               // larger than vina1.1, max(XS_TYPE_SIZE, AD_TYPE_SIZE + 2)
+#define GRIDS_SIZE 37               // larger than vina1.1, max(XS_TYPE_SIZE, AD_TYPE_SIZE + 2)
 
 #define MAX_NUM_OF_GRID_MI 128  // 55
 #define MAX_NUM_OF_GRID_MJ 128  // 55

unidock/src/cuda/precalculate_gpu.cuh

@@ -144,8 +144,11 @@ __device__ __forceinline__ fl vina_non_dir_h_bond_cuda_eval(sz t1, sz t2, fl r,
                                                             fl cutoff) {
     if (r >= cutoff) return 0.0;
     if ((t1 >= XS_TYPE_SIZE) || (t2 >= XS_TYPE_SIZE)) return 0.0;
-    if (xs_h_bond_possible_gpu(t1, t2))
-        return slope_step_gpu(bad, good, r - optimal_distance_gpu(t1, t2));
+    if (xs_h_bond_possible_gpu(t1, t2)){        
+        if  ((t1 >= 32 && t1 <= 35) || (t2 >= 32 && t2 <= 35))
+            return 10.0*slope_step_gpu(bad, good, r - optimal_distance_gpu(t1, t2));
+        else 
+            return slope_step_gpu(bad, good, r - optimal_distance_gpu(t1, t2));}
     return 0.0;
 };
 

unidock/src/lib/ad4cache.cpp

@@ -412,7 +412,7 @@ void ad4cache::read(const std::string& map_prefix) {
             read_ad4_map(p, gds, m_grids[t]);
 
         }  // if file exists
-    }      // map loop
+    }  // map loop
 
     //  elec map
     filename = map_prefix + ".e.map";
@@ -506,9 +506,9 @@ void ad4cache::write(const std::string& out_prefix, const szv& atom_types,
                 VINA_FOR(y, m_grids[t].m_data.dim1()) {
                     VINA_FOR(x, m_grids[t].m_data.dim0()) {
                         out << std::setprecision(4) << m_grids[t].m_data(x, y, z) << "\n";  // slow?
-                    }                                                                       // x
-                }                                                                           // y
-            }                                                                               // z
+                    }  // x
+                }  // y
+            }  // z
         }  // map initialized
-    }      // map atom type
+    }  // map atom type
 }  // cache::write

unidock/src/lib/array3d.h

@@ -41,10 +41,10 @@ template <typename T> class array3d {
 
     friend class boost::serialization::access;
     template <typename Archive> void serialize(Archive& ar, const unsigned version) {
-        ar& m_i;
-        ar& m_j;
-        ar& m_k;
-        ar& m_data;
+        ar & m_i;
+        ar & m_j;
+        ar & m_k;
+        ar & m_data;
     }
 
 public:

unidock/src/lib/atom.h

@@ -34,8 +34,8 @@ struct atom_index {
 private:
     friend class boost::serialization::access;
     template <class Archive> void serialize(Archive& ar, const unsigned version) {
-        ar& i;
-        ar& in_grid;
+        ar & i;
+        ar & in_grid;
     }
 };
 
@@ -54,9 +54,9 @@ struct bond {
 private:
     friend class boost::serialization::access;
     template <class Archive> void serialize(Archive& ar, const unsigned version) {
-        ar& connected_atom_index;
-        ar& length;
-        ar& rotatable;
+        ar & connected_atom_index;
+        ar & length;
+        ar & rotatable;
     }
 };
 
@@ -70,8 +70,8 @@ struct atom : public atom_base {
     friend class boost::serialization::access;
     template <class Archive> void serialize(Archive& ar, const unsigned version) {
         ar& boost::serialization::base_object<atom_base>(*this);
-        ar& coords;
-        ar& bonds;
+        ar & coords;
+        ar & bonds;
     }
 };
 

unidock/src/lib/atom_base.h

@@ -33,7 +33,7 @@ struct atom_base : public atom_type {
     friend class boost::serialization::access;
     template <class Archive> void serialize(Archive& ar, const unsigned version) {
         ar& boost::serialization::base_object<atom_type>(*this);
-        ar& charge;
+        ar & charge;
     }
 };
 

unidock/src/lib/atom_constants.h

@@ -73,8 +73,12 @@ const sz AD_TYPE_CG0 = 26;
 const sz AD_TYPE_CG1 = 27;
 const sz AD_TYPE_CG2 = 28;
 const sz AD_TYPE_CG3 = 29;
-const sz AD_TYPE_W = 30;  // hydrated ligand
-const sz AD_TYPE_SIZE = 31;
+const sz AD_TYPE_W = 30;    // hydrated ligand
+const sz AD_TYPE_OXA = 31;  // biased protein hydrogen bond acceptor
+const sz AD_TYPE_NXA = 32;  // biased protein hydrogen bond acceptor
+const sz AD_TYPE_OXD = 33;  // biased protein hydrogen bond donor
+const sz AD_TYPE_NXD = 34;  // biased protein hydrogen bond donor
+const sz AD_TYPE_SIZE = 35;
 
 // X-Score
 const sz XS_TYPE_C_H = 0;
@@ -109,7 +113,11 @@ const sz XS_TYPE_C_H_CG3 = 28;
 const sz XS_TYPE_C_P_CG3 = 29;
 const sz XS_TYPE_G3 = 30;
 const sz XS_TYPE_W = 31;  // hydrated ligand
-const sz XS_TYPE_SIZE = 32;
+const sz XS_TYPE_O_XA = 32;
+const sz XS_TYPE_N_XA = 33;
+const sz XS_TYPE_O_XD = 34;
+const sz XS_TYPE_N_XD = 35;
+const sz XS_TYPE_SIZE = 36;
 
 // DrugScore-CSD
 const sz SY_TYPE_C_3 = 0;
@@ -178,7 +186,11 @@ const atom_kind atom_kind_data[] = {
     {"CG1", 2.00000, 0.15000, 0.0, 0.0, -0.00143, 33.51030, 0.77},  // 27
     {"CG2", 2.00000, 0.15000, 0.0, 0.0, -0.00143, 33.51030, 0.77},  // 28
     {"CG3", 2.00000, 0.15000, 0.0, 0.0, -0.00143, 33.51030, 0.77},  // 29
-    {"W", 0.00000, 0.00000, 0.0, 0.0, 0.00000, 0.00000, 0.00}       // 30
+    {"W", 0.00000, 0.00000, 0.0, 0.0, 0.00000, 0.00000, 0.00},      // 30
+    {"OXA", 1.60000, 0.20000, 0.0, 0.0, -0.00251, 17.15730, 0.73},  //  31
+    {"NXA", 1.75000, 0.16000, 0.0, 0.0, -0.00162, 22.44930, 0.75},  //  32
+    {"OXD", 1.60000, 0.20000, 0.0, 0.0, -0.00251, 17.15730, 0.73},  //  33
+    {"NXD", 1.75000, 0.16000, 0.0, 0.0, -0.00162, 22.44930, 0.75}   //  34
 };
 
 const fl metal_solvation_parameter = -0.00110;
@@ -281,6 +293,14 @@ inline sz ad_type_to_el_type(sz t) {
             return EL_TYPE_Dummy;
         case AD_TYPE_W:
             return EL_TYPE_Dummy;
+        case AD_TYPE_OXA:
+            return EL_TYPE_O;
+        case AD_TYPE_NXA:
+            return EL_TYPE_N;
+        case AD_TYPE_OXD:
+            return EL_TYPE_O;
+        case AD_TYPE_NXD:
+            return EL_TYPE_N;
         case AD_TYPE_SIZE:
             return EL_TYPE_SIZE;
         default:
@@ -322,7 +342,11 @@ const fl xs_vdw_radii[] = {
     0.0,  // G1
     0.0,  // G2
     0.0,  // G3
-    0.0   // W
+    0.0,  // W
+    1.7,  // O_XA
+    1.8,  // N_XA
+    1.7,  // O_XD
+    1.8   // N_XD
 };
 
 const fl xs_vinardo_vdw_radii[] = {
@@ -357,7 +381,11 @@ const fl xs_vinardo_vdw_radii[] = {
     0.0,  // G1
     0.0,  // G2
     0.0,  // G3
-    0.0   // W
+    0.0,  // W
+    1.6,  // O_XA
+    1.7,  // N_XA
+    1.6,  // O_XD
+    1.7   // N_XD
 };
 
 inline fl xs_radius(sz t) {
@@ -388,12 +416,13 @@ inline bool xs_is_hydrophobic(sz xs) {
 }
 
 inline bool xs_is_acceptor(sz xs) {
-    return xs == XS_TYPE_N_A || xs == XS_TYPE_N_DA || xs == XS_TYPE_O_A || xs == XS_TYPE_O_DA;
+    return xs == XS_TYPE_N_A || xs == XS_TYPE_N_DA || xs == XS_TYPE_O_A || xs == XS_TYPE_O_DA
+           || xs == XS_TYPE_O_XA || xs == XS_TYPE_N_XA;
 }
 
 inline bool xs_is_donor(sz xs) {
     return xs == XS_TYPE_N_D || xs == XS_TYPE_N_DA || xs == XS_TYPE_O_D || xs == XS_TYPE_O_DA
-           || xs == XS_TYPE_Met_D;
+           || xs == XS_TYPE_Met_D || xs == XS_TYPE_N_XD || xs == XS_TYPE_O_XD;
 }
 
 inline bool xs_donor_acceptor(sz t1, sz t2) { return xs_is_donor(t1) && xs_is_acceptor(t2); }

unidock/src/lib/atom_type.h

@@ -73,10 +73,10 @@ struct atom_type {
 private:
     friend class boost::serialization::access;
     template <class Archive> void serialize(Archive& ar, const unsigned version) {
-        ar& el;
-        ar& ad;
-        ar& xs;
-        ar& sy;
+        ar & el;
+        ar & ad;
+        ar & xs;
+        ar & sy;
     }
 };
 

unidock/src/lib/cache.cpp

@@ -68,6 +68,14 @@ std::string convert_XS_to_string(sz t) {
             return "Met_D";
         case XS_TYPE_W:
             return "W";
+        case XS_TYPE_O_XA:
+            return "O_XA";
+        case XS_TYPE_N_XA:
+            return "N_XA";
+        case XS_TYPE_O_XD:
+            return "O_XD";
+        case XS_TYPE_N_XD:
+            return "N_XD";
         default:
             VINA_CHECK(false);
     }
@@ -445,12 +453,12 @@ void cache::write(const std::string& out_prefix, const szv& atom_types,
                         VINA_FOR(x, m_grids[t].m_data.dim0()) {
                             out << std::setprecision(4) << m_grids[t].m_data(x, y, z)
                                 << "\n";  // slow?
-                        }                 // x
-                    }                     // y
-                }                         // z
-            }                             // even voxels
-        }                                 // map initialized
-    }                                     // map atom type
+                        }  // x
+                    }  // y
+                }  // z
+            }  // even voxels
+        }  // map initialized
+    }  // map atom type
 }  // cache::write
 
 float cache::get_slope() const { return this->m_slope; }

unidock/src/lib/common.h

@@ -115,7 +115,7 @@ struct vec {
 
 private:
     friend class boost::serialization::access;
-    template <class Archive> void serialize(Archive& ar, const unsigned version) { ar& data; }
+    template <class Archive> void serialize(Archive& ar, const unsigned version) { ar & data; }
 };
 
 inline vec operator*(fl s, const vec& v) { return vec(s * v[0], s * v[1], s * v[2]); }