VinaLC: Parallel Molecular Docking Program

Biochemical and Biophysical Systems Group
VinaLC version: 1.1.2
precalculate.h
Go to the documentation of this file.
1 /*
2 
3  Copyright (c) 2006-2010, The Scripps Research Institute
4 
5  Licensed under the Apache License, Version 2.0 (the "License");
6  you may not use this file except in compliance with the License.
7  You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11  Unless required by applicable law or agreed to in writing, software
12  distributed under the License is distributed on an "AS IS" BASIS,
13  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  See the License for the specific language governing permissions and
15  limitations under the License.
16 
17  Author: Dr. Oleg Trott <ot14@columbia.edu>,
18  The Olson Lab,
19  The Scripps Research Institute
20 
21 */
22 
23 #ifndef VINA_PRECALCULATE_H
24 #define VINA_PRECALCULATE_H
25 
26 #include "scoring_function.h"
27 #include "matrix.h"
28 
29 struct precalculate_element {
30  precalculate_element(sz n, fl factor_) : fast(n, 0), smooth(n, pr(0, 0)), factor(factor_) {}
31  fl eval_fast(fl r2) const {
32  assert(r2 * factor < fast.size());
33  sz i = sz(factor * r2); // r2 is expected < cutoff_sqr, and cutoff_sqr * factor + 1 < n, so no overflow
34  assert(i < fast.size());
35  return fast[i];
36  }
37  pr eval_deriv(fl r2) const {
38  fl r2_factored = factor * r2;
39  assert(r2_factored + 1 < smooth.size());
40  sz i1 = sz(r2_factored);
41  sz i2 = i1 + 1; // r2 is expected < cutoff_sqr, and cutoff_sqr * factor + 1 < n, so no overflow
42  assert(i1 < smooth.size());
43  assert(i2 < smooth.size());
44  fl rem = r2_factored - i1;
45  assert(rem >= -epsilon_fl);
46  assert(rem < 1 + epsilon_fl);
47  const pr& p1 = smooth[i1];
48  const pr& p2 = smooth[i2];
49  fl e = p1.first + rem * (p2.first - p1.first);
50  fl dor = p1.second + rem * (p2.second - p1.second);
51  return pr(e, dor);
52  }
53  void init_from_smooth_fst(const flv& rs) {
54  sz n = smooth.size();
55  VINA_CHECK(rs.size() == n);
56  VINA_CHECK(fast.size() == n);
57  VINA_FOR(i, n) {
58  // calculate dor's
59  fl& dor = smooth[i].second;
60  if(i == 0 || i == n-1)
61  dor = 0;
62  else {
63  fl delta = rs[i+1] - rs[i-1];
64  fl r = rs[i];
65  dor = (smooth[i+1].first - smooth[i-1].first) / (delta * r);
66  }
67  // calculate fast's from smooth.first's
68  fl f1 = smooth[i].first;
69  fl f2 = (i+1 >= n) ? 0 : smooth[i+1].first;
70  fast[i] = (f2 + f1) / 2;
71  }
72  }
73  sz min_smooth_fst() const {
74  sz tmp = 0; // returned if smooth.empty()
75  VINA_FOR_IN(i_inv, smooth) {
76  sz i = smooth.size() - i_inv - 1; // i_inv < smooth.size() => i_inv + 1 <= smooth.size()
77  if(i_inv == 0 || smooth[i].first < smooth[tmp].first)
78  tmp = i;
79  }
80  return tmp;
81  }
82  void widen_smooth_fst(const flv& rs, fl left, fl right) {
83  flv tmp(smooth.size(), 0); // the new smooth[].first's
84  sz min_index = min_smooth_fst();
85  VINA_CHECK(min_index < rs.size()); // won't hold for n == 0
86  VINA_CHECK(rs.size() == smooth.size());
87  fl optimal_r = rs[min_index];
88  VINA_FOR_IN(i, smooth) {
89  fl r = rs[i];
90  if (r < optimal_r - left ) r += left;
91  else if(r > optimal_r + right) r -= right;
92  else r = optimal_r;
93 
94  if(r < 0) r = 0;
95  if(r > rs.back()) r = rs.back();
96 
97  tmp[i] = eval_deriv(sqr(r)).first;
98  }
99  VINA_FOR_IN(i, smooth)
100  smooth[i].first = tmp[i];
101  }
102  void widen(const flv& rs, fl left, fl right) {
103  widen_smooth_fst(rs, left, right);
104  init_from_smooth_fst(rs);
105  }
106  flv fast;
107  prv smooth; // [(e, dor)]
108  fl factor;
109 };
110 
111 struct precalculate {
112  precalculate(const scoring_function& sf, fl v = max_fl, fl factor_ = 32) : // sf should not be discontinuous, even near cutoff, for the sake of the derivatives
113  m_cutoff_sqr(sqr(sf.cutoff())),
114  n(sz(factor_ * m_cutoff_sqr) + 3), // sz(factor * r^2) + 1 <= sz(factor * cutoff_sqr) + 2 <= n-1 < n // see assert below
115  factor(factor_),
116 
117  data(num_atom_types(sf.atom_typing_used()), precalculate_element(n, factor_)),
118  m_atom_typing_used(sf.atom_typing_used()) {
119 
120  VINA_CHECK(factor > epsilon_fl);
121  VINA_CHECK(sz(m_cutoff_sqr*factor) + 1 < n); // cutoff_sqr * factor is the largest float we may end up converting into sz, then 1 can be added to the result
122  VINA_CHECK(m_cutoff_sqr*factor + 1 < n);
123 
124  flv rs = calculate_rs();
125 
126  VINA_FOR(t1, data.dim())
127  VINA_RANGE(t2, t1, data.dim()) {
128  precalculate_element& p = data(t1, t2);
129  // init smooth[].first
130  VINA_FOR_IN(i, p.smooth)
131  p.smooth[i].first = (std::min)(v, sf.eval(t1, t2, rs[i]));
132 
133  // init the rest
134  p.init_from_smooth_fst(rs);
135  }
136  }
137  fl eval_fast(sz type_pair_index, fl r2) const {
138  assert(r2 <= m_cutoff_sqr);
139  return data(type_pair_index).eval_fast(r2);
140  }
141  pr eval_deriv(sz type_pair_index, fl r2) const {
142  assert(r2 <= m_cutoff_sqr);
143  return data(type_pair_index).eval_deriv(r2);
144  }
145  sz index_permissive(sz t1, sz t2) const { return data.index_permissive(t1, t2); }
146  atom_type::t atom_typing_used() const { return m_atom_typing_used; }
147  fl cutoff_sqr() const { return m_cutoff_sqr; }
148  void widen(fl left, fl right) {
149  flv rs = calculate_rs();
150  VINA_FOR(t1, data.dim())
151  VINA_RANGE(t2, t1, data.dim())
152  data(t1, t2).widen(rs, left, right);
153  }
154 private:
155  flv calculate_rs() const {
156  flv tmp(n, 0);
157  VINA_FOR(i, n)
158  tmp[i] = std::sqrt(i / factor);
159  return tmp;
160  }
161  fl m_cutoff_sqr;
162  sz n;
163  fl factor;
164  atom_type::t m_atom_typing_used;
165 
166  triangular_matrix<precalculate_element> data;
167 };
168 
169 #endif

LLNL

LLNL-WEB-613932
LLNL-SM-577132

Lawrence Livermore National Laboratory
7000 East Avenue. Livermore, CA 94550

Operated by Lawrence Livermore National Security, LLC, for the
Department of Energy's National Nuclear Security Administration.

NNSA U.S. DOE LLNS