NAMD
ReadAmberParm.C
Go to the documentation of this file.
1 #include "ReadAmberParm.h"
2 #include "InfoStream.h"
3 #include "common.h"
4 
5 #include <algorithm>
6 #include <cctype>
7 #include <cstring>
8 #include <fstream>
9 #include <sstream>
10 #include <stdexcept>
11 #include <utility>
12 #include <iomanip>
13 
14 namespace AmberParm7Reader {
15 
16 void parse_fortran_format(const std::string &str,
17  FortranFormatSpecifier &specifier) {
18  size_t left_quotation_mark = str.find("'");
19  if (left_quotation_mark == std::string::npos) {
20  std::string num_fields, type, width, dot, precision;
21  for (size_t i = 0; i < str.size(); ++i) {
22  // check if this is a digit at first.
23  if (std::isdigit(str[i])) {
24  // then check if the digit is before or after the type.
25  if (type.size() > 0) {
26  // if type is already read,
27  // then append this character to "width" or "precision".
28  if (dot.size() > 0) {
29  // the digit after "." denotes the precision.
30  precision += str[i];
31  } else {
32  // the digit is after the type field but before ".".
33  width += str[i];
34  }
35  } else {
36  // the digit before the type represents the number of fields.
37  num_fields += str[i];
38  }
39  } else if (std::isalpha(str[i])) {
40  // if it is [A-Za-z] then it denotes the type.
41  type += str[i];
42  } else if (str[i] == '.') {
43  dot += str[i];
44  } else {
45  // skip other characters
46  continue;
47  }
48  }
49  // populate the FortranFormatSpecifier struct
50  if (num_fields.size() > 0) {
51  specifier.NumOfFields = std::stoul(num_fields.c_str());
52  } else {
53  // ("i10", "E6.2")
54  specifier.NumOfFields = 1;
55  }
56  // must have "type"
57  specifier.Type = type.front();
58  // must have "width"
59  specifier.Width = std::stoul(width.c_str());
60  if (precision.size() > 0) {
61  specifier.Precision = std::stoul(precision.c_str());
62  } else {
63  // integers have no precision fields
64  specifier.Precision = -1;
65  }
66  } else {
67  size_t right_quotation_mark = str.find("'", left_quotation_mark + 1);
68  specifier.IsPadding = true;
69  specifier.Padding =
70  str.substr(left_quotation_mark + 1,
71  right_quotation_mark - (left_quotation_mark + 1));
72  specifier.NumOfFields = 1;
73  }
74 }
75 
76 // helper function to split a string by delimeters
77 void split_string(const std::string& data,
78  const std::string& delim,
79  std::vector<std::string>& dest) {
80  size_t index = 0, new_index = 0;
81  std::string tmpstr;
82  while (index != data.length()) {
83  new_index = data.find(delim, index);
84  if (new_index != std::string::npos) tmpstr = data.substr(index, new_index - index);
85  else tmpstr = data.substr(index, data.length());
86  if (!tmpstr.empty()) {
87  dest.push_back(tmpstr);
88  }
89  if (new_index == std::string::npos) break;
90  index = new_index + 1;
91  }
92 }
93 
94 void parse_fortran_format(const std::string &str,
95  vector<FortranFormatSpecifier> &specifiers) {
96  // split the string by comma
97  vector<string> sub_format_strings;
98  split_string(str, ",", sub_format_strings);
99  specifiers.resize(sub_format_strings.size());
100  for (size_t i = 0; i < sub_format_strings.size(); ++i) {
101  parse_fortran_format(sub_format_strings[i], specifiers[i]);
102  }
103 }
104 
106  const string &source, const vector<FortranFormatSpecifier> &specifiers,
107  vector<FortranData> &destination) {
108  size_t start_pos = 0;
109  for (const auto &spec : specifiers) {
110  if (spec.IsPadding == true) {
111  // I should check padding here but I don't know how to handle the it
112  // Anyway, AMBER parm7 doesn't use padding
113  } else {
114  for (int i = 0; i < spec.NumOfFields; ++i) {
115  if (start_pos < source.size()) {
116  const string tmp = source.substr(start_pos, spec.Width);
117  try {
118  switch (toupper(spec.Type)) {
119  case 'I':
120  destination.push_back(std::stoi(tmp));
121  break;
122  case 'A':
123  destination.push_back(tmp);
124  break;
125  case 'D':
126  case 'E':
127  case 'G':
128  case 'F':
129  // fall through for all floating point specifiers
130 #ifdef DOUBLE
131  destination.push_back(std::stod(tmp));
132 #else
133  destination.push_back(std::stof(tmp));
134 #endif
135  break;
136  }
137  start_pos += spec.Width;
138  } catch (...) {
139  std::cerr << "Error on parsing field \"" << tmp
140  << "\" with specifier " << spec.Type << std::endl;
141  throw;
142  }
143  }
144  }
145  }
146  }
147 }
148 
149 bool read_amber_parm_stage1(const char *filename, AmberTopparMap &toppar_map) {
150  using std::ifstream;
151  using std::istringstream;
152  using std::string;
153  bool success = true;
154  string error_msg;
155  ifstream ifs_parm(filename);
156  if (!ifs_parm.is_open()) return false;
157  string line;
158  string current_flag;
159  bool match_flag = false;
160  bool match_format = false;
161  // vector<FortranData> section_data;
162  vector<FortranFormatSpecifier> specifiers;
163  tuple<bool, vector<FortranData>> *p_map = nullptr;
164  istringstream iss;
165  size_t line_no = 0;
166  while (std::getline(ifs_parm, line)) {
167  // Remove any \r for compatibility with windows CR LF files
168  // Code from: https://stackoverflow.com/a/72365365/8807519
169  line.erase(std::remove(line.begin(), line.end(), '\r'), line.end());
170  /* AMBER format parser requirement 3:
171  * Parsers should not expect or require %COMMENT lines to exist,
172  * but should properly parse the file if any number of %COMMENT
173  * lines appear as indicated above.
174  */
175  if (line[0] == '%') {
176  if (line.find("%COMMENT") == 0) {
177  continue;
178  } else if (line.find("%VERSION") == 0) {
179  // parse the version line, read version_stamp and build time
180  char tmp_version_stamp[16];
181  char tmp_date[16];
182  char tmp_time[16];
183  if (sscanf(line.c_str(), "%%VERSION VERSION_STAMP = %9s DATE = %8s %8s",
184  tmp_version_stamp, tmp_date, tmp_time) == 3) {
185  continue;
186  } else {
187  return false;
188  }
189  } else if (line.find("%FLAG") == 0) {
190  // parse flags
191  if (match_flag && match_format) {
192  // insert the previous section into map
193  current_flag.clear();
194  specifiers.clear();
195  match_flag = false;
196  match_format = false;
197  }
198  istringstream iss(line.substr(strlen("%FLAG")));
199  iss >> current_flag;
200  p_map = &(toppar_map[current_flag]);
201  std::get<0>(*p_map) = false;
202  match_flag = true;
203  continue;
204  } else if (line.find("%FORMAT") == 0) {
205  // parse format
206  if (match_flag == false) {
207  success = false;
208  break;
209  }
210  const size_t format_start = line.find("(", strlen("%FORMAT"));
211  const size_t format_end = line.find(")", format_start + 1);
213  line.substr(format_start + 1, format_end - (format_start + 1)),
214  specifiers);
215  match_format = true;
216  continue;
217  }
218  }
219  // parse data fields
220  if (match_flag && match_format) {
221  try {
222  // read the according section
223  split_string_by_specifiers(line, specifiers, std::get<1>(*p_map));
224  } catch (...) {
225  const std::string error_msg =
226  "Error on parsing AMBER parameter file at section " +
227  current_flag + ", line " + std::to_string(line_no) + ":\"" + line + "\"\n";
228  NAMD_err(error_msg.c_str());
229  }
230  }
231  ++line_no;
232  }
233  return success;
234 }
235 
237  Ambertoppar &toppar_data) {
238  using std::get;
239  using std::ref;
240  bool success = true;
241  // TITLE
242  {
243  // check TITLE to see if this is a AMBER force field
244  auto current_section = toppar_map.find("TITLE");
245  if (current_section == toppar_map.end()) {
246  // check if this is a CHARMM force field in AMBER format
247  current_section = toppar_map.find("CTITLE");
248  if (current_section == toppar_map.end()) {
249  iout << iERROR << "TITLE/CTITLE section for AMBER/CHARMM force field is not "
250  "found!\n" << endi;
251  toppar_data.IsCharmmFF = false;
252  success = false;
253  iout << iERROR << "Failed to read the title section!\n" << endi;
254  } else {
255  iout << iWARN << "You are using CHARMM force field in AMBER format, which is not supported in NAMD!\n" << endi;
256  toppar_data.IsCharmmFF = true;
257  }
258  } else {
259  toppar_data.IsCharmmFF = false;
260  }
261  toppar_data.ititl = get<1>(current_section->second).size() > 0 ? get<1>(current_section->second)[0].String : string();
262  get<0>(current_section->second) = true;
263  }
264  // POINTERS
265  if (success) {
266  const string section_name{"POINTERS"};
267  const auto &current_section = toppar_map.find(section_name);
268  if (current_section != toppar_map.end()) {
269  success = parse_pointer(get<1>(current_section->second), toppar_data);
270  get<0>(current_section->second) = true;
271  } else {
272  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
273  }
274  }
275  // ATOM_NAME
276  if (success) {
277  const string section_name{"ATOM_NAME"};
278  const auto &current_section = toppar_map.find(section_name);
279  if (current_section != toppar_map.end()) {
280  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
281  toppar_data.AtomNames, section_name);
282  get<0>(current_section->second) = true;
283  } else {
284  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
285  }
286  }
287  // CHARGE
288  if (success) {
289  const string section_name{"CHARGE"};
290  const auto &current_section = toppar_map.find(section_name);
291  if (current_section != toppar_map.end()) {
292  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
293  toppar_data.Charges, section_name);
294  get<0>(current_section->second) = true;
295  } else {
296  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
297  }
298  }
299  // ATOMIC_NUMBER
300  if (success) {
301  const string section_name{"ATOMIC_NUMBER"};
302  const auto &current_section = toppar_map.find(section_name);
303  if (current_section != toppar_map.end()) {
304  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
305  toppar_data.AtomNumbers, section_name);
306  get<0>(current_section->second) = true;
307  } else {
308  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
309  }
310  }
311  // MASS
312  if (success) {
313  const string section_name{"MASS"};
314  const auto &current_section = toppar_map.find(section_name);
315  if (current_section != toppar_map.end()) {
316  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
317  toppar_data.Masses, section_name);
318  get<0>(current_section->second) = true;
319  } else {
320  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
321  }
322  }
323  // ATOM_TYPE_INDEX
324  if (success) {
325  const string section_name{"ATOM_TYPE_INDEX"};
326  const auto &current_section = toppar_map.find(section_name);
327  if (current_section != toppar_map.end()) {
328  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
329  toppar_data.Iac, section_name);
330  get<0>(current_section->second) = true;
331  } else {
332  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
333  }
334  }
335  // NUMBER_EXCLUDED_ATOMS
336  if (success) {
337  const string section_name{"NUMBER_EXCLUDED_ATOMS"};
338  const auto &current_section = toppar_map.find(section_name);
339  if (current_section != toppar_map.end()) {
340  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
341  toppar_data.Iblo, section_name);
342  get<0>(current_section->second) = true;
343  } else {
344  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
345  }
346  }
347  // NONBONDED_PARM_INDEX
348  if (success) {
349  const int Ntype2d = toppar_data.Ntypes * toppar_data.Ntypes;
350  const string section_name{"NONBONDED_PARM_INDEX"};
351  const auto &current_section = toppar_map.find(section_name);
352  if (current_section != toppar_map.end()) {
353  success = parse_section(get<1>(current_section->second), Ntype2d, toppar_data.Cno,
354  section_name);
355  get<0>(current_section->second) = true;
356  } else {
357  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
358  }
359  }
360  // RESIDUE_LABEL
361  if (success) {
362  const string section_name{"RESIDUE_LABEL"};
363  const auto &current_section = toppar_map.find(section_name);
364  if (current_section != toppar_map.end()) {
365  success = parse_section(get<1>(current_section->second), toppar_data.Nres,
366  toppar_data.ResNames, section_name);
367  get<0>(current_section->second) = true;
368  } else {
369  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
370  }
371  }
372  // RESIDUE_POINTER
373  if (success) {
374  const string section_name{"RESIDUE_POINTER"};
375  const auto &current_section = toppar_map.find(section_name);
376  if (current_section != toppar_map.end()) {
377  success = parse_section(get<1>(current_section->second), toppar_data.Nres,
378  toppar_data.Ipres, section_name);
379  get<0>(current_section->second) = true;
380  } else {
381  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
382  }
383  }
384  // BOND_FORCE_CONSTANT
385  if (success) {
386  const string section_name{"BOND_FORCE_CONSTANT"};
387  const auto &current_section = toppar_map.find(section_name);
388  if (current_section != toppar_map.end()) {
389  success = parse_section(get<1>(current_section->second), toppar_data.Numbnd,
390  toppar_data.Rk, section_name);
391  get<0>(current_section->second) = true;
392  } else {
393  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
394  }
395  }
396  // BOND_EQUIL_VALUE
397  if (success) {
398  const string section_name{"BOND_EQUIL_VALUE"};
399  const auto &current_section = toppar_map.find(section_name);
400  if (current_section != toppar_map.end()) {
401  success = parse_section(get<1>(current_section->second), toppar_data.Numbnd,
402  toppar_data.Req, section_name);
403  get<0>(current_section->second) = true;
404  } else {
405  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
406  }
407  }
408  // ANGLE_FORCE_CONSTANT
409  if (success) {
410  const string section_name{"ANGLE_FORCE_CONSTANT"};
411  const auto &current_section = toppar_map.find(section_name);
412  if (current_section != toppar_map.end()) {
413  success = parse_section(get<1>(current_section->second), toppar_data.Numang,
414  toppar_data.Tk, section_name);
415  get<0>(current_section->second) = true;
416  } else {
417  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
418  }
419  }
420  // ANGLE_EQUIL_VALUE
421  if (success) {
422  const string section_name{"ANGLE_EQUIL_VALUE"};
423  const auto &current_section = toppar_map.find(section_name);
424  if (current_section != toppar_map.end()) {
425  success = parse_section(get<1>(current_section->second), toppar_data.Numang,
426  toppar_data.Teq, section_name);
427  get<0>(current_section->second) = true;
428  } else {
429  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
430  }
431  }
432  // DIHEDRAL_FORCE_CONSTANT
433  if (success) {
434  const string section_name{"DIHEDRAL_FORCE_CONSTANT"};
435  const auto &current_section = toppar_map.find(section_name);
436  if (current_section != toppar_map.end()) {
437  success = parse_section(get<1>(current_section->second), toppar_data.Nptra,
438  toppar_data.Pk, section_name);
439  get<0>(current_section->second) = true;
440  } else {
441  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
442  }
443  }
444  // DIHEDRAL_PERIODICITY
445  if (success) {
446  const string section_name{"DIHEDRAL_PERIODICITY"};
447  const auto &current_section = toppar_map.find(section_name);
448  if (current_section != toppar_map.end()) {
449  success = parse_section(get<1>(current_section->second), toppar_data.Nptra,
450  toppar_data.Pn, section_name);
451  get<0>(current_section->second) = true;
452  } else {
453  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
454  }
455  }
456  // DIHEDRAL_PHASE
457  if (success) {
458  const string section_name{"DIHEDRAL_PHASE"};
459  const auto &current_section = toppar_map.find(section_name);
460  if (current_section != toppar_map.end()) {
461  success = parse_section(get<1>(current_section->second), toppar_data.Nptra,
462  toppar_data.Phase, section_name);
463  get<0>(current_section->second) = true;
464  } else {
465  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
466  }
467  }
468  // SCEE_SCALE_FACTOR
469  if (success) {
470  const string section_name{"SCEE_SCALE_FACTOR"};
471  const auto &current_section = toppar_map.find(section_name);
472  if (current_section != toppar_map.end()) {
473  success = parse_section(get<1>(current_section->second), toppar_data.Nptra,
474  toppar_data.SceeScaleFactor, section_name);
475  get<0>(current_section->second) = true;
476  } else {
477  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
478  }
479  }
480  // SCNB_SCALE_FACTOR
481  if (success) {
482  const string section_name{"SCNB_SCALE_FACTOR"};
483  const auto &current_section = toppar_map.find(section_name);
484  if (current_section != toppar_map.end()) {
485  success = parse_section(get<1>(current_section->second), toppar_data.Nptra,
486  toppar_data.ScnbScaleFactor, section_name);
487  get<0>(current_section->second) = true;
488  } else {
489  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
490  }
491  }
492  // SOLTY (UNUSED)
493  if (success) {
494  const string section_name{"SOLTY"};
495  const auto &current_section = toppar_map.find(section_name);
496  if (current_section != toppar_map.end()) {
497  // success = parse_section(get<1>(current_section->second),
498  // toppar_data.Natyp, toppar_data.Solty,
499  // section_name);
500  get<0>(current_section->second) = true;
501  } else {
502  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
503  }
504  }
505  // LENNARD_JONES_ACOEF
506  if (success) {
507  const int Nttyp = toppar_data.Ntypes * (toppar_data.Ntypes + 1) / 2;
508  const string section_name{"LENNARD_JONES_ACOEF"};
509  const auto &current_section = toppar_map.find(section_name);
510  if (current_section != toppar_map.end()) {
511  success = parse_section(get<1>(current_section->second), Nttyp, toppar_data.Cn1,
512  section_name);
513  get<0>(current_section->second) = true;
514  } else {
515  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
516  }
517  }
518  // LENNARD_JONES_BCOEF
519  if (success) {
520  const int Nttyp = toppar_data.Ntypes * (toppar_data.Ntypes + 1) / 2;
521  const string section_name{"LENNARD_JONES_BCOEF"};
522  const auto &current_section = toppar_map.find(section_name);
523  if (current_section != toppar_map.end()) {
524  success = parse_section(get<1>(current_section->second), Nttyp, toppar_data.Cn2,
525  section_name);
526  get<0>(current_section->second) = true;
527  } else {
528  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
529  }
530  }
531  // BONDS_INC_HYDROGEN
532  if (success) {
533  const string section_name{"BONDS_INC_HYDROGEN"};
534  const auto &current_section = toppar_map.find(section_name);
535  if (current_section != toppar_map.end()) {
536  success = parse_section(get<1>(current_section->second), toppar_data.Nbonh,
537  {ref(toppar_data.BondHAt1), ref(toppar_data.BondHAt2),
538  ref(toppar_data.BondHNum)},
539  section_name);
540  get<0>(current_section->second) = true;
541  } else {
542  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
543  }
544  }
545  // BONDS_WITHOUT_HYDROGEN
546  if (success) {
547  const string section_name{"BONDS_WITHOUT_HYDROGEN"};
548  const auto &current_section = toppar_map.find(section_name);
549  if (current_section != toppar_map.end()) {
550  success = parse_section(get<1>(current_section->second), toppar_data.Nbona,
551  {ref(toppar_data.BondAt1), ref(toppar_data.BondAt2),
552  ref(toppar_data.BondNum)},
553  section_name);
554  get<0>(current_section->second) = true;
555  } else {
556  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
557  }
558  }
559  // ANGLES_INC_HYDROGEN
560  if (success) {
561  const string section_name{"ANGLES_INC_HYDROGEN"};
562  const auto &current_section = toppar_map.find(section_name);
563  if (current_section != toppar_map.end()) {
564  success = parse_section(get<1>(current_section->second), toppar_data.Ntheth,
565  {ref(toppar_data.AngleHAt1), ref(toppar_data.AngleHAt2),
566  ref(toppar_data.AngleHAt3), ref(toppar_data.AngleHNum)},
567  section_name);
568  get<0>(current_section->second) = true;
569  } else {
570  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
571  }
572  }
573  // ANGLES_WITHOUT_HYDROGEN
574  if (success) {
575  const string section_name{"ANGLES_WITHOUT_HYDROGEN"};
576  const auto &current_section = toppar_map.find(section_name);
577  if (current_section != toppar_map.end()) {
578  success = parse_section(get<1>(current_section->second), toppar_data.Ntheta,
579  {ref(toppar_data.AngleAt1), ref(toppar_data.AngleAt2),
580  ref(toppar_data.AngleAt3), ref(toppar_data.AngleNum)},
581  section_name);
582  get<0>(current_section->second) = true;
583  } else {
584  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
585  }
586  }
587  // DIHEDRALS_INC_HYDROGEN
588  if (success) {
589  const string section_name{"DIHEDRALS_INC_HYDROGEN"};
590  const auto &current_section = toppar_map.find(section_name);
591  if (current_section != toppar_map.end()) {
592  success = parse_section(get<1>(current_section->second), toppar_data.Nphih,
593  {ref(toppar_data.DihHAt1), ref(toppar_data.DihHAt2),
594  ref(toppar_data.DihHAt3), ref(toppar_data.DihHAt4),
595  ref(toppar_data.DihHNum)},
596  section_name);
597  get<0>(current_section->second) = true;
598  } else {
599  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
600  }
601  }
602  // DIHEDRALS_WITHOUT_HYDROGEN
603  if (success) {
604  const string section_name{"DIHEDRALS_WITHOUT_HYDROGEN"};
605  const auto &current_section = toppar_map.find(section_name);
606  if (current_section != toppar_map.end()) {
607  success = parse_section(get<1>(current_section->second), toppar_data.Nphia,
608  {ref(toppar_data.DihAt1), ref(toppar_data.DihAt2),
609  ref(toppar_data.DihAt3), ref(toppar_data.DihAt4),
610  ref(toppar_data.DihNum)},
611  section_name);
612  get<0>(current_section->second) = true;
613  } else {
614  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
615  }
616  }
617  // EXCLUDED_ATOMS_LIST
618  if (success) {
619  const string section_name{"EXCLUDED_ATOMS_LIST"};
620  const auto &current_section = toppar_map.find(section_name);
621  if (current_section != toppar_map.end()) {
622  success = parse_section(get<1>(current_section->second), toppar_data.Nnb,
623  toppar_data.ExclAt, section_name);
624  get<0>(current_section->second) = true;
625  } else {
626  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
627  }
628  }
629  // HBOND_ACOEF
630  if (success) {
631  const string section_name{"HBOND_ACOEF"};
632  const auto &current_section = toppar_map.find(section_name);
633  if (current_section != toppar_map.end()) {
634  success = parse_section(get<1>(current_section->second), toppar_data.Nphb,
635  toppar_data.HB12, section_name);
636  get<0>(current_section->second) = true;
637  } else {
638  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
639  }
640  }
641  // HBOND_BCOEF
642  if (success) {
643  const string section_name{"HBOND_BCOEF"};
644  const auto &current_section = toppar_map.find(section_name);
645  if (current_section != toppar_map.end()) {
646  success = parse_section(get<1>(current_section->second), toppar_data.Nphb,
647  toppar_data.HB10, section_name);
648  get<0>(current_section->second) = true;
649  } else {
650  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
651  }
652  }
653  // HBCUT (UNUSED)
654  if (success) {
655  const string section_name{"HBCUT"};
656  const auto &current_section = toppar_map.find(section_name);
657  if (current_section != toppar_map.end()) {
658  // skip HBCUT
659  get<0>(current_section->second) = true;
660  } else {
661  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
662  }
663  }
664  // AMBER_ATOM_TYPE
665  if (success) {
666  const string section_name{"AMBER_ATOM_TYPE"};
667  const auto &current_section = toppar_map.find(section_name);
668  if (current_section != toppar_map.end()) {
669  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
670  toppar_data.AtomSym, section_name);
671  get<0>(current_section->second) = true;
672  } else {
673  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
674  }
675  }
676  // TREE_CHAIN_CLASSIFICATION
677  if (success) {
678  const string section_name{"TREE_CHAIN_CLASSIFICATION"};
679  const auto &current_section = toppar_map.find(section_name);
680  if (current_section != toppar_map.end()) {
681  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
682  toppar_data.AtomTree, section_name);
683  get<0>(current_section->second) = true;
684  } else {
685  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
686  }
687  }
688  // JOIN_ARRAY (UNUSED)
689  if (success) {
690  const string section_name{"JOIN_ARRAY"};
691  const auto &current_section = toppar_map.find(section_name);
692  if (current_section != toppar_map.end()) {
693  // success = parse_section(get<1>(current_section->second),
694  // toppar_data.Natom, toppar_data.TreeJoin,
695  // section_name);
696  get<0>(current_section->second) = true;
697  } else {
698  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
699  }
700  }
701  // IROTAT (UNUSED)
702  if (success) {
703  const string section_name{"IROTAT"};
704  const auto &current_section = toppar_map.find(section_name);
705  if (current_section != toppar_map.end()) {
706 // success = parse_section(get<1>(current_section->second),
707 // toppar_data.Natom, toppar_data.IrotatArray,
708 // section_name);
709  get<0>(current_section->second) = true;
710  } else {
711  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
712  }
713  }
714  // parse periodic boundary conditions only if IFBOX > 0
715  if (toppar_data.Ifbox > 0) {
716  // SOLVENT_POINTERS
717  if (success) {
718  const string section_name{"SOLVENT_POINTERS"};
719  const auto &current_section = toppar_map.find(section_name);
720  if (current_section != toppar_map.end()) {
721  if (get<1>(current_section->second).size() == 3) {
722  toppar_data.Iptres = get<1>(current_section->second)[0].Int;
723  toppar_data.Nspm = get<1>(current_section->second)[1].Int;
724  toppar_data.Nspsol = get<1>(current_section->second)[2].Int;
725  toppar_data.AtomsPerMolecule.reserve(toppar_data.Nspm);
726  toppar_data.BoxDimensions.reserve(4);
727  toppar_data.CapInfo.reserve(1);
728  toppar_data.CapInfo2.reserve(4);
729  get<0>(current_section->second) = true;
730  } else {
731  iout << iERROR << "Expect " << 3 << " fields but only "
732  << get<1>(current_section->second).size() << " are in "
733  << section_name.c_str() << "\n";
734  success = false;
735  }
736  } else {
737  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
738  }
739  }
740  // ATOMS_PER_MOLECULE
741  if (success) {
742  const string section_name{"ATOMS_PER_MOLECULE"};
743  const auto &current_section = toppar_map.find(section_name);
744  if (current_section != toppar_map.end()) {
745  success = parse_section(get<1>(current_section->second), toppar_data.Nspm,
746  toppar_data.AtomsPerMolecule, section_name);
747  get<0>(current_section->second) = true;
748  } else {
749  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
750  }
751  }
752  // BOX_DIMENSIONS
753  if (success) {
754  const string section_name{"BOX_DIMENSIONS"};
755  const auto &current_section = toppar_map.find(section_name);
756  if (current_section != toppar_map.end()) {
757  success = parse_section(get<1>(current_section->second), 4,
758  toppar_data.BoxDimensions, section_name);
759  get<0>(current_section->second) = true;
760  } else {
761  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
762  }
763  }
764  }
765  if (toppar_data.Ifcap > 0) {
766  // CAP_INFO
767  if (success) {
768  const string section_name{"CAP_INFO"};
769  const auto &current_section = toppar_map.find(section_name);
770  if (current_section != toppar_map.end()) {
771  success = parse_section(get<1>(current_section->second), 1, toppar_data.CapInfo,
772  section_name);
773  get<0>(current_section->second) = true;
774  } else {
775  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
776  }
777  }
778  // CAP_INFO2
779  if (success) {
780  const string section_name{"CAP_INFO2"};
781  const auto &current_section = toppar_map.find(section_name);
782  if (current_section != toppar_map.end()) {
783  success = parse_section(get<1>(current_section->second), 4, toppar_data.CapInfo2,
784  section_name);
785  get<0>(current_section->second) = true;
786  } else {
787  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
788  }
789  }
790  }
791  // RADIUS_SET
792  if (success) {
793  const string section_name{"RADIUS_SET"};
794  const auto &current_section = toppar_map.find(section_name);
795  if (current_section != toppar_map.end()) {
796  success = parse_section(get<1>(current_section->second), 1, toppar_data.RadiusSet,
797  section_name);
798  get<0>(current_section->second) = true;
799  } else {
800  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
801  }
802  }
803  // RADII
804  if (success) {
805  const string section_name{"RADII"};
806  const auto &current_section = toppar_map.find(section_name);
807  if (current_section != toppar_map.end()) {
808  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
809  toppar_data.Radii, section_name);
810  get<0>(current_section->second) = true;
811  } else {
812  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
813  }
814  }
815  // SCREEN
816  if (success) {
817  const string section_name{"SCREEN"};
818  const auto &current_section = toppar_map.find(section_name);
819  if (current_section != toppar_map.end()) {
820  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
821  toppar_data.Screen, section_name);
822  get<0>(current_section->second) = true;
823  } else {
824  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
825  }
826  }
827  // parse IPOL to determine the polarizability
828  if (success) {
829  const string section_name{"IPOL"};
830  const auto &current_section = toppar_map.find(section_name);
831  if (current_section != toppar_map.end()) {
832  toppar_data.Ipol = get<1>(current_section->second)[0].Int;
833  toppar_data.Polarizability.reserve(toppar_data.Natom);
834  get<0>(current_section->second) = true;
835  } else {
836  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
837  }
838  }
839  if (toppar_data.Ipol > 0) {
840  // POLARIZABILITY
841  if (success) {
842  const string section_name{"POLARIZABILITY"};
843  const auto &current_section = toppar_map.find(section_name);
844  if (current_section != toppar_map.end()) {
845  success = parse_section(get<1>(current_section->second), toppar_data.Natom,
846  toppar_data.Polarizability, section_name);
847  get<0>(current_section->second) = true;
848  } else {
849  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
850  }
851  }
852  }
853  if (toppar_data.IsCharmmFF) {
854  // I don't think supporting CHARMBER is a good idea since we have supported,
855  // CHARMM force field in PSF format, and the user should use PSF instead of
856  // AMBER format.
857  // CHARMM_UREY_BRADLEY_COUNT
858  if (success) {
859  const string section_name{"CHARMM_UREY_BRADLEY_COUNT"};
860  const auto &current_section = toppar_map.find(section_name);
861  if (current_section != toppar_map.end()) {
862  if (get<1>(current_section->second).size() == 2) {
863  toppar_data.Nub = get<1>(current_section->second)[0].Int;
864  toppar_data.Nubtypes = get<1>(current_section->second)[1].Int;
865  toppar_data.UreyBradleyAt1.reserve(toppar_data.Nub);
866  toppar_data.UreyBradleyAt2.reserve(toppar_data.Nub);
867  toppar_data.UreyBradleyNum.reserve(toppar_data.Nub);
868  toppar_data.UreyBradleyForceConstants.reserve(toppar_data.Nubtypes);
869  toppar_data.UreyBradleyEquilValues.reserve(toppar_data.Nubtypes);
870  get<0>(current_section->second) = true;
871  } else {
872  iout << iERROR << "Expect " << 2 << " fields but only "
873  << get<1>(current_section->second).size() << " are in "
874  << section_name.c_str() << "\n";
875  success = false;
876  }
877  } else {
878  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
879  }
880  }
881  // CHARMM_UREY_BRADLEY
882  if (success) {
883  const string section_name{"CHARMM_UREY_BRADLEY"};
884  const auto &current_section = toppar_map.find(section_name);
885  if (current_section != toppar_map.end()) {
886  success = parse_section(get<1>(current_section->second), toppar_data.Nub,
887  {ref(toppar_data.UreyBradleyAt1),
888  ref(toppar_data.UreyBradleyAt2),
889  ref(toppar_data.UreyBradleyNum)},
890  section_name);
891  get<0>(current_section->second) = true;
892  } else {
893  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
894  }
895  }
896  // CHARMM_UREY_BRADLEY_FORCE_CONSTANT
897  if (success) {
898  const string section_name{"CHARMM_UREY_BRADLEY_FORCE_CONSTANT"};
899  const auto &current_section = toppar_map.find(section_name);
900  if (current_section != toppar_map.end()) {
901  success = parse_section(get<1>(current_section->second), toppar_data.Nubtypes,
902  toppar_data.UreyBradleyForceConstants, section_name);
903  get<0>(current_section->second) = true;
904  } else {
905  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
906  }
907  }
908  // CHARMM_UREY_BRADLEY_EQUIL_VALUE
909  if (success) {
910  const string section_name{"CHARMM_UREY_BRADLEY_EQUIL_VALUE"};
911  const auto &current_section = toppar_map.find(section_name);
912  if (current_section != toppar_map.end()) {
913  success = parse_section(get<1>(current_section->second), toppar_data.Nubtypes,
914  toppar_data.UreyBradleyEquilValues, section_name);
915  get<0>(current_section->second) = true;
916  } else {
917  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
918  }
919  }
920  // CHARMM_NUM_IMPROPERS
921  if (success) {
922  const string section_name{"CHARMM_NUM_IMPROPERS"};
923  const auto &current_section = toppar_map.find(section_name);
924  if (current_section != toppar_map.end()) {
925  if (get<1>(current_section->second).size() == 1) {
926  toppar_data.Nimphi = get<1>(current_section->second)[0].Int;
927  toppar_data.ImproperAt1.reserve(toppar_data.Nimphi);
928  toppar_data.ImproperAt2.reserve(toppar_data.Nimphi);
929  toppar_data.ImproperAt3.reserve(toppar_data.Nimphi);
930  toppar_data.ImproperAt4.reserve(toppar_data.Nimphi);
931  toppar_data.ImproperNum.reserve(toppar_data.Nimphi);
932  get<0>(current_section->second) = true;
933  } else {
934  iout << iERROR << "Expect " << 1 << " fields but only "
935  << get<1>(current_section->second).size() << " are in "
936  << section_name.c_str() << "\n";
937  success = false;
938  }
939  } else {
940  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
941  }
942  }
943  // CHARMM_IMPROPERS
944  if (success) {
945  const string section_name{"CHARMM_IMPROPERS"};
946  const auto &current_section = toppar_map.find(section_name);
947  if (current_section != toppar_map.end()) {
948  success = parse_section(
949  get<1>(current_section->second), toppar_data.Nimphi,
950  {ref(toppar_data.ImproperAt1), ref(toppar_data.ImproperAt2),
951  ref(toppar_data.ImproperAt3), ref(toppar_data.ImproperAt4),
952  ref(toppar_data.ImproperNum)},
953  section_name);
954  get<0>(current_section->second) = true;
955  } else {
956  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
957  }
958  }
959  // CHARMM_NUM_IMPR_TYPES
960  // Not "CHARMM_NUM_IMPROPER_TYPES" in the documentation!
961  if (success) {
962  const string section_name{"CHARMM_NUM_IMPR_TYPES"};
963  const auto &current_section = toppar_map.find(section_name);
964  if (current_section != toppar_map.end()) {
965  toppar_data.NimprTypes = get<1>(current_section->second)[0].Int;
966  toppar_data.ImproperForceConstants.reserve(toppar_data.NimprTypes);
967  toppar_data.ImproperPhases.reserve(toppar_data.NimprTypes);
968  get<0>(current_section->second) = true;
969  } else {
970  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
971  }
972  }
973  // CHARMM_IMPROPER_FORCE_CONSTANT
974  if (success) {
975  const string section_name{"CHARMM_IMPROPER_FORCE_CONSTANT"};
976  const auto &current_section = toppar_map.find(section_name);
977  if (current_section != toppar_map.end()) {
978  success = parse_section(get<1>(current_section->second), toppar_data.NimprTypes,
979  toppar_data.ImproperForceConstants, section_name);
980  get<0>(current_section->second) = true;
981  } else {
982  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
983  }
984  }
985  // CHARMM_IMPROPER_PHASE
986  if (success) {
987  const string section_name{"CHARMM_IMPROPER_PHASE"};
988  const auto &current_section = toppar_map.find(section_name);
989  if (current_section != toppar_map.end()) {
990  success = parse_section(get<1>(current_section->second), toppar_data.NimprTypes,
991  toppar_data.ImproperPhases, section_name);
992  get<0>(current_section->second) = true;
993  } else {
994  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
995  }
996  }
997  // LENNARD_JONES_14_ACOEF
998  if (success) {
999  const int Nttyp = toppar_data.Ntypes * (toppar_data.Ntypes + 1) / 2;
1000  const string section_name{"LENNARD_JONES_14_ACOEF"};
1001  const auto &current_section = toppar_map.find(section_name);
1002  if (current_section != toppar_map.end()) {
1003  // this is not assigned before
1004  toppar_data.LJ14ACoefficients.reserve(Nttyp);
1005  success = parse_section(get<1>(current_section->second), Nttyp,
1006  toppar_data.LJ14ACoefficients, section_name);
1007  get<0>(current_section->second) = true;
1008  } else {
1009  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
1010  }
1011  }
1012  // LENNARD_JONES_14_BCOEF
1013  if (success) {
1014  const int Nttyp = toppar_data.Ntypes * (toppar_data.Ntypes + 1) / 2;
1015  const string section_name{"LENNARD_JONES_14_BCOEF"};
1016  const auto &current_section = toppar_map.find(section_name);
1017  if (current_section != toppar_map.end()) {
1018  // this is not assigned before
1019  toppar_data.LJ14BCoefficients.reserve(Nttyp);
1020  success = parse_section(get<1>(current_section->second), Nttyp,
1021  toppar_data.LJ14BCoefficients, section_name);
1022  get<0>(current_section->second) = true;
1023  } else {
1024  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
1025  }
1026  }
1027  }
1028  // Trying to parse CMAP terms
1029  // ff19SB also have CMAP terms, so CMAP is not CHARMM-specific.
1030  // However, in CHARMM FF they are prefixed with CHARMM_
1031  // CMAP_COUNT or CHARMM_CMAP_COUNT
1032  if (success) {
1033  const string section_name = toppar_data.IsCharmmFF
1034  ? string("CHARMM_CMAP_COUNT")
1035  : string("CMAP_COUNT");
1036  const auto &current_section = toppar_map.find(section_name);
1037  if (current_section != toppar_map.end()) {
1038  if (get<1>(current_section->second).size() == 2) {
1039  toppar_data.CMAPTermCount = get<1>(current_section->second)[0].Int;
1040  toppar_data.CMAPTypeCount = get<1>(current_section->second)[1].Int;
1041  toppar_data.HasCMAP = true;
1042  toppar_data.CMAPResolution.reserve(toppar_data.CMAPTypeCount);
1043  // CMAPParameter will be resized later after reading CMAPResolution
1044  toppar_data.CMAPParameter.reserve(toppar_data.CMAPTypeCount);
1045  toppar_data.CMAPIndex.reserve(toppar_data.CMAPTermCount * 6);
1046  get<0>(current_section->second) = true;
1047  } else {
1048  iout << iERROR << "Expect " << 2 << " fields but only "
1049  << get<1>(current_section->second).size() << " are in "
1050  << section_name.c_str() << "\n";
1051  success = false;
1052  }
1053  } else {
1054  toppar_data.HasCMAP = false;
1055  }
1056  }
1057  if (toppar_data.HasCMAP) {
1058  // CMAP_RESOLUTION or CHARMM_CMAP_RESOLUTION
1059  if (success) {
1060  const string section_name = toppar_data.IsCharmmFF
1061  ? string("CHARMM_CMAP_RESOLUTION")
1062  : string("CMAP_RESOLUTION");
1063  const auto &current_section = toppar_map.find(section_name);
1064  if (current_section != toppar_map.end()) {
1065  success = parse_section(get<1>(current_section->second),
1066  toppar_data.CMAPTypeCount, toppar_data.CMAPResolution,
1067  section_name);
1068  for (int i = 0; i < toppar_data.CMAPTypeCount; ++i) {
1069  toppar_data.CMAPParameter.push_back(vector<_REAL>());
1070  }
1071  get<0>(current_section->second) = true;
1072  } else {
1073  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
1074  }
1075  }
1076  // CHARMM_CMAP_PARAMETER_N or CMAP_PARAMETER_N
1077  if (success) {
1078  for (int i = 0; i < toppar_data.CMAPTypeCount; ++i) {
1079  string name_index = std::to_string(i + 1);
1080  if (name_index.size() < 2)
1081  name_index.insert(0, "0");
1082  const string section_name =
1083  (toppar_data.IsCharmmFF ? string("CHARMM_CMAP_PARAMETER_")
1084  : string("CMAP_PARAMETER_")) +
1085  name_index;
1086  const auto &current_section = toppar_map.find(section_name);
1087  if (current_section != toppar_map.end()) {
1088  const int resolution =
1089  toppar_data.CMAPResolution[i] * toppar_data.CMAPResolution[i];
1090  toppar_data.CMAPParameter[i].reserve(resolution);
1091  success = parse_section(get<1>(current_section->second), resolution,
1092  toppar_data.CMAPParameter[i], section_name);
1093  get<0>(current_section->second) = true;
1094  } else {
1095  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
1096  }
1097  }
1098  }
1099  // CMAP_INDEX
1100  if (success) {
1101  const string section_name = toppar_data.IsCharmmFF
1102  ? string("CHARMM_CMAP_INDEX")
1103  : string("CMAP_INDEX");
1104  const auto &current_section = toppar_map.find(section_name);
1105  if (current_section != toppar_map.end()) {
1106  success = parse_section(get<1>(current_section->second),
1107  toppar_data.CMAPTermCount * 6, toppar_data.CMAPIndex,
1108  section_name);
1109  get<0>(current_section->second) = true;
1110  } else {
1111  iout << iWARN << "Missing " << section_name.c_str() << " section!\n" << endi;
1112  }
1113  }
1114  }
1115  // NAMD requires the residue IDs for all atoms
1116  // so I need to fill AtomRes here before return
1117  if (success) {
1118  int res = 0;
1119  for (int i = 0; i < toppar_data.Natom; ++i) {
1120  if ((res+1 < toppar_data.Nres) && (i + 1 == toppar_data.Ipres[res+1]))
1121  ++res;
1122  toppar_data.AtomRes.push_back(res);
1123  }
1124  }
1125  // now we check and report any unknown sections
1126  for (auto it = toppar_map.begin(); it != toppar_map.end(); ++it) {
1127  if (get<0>(it->second) == false) {
1128  iout << iWARN << "Unknown section " << (it->first).c_str()
1129  << " found in the AMBER format paramter file!\n" << endi;
1130  }
1131  }
1132  return success;
1133 }
1134 
1135 bool parse_pointer(const vector<FortranData> &source, Ambertoppar &result) {
1136  bool success = true;
1137  if (source.size() > 30) {
1138  result.Natom = source[0].Int;
1139  result.Ntypes = source[1].Int;
1140  result.Nbonh = source[2].Int;
1141  result.Mbona = source[3].Int;
1142  result.Ntheth = source[4].Int;
1143  result.Mtheta = source[5].Int;
1144  result.Nphih = source[6].Int;
1145  result.Mphia = source[7].Int;
1146  result.Nhparm = source[8].Int;
1147  result.Nparm = source[9].Int;
1148  result.Nnb = source[10].Int;
1149  result.Nres = source[11].Int;
1150  result.Nbona = source[12].Int;
1151  result.Ntheta = source[13].Int;
1152  result.Nphia = source[14].Int;
1153  result.Numbnd = source[15].Int;
1154  result.Numang = source[16].Int;
1155  result.Nptra = source[17].Int;
1156  result.Natyp = source[18].Int;
1157  result.Nphb = source[19].Int;
1158  result.Ifpert = source[20].Int;
1159  result.Nbper = source[21].Int;
1160  result.Ngper = source[22].Int;
1161  result.Ndper = source[23].Int;
1162  result.Mbper = source[24].Int;
1163  result.Mgper = source[25].Int;
1164  result.Mdper = source[26].Int;
1165  result.Ifbox = source[27].Int;
1166  result.Nmxrs = source[28].Int;
1167  result.Ifcap = source[29].Int;
1168  result.Numextra = source[30].Int;
1169  // Ncopy may not be present
1170  result.Ncopy = 31 < source.size() ? source[31].Int : 0;
1171  // "allocate memory"
1172  const int Ntype2d = result.Ntypes * result.Ntypes;
1173  const int Nttyp = result.Ntypes * (result.Ntypes + 1) / 2;
1174  result.AtomNames.reserve(result.Natom);
1175  result.Charges.reserve(result.Natom);
1176  result.Masses.reserve(result.Natom);
1177  result.AtomNumbers.reserve(result.Natom);
1178  result.Iac.reserve(result.Natom);
1179  result.Iblo.reserve(result.Natom);
1180  result.Cno.reserve(Ntype2d);
1181  result.ResNames.reserve(result.Nres);
1182  result.Ipres.reserve(result.Nres);
1183  result.Rk.reserve(result.Numbnd);
1184  result.Req.reserve(result.Numbnd);
1185  result.Tk.reserve(result.Numang);
1186  result.Teq.reserve(result.Numang);
1187  result.Pk.reserve(result.Nptra);
1188  result.Pn.reserve(result.Nptra);
1189  result.Phase.reserve(result.Nptra);
1190  result.SceeScaleFactor.reserve(result.Nptra);
1191  result.ScnbScaleFactor.reserve(result.Nptra);
1192  result.Solty.reserve(result.Natyp);
1193  result.Cn1.reserve(Nttyp);
1194  result.Cn2.reserve(Nttyp);
1195  result.BondHAt1.reserve(result.Nbonh);
1196  result.BondHAt2.reserve(result.Nbonh);
1197  result.BondHNum.reserve(result.Nbonh);
1198  result.BondAt1.reserve(result.Nbona);
1199  result.BondAt2.reserve(result.Nbona);
1200  result.BondNum.reserve(result.Nbona);
1201  result.AngleHAt1.reserve(result.Ntheth);
1202  result.AngleHAt2.reserve(result.Ntheth);
1203  result.AngleHAt3.reserve(result.Ntheth);
1204  result.AngleHNum.reserve(result.Ntheth);
1205  result.AngleAt1.reserve(result.Ntheta);
1206  result.AngleAt2.reserve(result.Ntheta);
1207  result.AngleAt3.reserve(result.Ntheta);
1208  result.AngleNum.reserve(result.Ntheta);
1209  result.DihHAt1.reserve(result.Nphih);
1210  result.DihHAt2.reserve(result.Nphih);
1211  result.DihHAt3.reserve(result.Nphih);
1212  result.DihHAt4.reserve(result.Nphih);
1213  result.DihHNum.reserve(result.Nphih);
1214  result.DihAt1.reserve(result.Nphia);
1215  result.DihAt2.reserve(result.Nphia);
1216  result.DihAt3.reserve(result.Nphia);
1217  result.DihAt4.reserve(result.Nphia);
1218  result.DihNum.reserve(result.Nphia);
1219  result.ExclAt.reserve(result.Nnb);
1220  result.HB12.reserve(result.Nphb);
1221  result.HB10.reserve(result.Nphb);
1222  result.AtomSym.reserve(result.Natom);
1223  result.AtomTree.reserve(result.Natom);
1224  result.TreeJoin.reserve(result.Natom);
1225  result.AtomRes.reserve(result.Natom);
1226  result.RadiusSet.reserve(1);
1227  result.Radii.reserve(result.Natom);
1228  result.Screen.reserve(result.Natom);
1229  } else {
1230  success = false;
1231  }
1232  return success;
1233 }
1234 
1235 bool parse_section(const vector<FortranData> &source, const int &count,
1236  vector<string> &destination, const string &section_name) {
1237  bool success = true;
1238  if (int(source.size()) != count) {
1239  iout << iERROR << "Expect " << count << " fields but only "
1240  << source.size() << " are in "
1241  << section_name.c_str() << "\n";
1242  success = false;
1243  } else {
1244  for (int i = 0; i < count; ++i) {
1245  destination.push_back(source[i].String);
1246  }
1247  iout << iINFO << "Complete reading " << section_name.c_str() << " with " << count << " fields.\n" << endi;
1248  }
1249  return success;
1250 }
1251 
1252 bool parse_section(const vector<FortranData> &source, const int &count,
1253  vector<int> &destination, const string &section_name) {
1254  bool success = true;
1255  if (int(source.size()) != count) {
1256  iout << iERROR << "Expect " << count << " fields but only "
1257  << source.size() << " are in "
1258  << section_name.c_str() << "\n";
1259  success = false;
1260  } else {
1261  for (int i = 0; i < count; ++i) {
1262  destination.push_back(source[i].Int);
1263  }
1264  iout << iINFO << "Complete reading " << section_name.c_str() << " with " << count << " fields.\n" << endi;
1265  }
1266  return success;
1267 }
1268 
1269 bool parse_section(const vector<FortranData> &source, const int &count,
1270  vector<_REAL> &destination, const string &section_name) {
1271  bool success = true;
1272  if (int(source.size()) != count) {
1273  iout << iERROR << "Expect " << count << " fields but only "
1274  << source.size() << " are in "
1275  << section_name.c_str() << "\n";
1276  success = false;
1277  } else {
1278  for (int i = 0; i < count; ++i) {
1279  destination.push_back(source[i].Real);
1280  }
1281  iout << iINFO << "Complete reading " << section_name.c_str() << " with " << count << " fields.\n" << endi;
1282  }
1283  return success;
1284 }
1285 
1287  const vector<FortranData> &source, const int &count,
1288  std::initializer_list<std::reference_wrapper<vector<string>>> destination,
1289  const string &section_name) {
1290  bool success = true;
1291  const int total_size = destination.size() * count;
1292  if (int(source.size()) != total_size) {
1293  iout << iERROR << "Expect " << total_size << " fields but only "
1294  << source.size() << " are in "
1295  << section_name.c_str() << "\n";
1296  success = false;
1297  } else {
1298  const int stride = destination.size();
1299  for (int i = 0; i < count; ++i) {
1300  int j = 0;
1301  for (auto it_j = destination.begin(); it_j != destination.end(); ++it_j) {
1302  it_j->get().push_back(source[stride * i + j].String);
1303  ++j;
1304  }
1305  }
1306  iout << iINFO << "Complete reading " << section_name.c_str() << " with " << total_size << " fields.\n" << endi;
1307  }
1308  return success;
1309 }
1310 
1312  const vector<FortranData> &source, const int &count,
1313  std::initializer_list<std::reference_wrapper<vector<int>>> destination,
1314  const string &section_name) {
1315  bool success = true;
1316  const int total_size = destination.size() * count;
1317  if (int(source.size()) != total_size) {
1318  iout << iERROR << "Expect " << total_size << " fields but only "
1319  << source.size() << " are in "
1320  << section_name.c_str() << "\n";
1321  success = false;
1322  } else {
1323  const int stride = destination.size();
1324  for (int i = 0; i < count; ++i) {
1325  int j = 0;
1326  for (auto it_j = destination.begin(); it_j != destination.end(); ++it_j) {
1327  it_j->get().push_back(source[stride * i + j].Int);
1328  ++j;
1329  }
1330  }
1331  iout << iINFO << "Complete reading " << section_name.c_str() << " with " << total_size << " fields.\n" << endi;
1332  }
1333  return success;
1334 }
1335 
1337  const vector<FortranData> &source, const int &count,
1338  std::initializer_list<std::reference_wrapper<vector<_REAL>>> destination,
1339  const string &section_name) {
1340  bool success = true;
1341  const int total_size = destination.size() * count;
1342  if (int(source.size()) != total_size) {
1343  iout << iERROR << "Expect " << total_size << " fields but only "
1344  << source.size() << " are in "
1345  << section_name.c_str() << "\n";
1346  success = false;
1347  } else {
1348  const int stride = destination.size();
1349  for (int i = 0; i < count; ++i) {
1350  int j = 0;
1351  for (auto it_j = destination.begin(); it_j != destination.end(); ++it_j) {
1352  it_j->get().push_back(source[stride * i + j].Real);
1353  ++j;
1354  }
1355  }
1356  iout << iINFO << "Complete reading " << section_name.c_str() << " with " << total_size << " fields.\n" << endi;
1357  }
1358  return success;
1359 }
1360 
1361 Ambertoppar readparm(const char *filename) {
1362  AmberParm7Reader::AmberTopparMap Amber_toppar_map;
1364  iout << iINFO << "=== Start reading AMBER parm7 format file " << filename << " ===\n" << endi;
1365  if (AmberParm7Reader::read_amber_parm_stage1(filename, Amber_toppar_map)) {
1366  if (AmberParm7Reader::read_amber_parm_stage2(Amber_toppar_map, data) == true) {
1367  data.HasData = true;
1368  } else {
1369  iout << iERROR << "Failed to read AMBER parm7 file: " << filename << " at read_amber_parm_stage2\n" << endi;
1370  NAMD_die("Error in Ambertoppar readparm(const char *filename)\n");
1371  }
1372  } else {
1373  iout << iERROR << "Failed to read AMBER parm7 file: " << filename << " at read_amber_parm_stage1\n" << endi;
1374  NAMD_die("Error in Ambertoppar readparm(const char *filename)\n");
1375  }
1376  iout << iINFO << "=== Finish reading AMBER parm7 format file. ===\n" << endi;
1377  return data;
1378 }
1379 
1380 } // namespace AmberParm7Reader
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:81
void NAMD_err(const char *err_msg)
Definition: common.C:170
float Real
Definition: common.h:118
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
unordered_map< string, tuple< bool, vector< FortranData > > > AmberTopparMap
Definition: ReadAmberParm.h:87
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:82
void split_string_by_specifiers(const string &source, const vector< FortranFormatSpecifier > &specifiers, vector< FortranData > &destination)
This function splits a string by a set of fortran format specifiers.
#define iout
Definition: InfoStream.h:51
bool read_amber_parm_stage2(AmberTopparMap &toppar_map, Ambertoppar &toppar_data)
Read an AmberTopparMap into an Ambertoppar struct for NAMD.
void parse_fortran_format(const std::string &str, FortranFormatSpecifier &specifier)
Parse a single fortran format specifier.
Definition: ReadAmberParm.C:16
bool read_amber_parm_stage1(const char *filename, AmberTopparMap &toppar_map)
Read an AMBER topology file into an AmberTopparMap.
bool parse_pointer(const vector< FortranData > &source, Ambertoppar &result)
void split_string(const std::string &data, const std::string &delim, std::vector< std::string > &dest)
Definition: ReadAmberParm.C:77
vector< _REAL > LJ14BCoefficients
vector< _REAL > LJ14ACoefficients
void NAMD_die(const char *err_msg)
Definition: common.C:147
vector< vector< _REAL > > CMAPParameter
vector< _REAL > UreyBradleyEquilValues
Ambertoppar readparm(const char *filename)
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
vector< _REAL > ImproperForceConstants
vector< _REAL > UreyBradleyForceConstants
bool parse_section(const vector< FortranData > &source, const int &count, vector< string > &destination, const string &section_name)
Copy data from source to destination.