1 // SASAEVAL.H : write a short description here... 2 3 // Copyright (C) 2000 Tommi Hassinen. 4 5 // This package is free software; you can redistribute it and/or modify 6 // it under the terms of the GNU General Public License as published by 7 // the Free Software Foundation; either version 2 of the License, or 8 // (at your option) any later version. 9 10 // This package is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU General Public License for more details. 14 15 // You should have received a copy of the GNU General Public License 16 // along with this package; if not, write to the Free Software 17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 18 19 /*################################################################################################*/ 20 21 //#include "libghemicalconfig2.h" 22 23 #ifndef SASAEVAL_H 24 #define SASAEVAL_H 25 26 #include "typedef.h" 27 #include "notice.h" 28 29 #include <cstdlib> 30 using namespace std; 31 32 struct cg_nbt3_nl; // SASA neighbor list. 33 34 struct cg_nbt3_nd; // SASA neighbor data. 35 struct cg_nbt3_ipd; // SASA intersection point data. 36 37 struct cg_nbt3_coi; 38 struct cg_nbt3_ips; 39 struct cg_nbt3_arc; 40 41 #define SIZE_NLI 200 42 #define SIZE_NT 100 43 #define SIZE_COI 100 44 #define SIZE_IPD 50 45 #define SIZE_IPS 100 46 #define SIZE_ARC 100 47 48 #define INDEX_FLAG 0x8000000 // index of the point 49 #define ORDER_FLAG 0x4000000 // 0 = starting point, 1 = ending point 50 51 #define FLAG_MASK ~(INDEX_FLAG | ORDER_FLAG) 52 53 /*################################################################################################*/ 54 55 struct cg_nbt3_nl ///< SASA neighbor list. 56 { 57 i32s index_count; 58 i32s * index; 59 }; 60 61 struct cg_nbt3_nd ///< SASA neighbor data. 62 { 63 i32s index; 64 f64 distance; 65 66 // these are sorted in reverse order, from large to small... 67 // these are sorted in reverse order, from large to small... 68 // these are sorted in reverse order, from large to small... 69 70 bool operator<(const cg_nbt3_nd & p1) const 71 { 72 return (distance > p1.distance); 73 } 74 }; 75 76 struct cg_nbt3_ipd ///< SASA intersection point data. 77 { 78 f64 angle; 79 i32u ipdata; 80 81 bool operator<(const cg_nbt3_ipd & p1) const 82 { 83 return (angle < p1.angle); 84 } 85 }; 86 87 struct cg_nbt3_coi ///< SASA circle of intersection. 88 { 89 i32s index; 90 91 i32s ipd_count; 92 cg_nbt3_ipd ipdt[SIZE_IPD]; 93 94 f64 refv[3]; 95 96 f64 dist; 97 f64 dv[3]; f64 ddv[3][3]; 98 99 f64 g; f64 dg[3]; 100 f64 ct; f64 dct[3]; 101 102 bool flag; 103 AddIPDcg_nbt3_coi104 void AddIPD(f64 * p1, i32u p2) 105 { 106 ipdt[ipd_count].ipdata = p2; 107 108 if (!ipd_count) 109 { 110 f64 t1a[3]; 111 t1a[0] = dv[0] * p1[0]; 112 t1a[1] = dv[1] * p1[1]; 113 t1a[2] = dv[2] * p1[2]; 114 115 f64 t1b = t1a[0] + t1a[1] + t1a[2]; 116 117 refv[0] = p1[0] - dv[0] * t1b; 118 refv[1] = p1[1] - dv[1] * t1b; 119 refv[2] = p1[2] - dv[2] * t1b; 120 121 f64 t1c = sqrt(refv[0] * refv[0] + refv[1] * refv[1] + refv[2] * refv[2]); 122 refv[0] /= t1c; refv[1] /= t1c; refv[2] /= t1c; 123 124 ipdt[ipd_count].angle = 0.0; 125 } 126 else 127 { 128 f64 t1a[3]; 129 t1a[0] = dv[0] * p1[0]; 130 t1a[1] = dv[1] * p1[1]; 131 t1a[2] = dv[2] * p1[2]; 132 133 f64 t1b = t1a[0] + t1a[1] + t1a[2]; 134 135 f64 t2a[3]; 136 t2a[0] = p1[0] - dv[0] * t1b; 137 t2a[1] = p1[1] - dv[1] * t1b; 138 t2a[2] = p1[2] - dv[2] * t1b; 139 140 f64 t1c = sqrt(t2a[0] * t2a[0] + t2a[1] * t2a[1] + t2a[2] * t2a[2]); 141 t2a[0] /= t1c; t2a[1] /= t1c; t2a[2] /= t1c; 142 143 f64 t1d = refv[0] * t2a[0] + refv[1] * t2a[1] + refv[2] * t2a[2]; 144 if (t1d < -1.0) t1d = -1.0; // domain check... 145 if (t1d > +1.0) t1d = +1.0; // domain check... 146 147 f64 t9a = acos(t1d); 148 149 f64 t3a[3]; 150 t3a[0] = dv[1] * t2a[2] - dv[2] * t2a[1]; 151 t3a[1] = dv[2] * t2a[0] - dv[0] * t2a[2]; 152 t3a[2] = dv[0] * t2a[1] - dv[1] * t2a[0]; 153 154 f64 t9b = refv[0] * t3a[0] + refv[1] * t3a[1] + refv[2] * t3a[2]; 155 156 if (t9b < 0.0) ipdt[ipd_count].angle = -t9a; 157 else ipdt[ipd_count].angle = +t9a; 158 } 159 160 ipd_count++; 161 if (ipd_count >= SIZE_IPD) { assertion_failed(__FILE__, __LINE__, "ipd_count >= SIZE_IPD"); } 162 } 163 }; 164 165 struct cg_nbt3_ips ///< SASA intersection points. 166 { 167 i32s coi[2]; 168 169 f64 ipv[2][3]; 170 f64 dipv[2][2][3][3]; 171 }; 172 173 struct cg_nbt3_arc ///< SASA positively oriented arc. 174 { 175 i32s coi; 176 i32s index[2][2]; 177 178 i32u ipdata[2]; 179 180 f64 tv[2][3]; 181 f64 dtv[2][2][3][3]; 182 183 bool flag; 184 }; 185 186 class sasaeval 187 { 188 protected: 189 190 engine * eng; 191 192 i32s natm_GLOB; 193 i32s natm_loc; 194 195 f64 * radius_GLOB; // pre-reg ; using global indices 196 i32u * index_GLOB_2_LOC; // index global-to-local translation 197 198 //////////////////////////////////////// 199 200 // the rest will work using local (registered) atom indices... 201 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 202 203 i32u * index_l2g; 204 205 f64 * radius1; 206 f64 * radius2; 207 208 i32s * dist1; 209 f64 * dist2; 210 211 cg_nbt3_nl * nl; 212 213 f64 * sasa; 214 f64 * d_sasa; 215 216 public: 217 218 sasaeval(engine *); 219 ~sasaeval(void); 220 221 bool RegisterAtom(i32u, double); 222 void RegisterAtomsFinished(void); 223 224 void HandleNL(i32u, i32u, f64); 225 226 void Evaluate(i32s); 227 }; 228 229 /*################################################################################################*/ 230 231 #endif // SASAEVAL_H 232 233 // eof 234