RbtVdwSF.cxx 18.6 KB
Newer Older
sruizcarmona's avatar
sruizcarmona committed
1
/***********************************************************************
2 3 4 5 6 7 8 9 10 11
 * The rDock program was developed from 1998 - 2006 by the software team
 * at RiboTargets (subsequently Vernalis (R&D) Ltd).
 * In 2006, the software was licensed to the University of York for
 * maintenance and distribution.
 * In 2012, Vernalis and the University of York agreed to release the
 * program as Open Source software.
 * This version is licensed under GNU-LGPL version 3.0 with support from
 * the University of Barcelona.
 * http://rdock.sourceforge.net/
 ***********************************************************************/
pschmidtke's avatar
pschmidtke committed
12 13 14

#include "RbtVdwSF.h"
#include "RbtModel.h"
15
#include <spdlog/spdlog.h>
pschmidtke's avatar
pschmidtke committed
16

17
// Static data members
18 19 20 21 22 23
std::string RbtVdwSF::_CT("RbtVdwSF");
std::string RbtVdwSF::_USE_4_8("USE_4_8");
std::string RbtVdwSF::_USE_TRIPOS("USE_TRIPOS");
std::string RbtVdwSF::_RMAX("RMAX");
std::string RbtVdwSF::_ECUT("ECUT");
std::string RbtVdwSF::_E0("E0");
pschmidtke's avatar
pschmidtke committed
24

25 26 27
RbtVdwSF::RbtVdwSF()
    : m_use_4_8(true), m_use_tripos(false), m_rmax(1.5), m_ecut(1.0),
      m_e0(1.5) {
28
  spdlog::trace("RbtVdwSF default constructor");
29 30 31 32 33 34 35 36
  // Add parameters
  AddParameter(_USE_4_8, m_use_4_8);
  AddParameter(_USE_TRIPOS, m_use_tripos);
  AddParameter(_RMAX, m_rmax);
  AddParameter(_ECUT, m_ecut);
  AddParameter(_E0, m_e0);
  m_spVdwSource = RbtParameterFileSourcePtr(new RbtParameterFileSource(
      Rbt::GetRbtFileName("data/sf", "Tripos52_vdw.prm")));
pschmidtke's avatar
pschmidtke committed
37 38 39 40 41
  Setup();
  _RBTOBJECTCOUNTER_CONSTR_(_CT);
}

RbtVdwSF::~RbtVdwSF() {
42
  spdlog::trace("RbtVdwSF destructor");
pschmidtke's avatar
pschmidtke committed
43 44 45
  _RBTOBJECTCOUNTER_DESTR_(_CT);
}

46 47 48
// As this has a virtual base class we need a separate OwnParameterUpdated
// which can be called by concrete subclass ParameterUpdated methods
// See Stroustrup C++ 3rd edition, p395, on programming virtual base classes
49
void RbtVdwSF::OwnParameterUpdated(const std::string &strName) {
50
  // DM 25 Oct 2000 - heavily used params
pschmidtke's avatar
pschmidtke committed
51
  if (strName == _USE_4_8) {
52
    m_use_4_8 = GetParameter(_USE_4_8);
pschmidtke's avatar
pschmidtke committed
53
    Setup();
54
  } else if (strName == _USE_TRIPOS) {
pschmidtke's avatar
pschmidtke committed
55 56
    m_use_tripos = GetParameter(_USE_TRIPOS);
    Setup();
57
  } else if (strName == _RMAX) {
pschmidtke's avatar
pschmidtke committed
58 59
    m_rmax = GetParameter(_RMAX);
    Setup();
60
  } else if (strName == _ECUT) {
pschmidtke's avatar
pschmidtke committed
61 62
    m_ecut = GetParameter(_ECUT);
    SetupCloseRange();
63
  } else if (strName == _E0) {
pschmidtke's avatar
pschmidtke committed
64 65 66 67 68
    m_e0 = GetParameter(_E0);
    SetupCloseRange();
  }
}

69 70
// Used by subclasses to calculate vdW potential between pAtom and all atoms in
// atomList
71 72 73
double RbtVdwSF::VdwScore(const RbtAtom *pAtom,
                          const RbtAtomRList &atomList) const {
  double score = 0.0;
pschmidtke's avatar
pschmidtke committed
74 75 76 77
  if (atomList.empty()) {
    return score;
  }

78 79 80
  const RbtCoord &c1 = pAtom->GetCoords();
  // Get the iterator into the appropriate row of the vdw table for this atom
  // type
pschmidtke's avatar
pschmidtke committed
81 82 83
  RbtTriposAtomType::eType type1 = pAtom->GetTriposType();
  RbtVdwTableConstIter iter1 = m_vdwTable.begin() + type1;

84
  // 4-8 potential, never annotated
pschmidtke's avatar
pschmidtke committed
85
  if (m_use_4_8) {
86 87 88
    for (RbtAtomRListConstIter iter = atomList.begin(); iter != atomList.end();
         iter++) {
      const RbtCoord &c2 = (*iter)->GetCoords();
89
      double R_sq = Rbt::Length2(c1, c2); // Distance squared
pschmidtke's avatar
pschmidtke committed
90
      RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
91
      // iter2 points to the vdw params for this atom type pair
pschmidtke's avatar
pschmidtke committed
92
      RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
93
      double s = f4_8(R_sq, *iter2);
94
      spdlog::debug("4-8vdw {}  weight: {} score {} score*w {}", (*iter)->GetName(), (*iter)->GetReweight(), s, s * (*iter)->GetReweight());
pschmidtke's avatar
pschmidtke committed
95 96 97
      score += s;
    }
  }
98
  // 6-12 with annotation
pschmidtke's avatar
pschmidtke committed
99
  else if (isAnnotationEnabled()) {
100 101 102
    for (RbtAtomRListConstIter iter = atomList.begin(); iter != atomList.end();
         iter++) {
      const RbtCoord &c2 = (*iter)->GetCoords();
103
      double R_sq = Rbt::Length2(c1, c2); // Distance squared
pschmidtke's avatar
pschmidtke committed
104
      RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
105
      // iter2 points to the vdw params for this atom type pair
pschmidtke's avatar
pschmidtke committed
106
      RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
107
      double s = f6_12(R_sq, *iter2);
108
      spdlog::debug("6-12vdw {}  weight: {} score {} score*w {}", (*iter)->GetName(), (*iter)->GetReweight(), s, s * (*iter)->GetReweight());
pschmidtke's avatar
pschmidtke committed
109
      if (s != 0.0) {
110 111
        score += s;
        RbtAnnotationPtr spAnnotation(
112
            new RbtAnnotation(pAtom, *iter, std::sqrt(R_sq), s));
113
        AddAnnotation(spAnnotation);
pschmidtke's avatar
pschmidtke committed
114 115 116
      }
    }
  }
117
  // 6-12 without annotation
pschmidtke's avatar
pschmidtke committed
118
  else {
119 120 121
    for (RbtAtomRListConstIter iter = atomList.begin(); iter != atomList.end();
         iter++) {
      const RbtCoord &c2 = (*iter)->GetCoords();
122
      double R_sq = Rbt::Length2(c1, c2); // Distance squared
pschmidtke's avatar
pschmidtke committed
123
      RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
124
      // iter2 points to the vdw params for this atom type pair
pschmidtke's avatar
pschmidtke committed
125
      RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
126
      double s = f6_12(R_sq, *iter2);
127
      spdlog::debug("6-12vdw {}  weight: {} score {} score*w {}", (*iter)->GetName(), (*iter)->GetReweight(), s, s * (*iter)->GetReweight());
pschmidtke's avatar
pschmidtke committed
128 129 130 131 132 133
      score += s;
    }
  }
  return score;
}

134
// As above, but score is calculated only between enabled atoms
135 136 137
double RbtVdwSF::VdwScoreEnabledOnly(const RbtAtom *pAtom,
                                     const RbtAtomRList &atomList) const {
  double score = 0.0;
pschmidtke's avatar
pschmidtke committed
138 139 140 141
  if (!pAtom->GetEnabled() || atomList.empty()) {
    return score;
  }

142 143 144
  const RbtCoord &c1 = pAtom->GetCoords();
  // Get the iterator into the appropriate row of the vdw table for this atom
  // type
pschmidtke's avatar
pschmidtke committed
145 146 147
  RbtTriposAtomType::eType type1 = pAtom->GetTriposType();
  RbtVdwTableConstIter iter1 = m_vdwTable.begin() + type1;

148
  // 4-8 potential, never annotated
pschmidtke's avatar
pschmidtke committed
149
  if (m_use_4_8) {
150 151
    for (RbtAtomRListConstIter iter = atomList.begin(); iter != atomList.end();
         iter++) {
pschmidtke's avatar
pschmidtke committed
152
      if ((*iter)->GetEnabled()) {
153
        const RbtCoord &c2 = (*iter)->GetCoords();
154
        double R_sq = Rbt::Length2(c1, c2); // Distance squared
155 156 157
        RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
        // iter2 points to the vdw params for this atom type pair
        RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
158
        double s = f4_8(R_sq, *iter2);
159
        score += s;
pschmidtke's avatar
pschmidtke committed
160 161 162
      }
    }
  }
163
  // 6-12 with annotation
pschmidtke's avatar
pschmidtke committed
164
  else if (isAnnotationEnabled()) {
165 166
    for (RbtAtomRListConstIter iter = atomList.begin(); iter != atomList.end();
         iter++) {
pschmidtke's avatar
pschmidtke committed
167
      if ((*iter)->GetEnabled()) {
168
        const RbtCoord &c2 = (*iter)->GetCoords();
169
        double R_sq = Rbt::Length2(c1, c2); // Distance squared
170 171 172
        RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
        // iter2 points to the vdw params for this atom type pair
        RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
173
        double s = f6_12(R_sq, *iter2);
174 175 176
        if (s != 0.0) {
          score += s;
          RbtAnnotationPtr spAnnotation(
177
              new RbtAnnotation(pAtom, *iter, std::sqrt(R_sq), s));
178 179
          AddAnnotation(spAnnotation);
        }
pschmidtke's avatar
pschmidtke committed
180 181 182
      }
    }
  }
183
  // 6-12 without annotation
pschmidtke's avatar
pschmidtke committed
184
  else {
185 186
    for (RbtAtomRListConstIter iter = atomList.begin(); iter != atomList.end();
         iter++) {
pschmidtke's avatar
pschmidtke committed
187
      if ((*iter)->GetEnabled()) {
188
        const RbtCoord &c2 = (*iter)->GetCoords();
189
        double R_sq = Rbt::Length2(c1, c2); // Distance squared
190 191 192
        RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
        // iter2 points to the vdw params for this atom type pair
        RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
193
        double s = f6_12(R_sq, *iter2);
194
        score += s;
pschmidtke's avatar
pschmidtke committed
195 196 197 198 199 200
      }
    }
  }
  return score;
}

201 202 203
// XB This is the old  VdwScore, without reweighting factors
// RbtDouble RbtVdwSF::VdwScoreIntra(const RbtAtom* pAtom, const RbtAtomRList&
// atomList) const {
204 205 206 207 208 209
//  RbtDouble score = 0.0;
//  if (atomList.empty()) {
//    return score;
//  }
//
//  const RbtCoord& c1 = pAtom->GetCoords();
210 211
//  //Get the iterator into the appropriate row of the vdw table for this atom
//  type RbtTriposAtomType::eType type1 = pAtom->GetTriposType();
212 213 214 215
//  RbtVdwTableConstIter iter1 = m_vdwTable.begin() + type1;
//
//  //4-8 potential, never annotated
//  if (m_use_4_8) {
216 217
//    for (RbtAtomRListConstIter iter = atomList.begin(); iter !=
//    atomList.end(); iter++) {
218 219 220 221 222 223 224 225 226 227 228
//      const RbtCoord& c2 = (*iter)->GetCoords();
//      RbtDouble R_sq = Rbt::Length2(c1,c2);//Distance squared
//      RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
//      //iter2 points to the vdw params for this atom type pair
//      RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
//      RbtDouble s = f4_8(R_sq,*iter2);
//      score += s;
//    }
//  }
//  //6-12 with annotation
//  else if (isAnnotationEnabled()) {
229 230
//    for (RbtAtomRListConstIter iter = atomList.begin(); iter !=
//    atomList.end(); iter++) {
231 232 233 234 235 236 237 238
//      const RbtCoord& c2 = (*iter)->GetCoords();
//      RbtDouble R_sq = Rbt::Length2(c1,c2);//Distance squared
//      RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
//      //iter2 points to the vdw params for this atom type pair
//      RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
//      RbtDouble s = f6_12(R_sq,*iter2);
//      if (s != 0.0) {
//	score += s;
239
//	RbtAnnotationPtr spAnnotation(new
240 241
// RbtAnnotation(pAtom,*iter,std::sqrt(R_sq),s));
// AddAnnotation(spAnnotation);
242 243 244 245 246
//      }
//    }
//  }
//  //6-12 without annotation
//  else {
247 248
//    for (RbtAtomRListConstIter iter = atomList.begin(); iter !=
//    atomList.end(); iter++) {
249 250 251 252 253 254
//      const RbtCoord& c2 = (*iter)->GetCoords();
//      RbtDouble R_sq = Rbt::Length2(c1,c2);//Distance squared
//      RbtTriposAtomType::eType type2 = (*iter)->GetTriposType();
//      //iter2 points to the vdw params for this atom type pair
//      RbtVdwRowConstIter iter2 = (*iter1).begin() + type2;
//      RbtDouble s = f6_12(R_sq,*iter2);
255
//     score += s;
256 257 258 259
//    }
//  }
//  return score;
//}
260
double RbtVdwSF::MaxVdwRange(const RbtAtom *pAtom) const {
pschmidtke's avatar
pschmidtke committed
261 262 263
  return m_maxRange[pAtom->GetTriposType()];
}

264
double RbtVdwSF::MaxVdwRange(RbtTriposAtomType::eType t) const {
pschmidtke's avatar
pschmidtke committed
265 266 267
  return m_maxRange[t];
}

268
// Initialise m_vdwTable with appropriate params for each atom type pair
pschmidtke's avatar
pschmidtke committed
269 270
void RbtVdwSF::Setup() {
  RbtTriposAtomType triposType;
271 272
  spdlog::trace("RbtVdwSF::Setup()");
  spdlog::debug("TYPE1,TYPE2,A,B,RMIN,KIJ,RMA");
273 274
  // Dummy read to force parsing of file, otherwise the first SetSection is
  // overridden
pschmidtke's avatar
pschmidtke committed
275
  RbtStringList secList = m_spVdwSource->GetSectionList();
276 277 278 279 280 281
  std::string _R("R");
  std::string _K("K");
  std::string _IP("IP");
  std::string _POL("POL");
  std::string _ISHBD("isHBD");
  std::string _ISHBA("isHBA");
282 283 284
  m_vdwTable = RbtVdwTable(RbtTriposAtomType::MAXTYPES,
                           RbtVdwRow(RbtTriposAtomType::MAXTYPES));
  m_maxRange = RbtDoubleList(RbtTriposAtomType::MAXTYPES, 0.0);
285
  for (int i = RbtTriposAtomType::UNDEFINED; i < RbtTriposAtomType::MAXTYPES;
286 287
       i++) {
    // Read the params for atom type i
288
    std::string stri = triposType.Type2Str(RbtTriposAtomType::eType(i));
pschmidtke's avatar
pschmidtke committed
289
    m_spVdwSource->SetSection(stri);
290 291 292 293 294 295 296
    double Ri = m_spVdwSource->GetParameterValue(_R); // vdw radius
    double Ki = m_spVdwSource->GetParameterValue(_K); // Tripos 5.2 well
                                                      // depth
    bool hasIPi = m_spVdwSource->isParameterPresent(_IP);
    bool isHBDi = m_spVdwSource->isParameterPresent(_ISHBD);
    bool isHBAi = m_spVdwSource->isParameterPresent(_ISHBA);
    double Ii, alphai;
pschmidtke's avatar
pschmidtke committed
297
    if (hasIPi) {
298 299
      Ii = m_spVdwSource->GetParameterValue(_IP); // Ionisation potential, I
      alphai = m_spVdwSource->GetParameterValue(_POL); // Polarisability, alpha
pschmidtke's avatar
pschmidtke committed
300
    }
301 302
    // Symmetric matrix, so only need to generate one half (i.e.
    // m_vdwTable[i][j] = m_vdwTable[j][i])
303
    for (int j = i; j < RbtTriposAtomType::MAXTYPES; j++) {
304
      // Read the params for atom type j
305
      std::string strj = triposType.Type2Str(RbtTriposAtomType::eType(j));
pschmidtke's avatar
pschmidtke committed
306
      m_spVdwSource->SetSection(strj);
307 308 309 310 311 312
      double Rj = m_spVdwSource->GetParameterValue(_R); // vdw radius
      double Kj = m_spVdwSource->GetParameterValue(_K); // Tripos 5.2 well depth
      bool hasIPj = m_spVdwSource->isParameterPresent(_IP);
      bool isHBDj = m_spVdwSource->isParameterPresent(_ISHBD);
      bool isHBAj = m_spVdwSource->isParameterPresent(_ISHBA);
      double Ij, alphaj;
pschmidtke's avatar
pschmidtke committed
313
      if (hasIPj) {
314 315 316
        Ij = m_spVdwSource->GetParameterValue(_IP); // Ionisation potential, I
        alphaj = m_spVdwSource->GetParameterValue(_POL); // Polarisability,
                                                         // alpha
pschmidtke's avatar
pschmidtke committed
317 318
      }

319
      // Create the vdw params for this atom type pair
pschmidtke's avatar
pschmidtke committed
320
      vdwprms prms;
321
      prms.rmin = Ri + Rj; // rmin is the sum of Ri and Rj
322
      double rmax = prms.rmin * m_rmax;
323 324 325 326 327 328 329 330
      m_maxRange[i] = std::max(m_maxRange[i],
                               rmax); // Keep track of max range for atom type i
      m_maxRange[j] = std::max(m_maxRange[j],
                               rmax); // Keep track of max range for atom type j
      prms.rmax_sq = rmax * rmax;     // Max range**2
      // EITHER: Well depth is zero between donor Hs and acceptors
      if ((isHBDi && isHBAj) || (isHBDj && isHBAi)) {
        prms.kij = 0.0;
pschmidtke's avatar
pschmidtke committed
331
      }
332 333
      // OR, use Standard L-J combination rules for well depth
      // switch to this mode if either atom is missing IP values
pschmidtke's avatar
pschmidtke committed
334
      else if (m_use_tripos || !hasIPi || !hasIPj) {
335
        prms.kij = std::sqrt(Ki * Kj);
pschmidtke's avatar
pschmidtke committed
336
      }
337
      // OR, use GOLD rules based on ionisation potential and polarisability
pschmidtke's avatar
pschmidtke committed
338
      else {
339
        double D = 0.345 * Ii * Ij * alphai * alphaj / (Ii + Ij);
340
        double C = 0.5 * D * std::pow(prms.rmin, 6);
341
        prms.kij = D * D / (4.0 * C);
pschmidtke's avatar
pschmidtke committed
342
      }
343 344
      // Having got the well depth, we can now generate A and B for either 4-8
      // or 6-12
345 346
      double rmin_pwr =
          (m_use_4_8) ? std::pow(prms.rmin, 4) : std::pow(prms.rmin, 6);
pschmidtke's avatar
pschmidtke committed
347 348 349 350
      prms.A = prms.kij * rmin_pwr * rmin_pwr;
      prms.B = 2.0 * prms.kij * rmin_pwr;
      m_vdwTable[i][j] = prms;
      m_vdwTable[j][i] = prms;
351
      spdlog::debug("{},{},{},{},{},{},{}", triposType.Type2Str(RbtTriposAtomType::eType(i)), triposType.Type2Str(RbtTriposAtomType::eType(j)), prms.A, prms.B, prms.rmin, prms.kij, std::sqrt(prms.rmax_sq));
pschmidtke's avatar
pschmidtke committed
352 353
    }
  }
354
  // Now we can regenerate the close range params
pschmidtke's avatar
pschmidtke committed
355
  SetupCloseRange();
356
  // Set the overall range automatically from the max range for each atom type
357
  double range = *(std::max_element(m_maxRange.begin(), m_maxRange.end()));
pschmidtke's avatar
pschmidtke committed
358
  SetRange(range);
359
  spdlog::info("Overall max range for scoring function = {}", range);
pschmidtke's avatar
pschmidtke committed
360 361
}

362
// Regenerate the short-range params only (called more frequently)
pschmidtke's avatar
pschmidtke committed
363 364
void RbtVdwSF::SetupCloseRange() {
  RbtTriposAtomType triposType;
365 366
  spdlog::trace("RbtVdwSF::SetupCloseRange()");
  spdlog::debug("TYPE1,TYPE2,RCUT,ECUT,SLOPE,E0");
367 368
  double x = 1.0 + std::sqrt(1.0 + m_ecut);
  double p = (m_use_4_8) ? std::pow(x, 1.0 / 4.0) : std::pow(x, 1.0 / 6.0);
369
  double c = 1.0 / p;
370 371 372 373
  for (RbtVdwTableIter iter1 = m_vdwTable.begin(); iter1 != m_vdwTable.end();
       iter1++) {
    for (RbtVdwRowIter iter2 = (*iter1).begin(); iter2 != (*iter1).end();
         iter2++) {
374
      (*iter2).rcutoff_sq = std::pow((*iter2).rmin * c, 2);
pschmidtke's avatar
pschmidtke committed
375 376 377
      (*iter2).ecutoff = (*iter2).kij * m_ecut;
      (*iter2).e0 = (*iter2).ecutoff * m_e0;
      (*iter2).slope = ((*iter2).e0 - (*iter2).ecutoff) / (*iter2).rcutoff_sq;
378 379 380 381 382 383
      spdlog::debug("{},{},{},{},{},{}", triposType.Type2Str(RbtTriposAtomType::eType(iter1 - m_vdwTable.begin())),
              triposType.Type2Str(RbtTriposAtomType::eType(iter2 - (*iter1).begin())),
              std::sqrt((*iter2).rcutoff_sq),
              (*iter2).ecutoff,
              (*iter2).slope,
              (*iter2).e0);
pschmidtke's avatar
pschmidtke committed
384 385 386 387
    }
  }
}

388 389 390 391 392 393 394
// Index the flexible interactions between two atom lists.
// Foreach atom in the first list, compile a list of the atoms in the second
// list whose distance can vary to that atom intns container should have been
// initialised to a vector of empty atom lists prior to calling BuildIntraMap
void RbtVdwSF::BuildIntraMap(const RbtAtomRList &atomList1,
                             const RbtAtomRList &atomList2,
                             RbtAtomRListList &intns) const {
pschmidtke's avatar
pschmidtke committed
395 396
  Rbt::SelectAtom selectAtom(true);
  Rbt::isAtomSelected isSelected;
397 398
  bool bSingleList = atomList2.empty(); // If true, then index the flexible
                                        // interactions within atomList1
399 400 401 402 403
  for (RbtAtomRListConstIter iter = atomList1.begin(); iter != atomList1.end();
       iter++) {
    // 1. First select all atoms in the list
    std::for_each(atomList1.begin(), atomList1.end(), selectAtom);
    RbtAtom *pAtom = (*iter);
404
    unsigned int id = pAtom->GetAtomId() - 1;
pschmidtke's avatar
pschmidtke committed
405
    Assert<RbtAssert>(id >= 0 && id < intns.size());
406 407
    RbtModel *pModel = pAtom->GetModelPtr();
    // 2. Now deselect all atoms in the same model as the current atom
pschmidtke's avatar
pschmidtke committed
408
    pModel->SetAtomSelectionFlags(false);
409 410 411 412
    // 3. Now select all the flexible interactions to this atom within the
    // current model (Atoms in different models are not affected, therefore
    // remain selected from step 1) This way, we deal correctly with situations
    // where the atoms in atomList may belong to different models
pschmidtke's avatar
pschmidtke committed
413 414
    pModel->SelectFlexAtoms(pAtom);
    if (bSingleList) {
415 416 417 418 419
      std::copy_if(iter + 1, atomList1.end(), std::back_inserter(intns[id]),
                   isSelected);
    } else {
      std::copy_if(atomList2.begin(), atomList2.end(),
                   std::back_inserter(intns[id]), isSelected);
pschmidtke's avatar
pschmidtke committed
420
    }
421 422 423 424 425

    if (intns[id].empty()) {
        spdlog::debug("RbtVdwSF::BuildIntraMap: {}; N = 0", pAtom->GetFullAtomName());
    } else {
        spdlog::debug("RbtVdwSF::BuildIntraMap: {}; {} to {}; N = {}", pAtom->GetFullAtomName(), intns[id].front()->GetFullAtomName(), intns[id].back()->GetFullAtomName(), intns[id].size());
pschmidtke's avatar
pschmidtke committed
426 427 428 429
    }
  }
}

430 431 432 433 434 435
// Convenience method for above.
// Indexes the flexible interactions within a single atom list
void RbtVdwSF::BuildIntraMap(const RbtAtomRList &atomList,
                             RbtAtomRListList &intns) const {
  RbtAtomRList emptyList;
  BuildIntraMap(atomList, emptyList, intns);
pschmidtke's avatar
pschmidtke committed
436 437
}

438 439
void RbtVdwSF::Partition(const RbtAtomRList &atomList,
                         const RbtAtomRListList &intns,
440
                         RbtAtomRListList &prtIntns, double dist) const {
pschmidtke's avatar
pschmidtke committed
441

442 443 444
  // Clear the existing partitioned lists
  for (RbtAtomRListListIter iter = prtIntns.begin(); iter != prtIntns.end();
       iter++)
pschmidtke's avatar
pschmidtke committed
445
    (*iter).clear();
446
  spdlog::trace("RbtVdwSF: Partition (dist={})", dist);
pschmidtke's avatar
pschmidtke committed
447

448 449
  for (RbtAtomRListConstIter iter = atomList.begin(); iter != atomList.end();
       iter++) {
450
    int id = (*iter)->GetAtomId() - 1;
pschmidtke's avatar
pschmidtke committed
451
    if (dist > 0.0) {
452 453 454 455 456 457
      isD_lt bInRange(*iter, dist);
      // This copies all interactions which are within dist A of atom
      std::copy_if(intns[id].begin(), intns[id].end(),
                   std::back_inserter(prtIntns[id]), bInRange);
    } else {
      // Remove partition - simply copy from the master list of vdW interactions
pschmidtke's avatar
pschmidtke committed
458 459
      prtIntns[id] = intns[id];
    }
460
    spdlog::debug("RbtVdwSF: {}; $vdw={}: #prtn={}", (*iter)->GetFullAtomName(), intns[id].size(), prtIntns[id].size());
pschmidtke's avatar
pschmidtke committed
461 462
  }
}