dlpolytopologyreader.cc 12.6 KB
Newer Older
1
/*
2
 * Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

18
#include <boost/algorithm/string.hpp>
19 20 21 22 23 24
#include <boost/filesystem/convenience.hpp>
#include <boost/lexical_cast.hpp>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <votca/csg/topology.h>
25
#include <votca/tools/getline.h>
26
#include <votca/tools/tokenizer.h>
27 28 29 30

#ifndef HAVE_NO_CONFIG
#include <votca_config.h>
#endif
31 32 33

#include "dlpolytopologyreader.h"

34
using namespace votca::tools;
35 36
using namespace std;

37 38
namespace votca {
namespace csg {
39

40 41 42 43 44
string DLPOLYTopologyReader::_NextKeyline(ifstream &fs, const char *wspace)
// function to find and read the next line starting with a keyword/directive
// (skipping comments starting with "#" or ";") NOTE: the line is returned
// case-preserved, not to alter molecule/atom names (consider if --no-map is
// used)
45 46
{
  string line;
47
  size_t i_nws = 0;
48 49 50

  do {
    getline(fs, line);
51

52
    if (fs.eof()) {
53 54
      throw std::runtime_error("Error: unexpected end of dlpoly file '" +
                               _fname + "'");
55
    }
56

57
    i_nws = line.find_first_not_of(wspace);
58
  } while (line.substr(i_nws, 1) == "#" || line.substr(i_nws, 1) == ";");
59

60
  return line.substr(i_nws, line.size() - i_nws);
61 62
}

63
string DLPOLYTopologyReader::_NextKeyInt(ifstream &fs, const char *wspace,
Christoph Junghans's avatar
Christoph Junghans committed
64
                                         const string &word, long int &ival)
65 66 67 68
// function to read the next line containing only a given keyword and an integer
// value after it (only skipping comments!) NOTE: this function must only be
// called when the next directive line has to contain the given keyword and an
// integer value
69
{
70 71
  stringstream sl(_NextKeyline(fs, wspace));
  string line, sval;
72

73
  sl >> line;  // allow user not to bother about the case
74
  boost::to_upper(line);
75

76
  if (line.substr(0, word.size()) != word) {
77 78 79
    throw std::runtime_error("Error: unexpected line from dlpoly file '" +
                             _fname + "', expected '" + word + "' but got '" +
                             line + "'");
80
  }
81

82 83
  sl >> sval;

84
  size_t i_num = sval.find_first_of(
85
      "0123456789");  // assume integer number straight after the only keyword
86

87
  if (i_num > 0) {
88 89
    throw std::runtime_error("Error: missing integer number in directive '" +
                             line + "' in topology file '" + _fname + "'");
90
  }
91 92

  ival = boost::lexical_cast<int>(sval);
93 94 95 96

  return sl.str();
}

97 98 99 100 101
bool DLPOLYTopologyReader::_isKeyInt(const string &line, const char *wspace,
                                     const string &word, int &ival)
// function to check if the given (last read) directive line starts with a given
// keyword and has an integer value at the end NOTE: comments are allowed beyond
// the integer (anything after the first integer is ignored)
102
{
103 104 105
  // split directives consisting of a few words: the sought keyword must be the
  // first one!
  Tokenizer tok(line, wspace);
106 107 108 109 110
  vector<string> fields;
  tok.ToVector(fields);

  ival = 0;

111 112 113
  if (fields.size() < 2) {
    return false;
  }
114 115 116

  boost::to_upper(fields[0]);

117
  if (fields[0].substr(0, word.size()) != word) {
118 119 120
    throw std::runtime_error("Error: unexpected directive from dlpoly file '" +
                             _fname + "', expected keyword '" + word +
                             "' but got '" + fields[0] + "'");
121
  }
122

123
  size_t i_num = string::npos;
124

125
  unsigned int i = 1;
126 127
  do {  // find integer number in the field with the lowest index (closest to
        // the keyword)
128
    i_num = fields[i++].find_first_of("0123456789");
129
  } while (i_num > 0 && i < fields.size());
130

131 132 133
  if (i_num > 0) {
    return false;
  }
134

135
  ival = boost::lexical_cast<int>(fields[i - 1]);
136 137 138 139

  return true;
}

140 141
bool DLPOLYTopologyReader::ReadTopology(string file, Topology &top) {
  const char *WhiteSpace = " \t";
142

Christoph Junghans's avatar
Christoph Junghans committed
143 144
  long int matoms = 0;
  long int natoms = 0;
145

146 147
  std::ifstream fl;
  boost::filesystem::path filepath(file.c_str());
148

149
  string line;
150

151 152 153
  // TODO: fix residue naming / assignment - DL_POLY has no means to recognise
  // residues!
  Residue *res = top.CreateResidue("no");
154

155 156
  if (boost::filesystem::basename(filepath).size() == 0) {
    if (filepath.parent_path().string().size() == 0) {
157 158
      _fname = "FIELD";  // DL_POLY uses fixed file names in current/working
                         // directory
159
    } else {
160
      _fname = filepath.parent_path().string() + "/FIELD";
161
    }
162 163 164
  } else {
    _fname = file;
  }
165

166
  fl.open(_fname.c_str());
167

168 169 170
  if (!fl.is_open()) {
    throw std::runtime_error("Error on opening dlpoly file '" + _fname + "'");
  } else {
171

172 173
    line = _NextKeyline(fl, WhiteSpace);  // read title line and skip it
    line = _NextKeyline(fl, WhiteSpace);  // read next directive line
174
    boost::to_upper(line);
175

176 177
    if (line.substr(0, 4) == "UNIT") {      // skip 'unit' line
      line = _NextKeyline(fl, WhiteSpace);  // read next directive line
178 179
      boost::to_upper(line);
    }
180

181 182 183
    if (line.substr(0, 4) == "NEUT") {  // skip 'neutral groups' line (DL_POLY
                                        // Classic FIELD format)
      line = _NextKeyline(fl, WhiteSpace);  // look for next directive line
184 185
      boost::to_upper(line);
    }
186

187
    int nmol_types;
188

189
    if (!_isKeyInt(line, WhiteSpace, "MOLEC", nmol_types)) {
190 191
      throw std::runtime_error("Error: missing integer number in directive '" +
                               line + "' in topology file '" + _fname + "'");
192
    }
193

194
#ifdef DEBUG
195 196
    cout << "Read from dlpoly file '" << _fname << "' : '" << line << "' - "
         << nmol_types << endl;
197
#endif
198

199
    string mol_name;
200

201
    for (int nmol_type = 0; nmol_type < nmol_types; nmol_type++) {
202

203 204
      mol_name = _NextKeyline(fl, WhiteSpace);
      Molecule *mi = top.CreateMolecule(mol_name);
205

Christoph Junghans's avatar
Christoph Junghans committed
206
      long int nreplica = 1;
207
      line = _NextKeyInt(fl, WhiteSpace, "NUMMOL", nreplica);
208

209
#ifdef DEBUG
210 211
      cout << "Read from dlpoly file '" << _fname << "' : '" << mol_name
           << "' - '" << line << "' - " << nreplica << endl;
212
#endif
213

214
      line = _NextKeyInt(fl, WhiteSpace, "ATOMS", natoms);
215

216
#ifdef DEBUG
217 218
      cout << "Read from dlpoly file '" << _fname << "' : '" << line << "' - "
           << natoms << endl;
219
#endif
220

221
      // read molecule
Christoph Junghans's avatar
Christoph Junghans committed
222 223
      unique_ptr<long int[]> id_map =
          unique_ptr<long int[]>{new long int[natoms]};
224
      for (int i = 0; i < natoms;) {  // i is altered in repeater loop
225
        stringstream sl(_NextKeyline(fl, WhiteSpace));
226 227

#ifdef DEBUG
228 229
        cout << "Read atom specs from dlpoly topology : '" << sl.str() << "'"
             << endl;
230
#endif
231 232 233 234 235 236 237 238 239 240 241
        string beadtype;
        sl >> beadtype;
        if (!top.BeadTypeExist(beadtype)) {
          top.RegisterBeadType(beadtype);
        }
        double mass;
        sl >> mass;
        double charge;
        sl >> charge;

        line = " ";
242
        sl >> line;  // rest of the atom line
243 244 245 246

        Tokenizer tok(line, WhiteSpace);
        vector<string> fields;
        tok.ToVector(fields);
247

248
#ifdef DEBUG
249 250
        cout << "Rest atom specs from dlpoly topology : '" << line << "'"
             << endl;
251 252
#endif

253
        int repeater = 1;
254 255 256
        if (fields.size() > 1) {
          repeater = boost::lexical_cast<int>(fields[0]);
        }
257

258
        for (int j = 0; j < repeater; j++) {
259

260 261 262
          string beadname = beadtype + "#" + boost::lexical_cast<string>(i + 1);
          Bead *bead =
              top.CreateBead(1, beadname, beadtype, res->getId(), mass, charge);
263

264 265 266 267 268 269 270
          stringstream nm;
          nm << bead->getResnr() + 1 << ":"
             << top.getResidue(bead->getResnr())->getName() << ":"
             << bead->getName();
          mi->AddBead(bead, nm.str());
          id_map[i] = bead->getId();
          i++;
271
#ifdef DEBUG
272
          cout << "Atom identification in maps '" << nm.str() << "'" << endl;
273
#endif
274 275 276
        }
        matoms += repeater;
      }
277

278
      while (line != "FINISH") {
279

280 281
        stringstream nl(_NextKeyline(fl, WhiteSpace));
        nl >> line;
282 283

#ifdef DEBUG
284 285
        cout << "Read unit type# from dlpoly topology : '" << nl.str() << "'"
             << endl;
286
#endif
287

288 289 290 291 292 293 294
        boost::to_upper(line);
        line = line.substr(0, 6);
        if ((line == "BONDS") || (line == "ANGLES") || (line == "DIHEDR")) {
          string type = line;
          int count;
          nl >> count;
          for (int i = 0; i < count; i++) {
295

296
            stringstream sl(_NextKeyline(fl, WhiteSpace));
297
#ifdef DEBUG
298 299
            cout << "Read unit specs from dlpoly topology : '" << sl.str()
                 << "'" << endl;
300
#endif
301 302
            sl >> line;  // internal dlpoly bond/angle/dihedral function types
                         // are merely skipped (ignored)
303
            int ids[4];
304
            Interaction *ic = nullptr;
305 306 307 308
            sl >> ids[0];
            sl >> ids[1];
            if (type == "BONDS") {
              ic = new IBond(id_map[ids[0] - 1],
309
                             id_map[ids[1] - 1]);  // -1 due to fortran vs c
310 311 312
            } else if (type == "ANGLES") {
              sl >> ids[2];
              ic = new IAngle(id_map[ids[0] - 1], id_map[ids[1] - 1],
313
                              id_map[ids[2] - 1]);  // -1 due to fortran vs c
314 315 316 317 318 319
            } else if (type.substr(0, 6) == "DIHEDR") {
              type = "DIHEDRALS";
              sl >> ids[2];
              sl >> ids[3];
              ic = new IDihedral(id_map[ids[0] - 1], id_map[ids[1] - 1],
                                 id_map[ids[2] - 1],
320
                                 id_map[ids[3] - 1]);  // -1 due to fortran vs c
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335
            } else {
              throw std::runtime_error(
                  "Error: type should be BONDS, ANGLES or DIHEDRALS");
            }
            // could one use bond/angle/dihedral function types for 1:1 mapping?
            // (CG map overwrites ic->Group anyway)
            // ic->setGroup(line);
            ic->setGroup(type);
            ic->setIndex(i);
            ic->setMolecule(mi->getId());
            top.AddBondedInteraction(ic);
            mi->AddInteraction(ic);
          }
        }
      }
336

337
#ifdef DEBUG
338 339
      cout << "Read from dlpoly file '" << _fname << "' : '" << line
           << "' - done with '" << mol_name << "'" << endl;
340
#endif
341

342 343 344 345 346 347 348 349 350 351 352 353 354
      // replicate molecule
      for (int replica = 1; replica < nreplica; replica++) {
        Molecule *mi_replica = top.CreateMolecule(mol_name);
        for (int i = 0; i < mi->BeadCount(); i++) {
          Bead *bead = mi->getBead(i);
          string type = bead->getType();
          Bead *bead_replica =
              top.CreateBead(1, bead->getName(), type, res->getId(),
                             bead->getMass(), bead->getQ());
          mi_replica->AddBead(bead_replica, bead->getName());
        }
        matoms += mi->BeadCount();
        InteractionContainer ics = mi->Interactions();
355
        for (auto &ic : ics) {
356
          Interaction *ic_replica = nullptr;
Christoph Junghans's avatar
Christoph Junghans committed
357
          long int offset =
358
              mi_replica->getBead(0)->getId() - mi->getBead(0)->getId();
359 360 361 362 363 364 365 366
          if (ic->BeadCount() == 2) {
            ic_replica =
                new IBond(ic->getBeadId(0) + offset, ic->getBeadId(1) + offset);
          } else if (ic->BeadCount() == 3) {
            ic_replica =
                new IAngle(ic->getBeadId(0) + offset, ic->getBeadId(1) + offset,
                           ic->getBeadId(2) + offset);
          } else if (ic->BeadCount() == 4) {
367
            ic_replica = new IDihedral(
368 369
                ic->getBeadId(0) + offset, ic->getBeadId(1) + offset,
                ic->getBeadId(2) + offset, ic->getBeadId(3) + offset);
370 371 372
          } else {
            throw std::runtime_error("Error: BeadCount not equal 2, 3 or 4");
          }
373 374
          ic_replica->setGroup(ic->getGroup());
          ic_replica->setIndex(ic->getIndex());
375 376 377 378
          ic_replica->setMolecule(mi_replica->getId());
          top.AddBondedInteraction(ic_replica);
          mi_replica->AddInteraction(ic_replica);
        }
379 380
      }
    }
381 382
    top.RebuildExclusions();
  }
383

384
#ifdef DEBUG
385
  getline(fl, line);  // is "close" found?
386 387 388 389 390 391 392 393
  if (line == "close") {
    cout << "Read from dlpoly file '" << _fname << "' : '" << line
         << "' - done with topology" << endl;
  } else {
    cout << "Read from dlpoly file '" << _fname
         << "' : 'EOF' - done with topology (directive 'close' not read!)"
         << endl;
  }
394 395
#endif

396 397
  // we don't need the rest
  fl.close();
398

399
  return true;
400 401
}

402 403
}  // namespace csg
}  // namespace votca