2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/gmxpreprocess/add_par.h"
53 #include "gromacs/gmxpreprocess/fflibutil.h"
54 #include "gromacs/gmxpreprocess/gen_ad.h"
55 #include "gromacs/gmxpreprocess/gen_vsite.h"
56 #include "gromacs/gmxpreprocess/grompp_impl.h"
57 #include "gromacs/gmxpreprocess/h_db.h"
58 #include "gromacs/gmxpreprocess/notset.h"
59 #include "gromacs/gmxpreprocess/pgutil.h"
60 #include "gromacs/gmxpreprocess/specbond.h"
61 #include "gromacs/gmxpreprocess/topdirs.h"
62 #include "gromacs/gmxpreprocess/topio.h"
63 #include "gromacs/gmxpreprocess/toputil.h"
64 #include "gromacs/math/functions.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/topology/residuetypes.h"
67 #include "gromacs/topology/symtab.h"
68 #include "gromacs/utility/binaryinformation.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/dir_separator.h"
71 #include "gromacs/utility/exceptions.h"
72 #include "gromacs/utility/fatalerror.h"
73 #include "gromacs/utility/futil.h"
74 #include "gromacs/utility/niceheader.h"
75 #include "gromacs/utility/path.h"
76 #include "gromacs/utility/programcontext.h"
77 #include "gromacs/utility/smalloc.h"
78 #include "gromacs/utility/strconvert.h"
79 #include "gromacs/utility/strdb.h"
80 #include "gromacs/utility/stringutil.h"
81 #include "gromacs/utility/textwriter.h"
83 #include "hackblock.h"
86 /* this must correspond to enum in pdb2top.h */
87 const char* hh
[ehisNR
] = { "HISD", "HISE", "HISH", "HIS1" };
89 static int missing_atoms(const PreprocessResidue
* rp
, int resind
, t_atoms
* at
, int i0
, int i
)
92 for (int j
= 0; j
< rp
->natom(); j
++)
94 const char* name
= *(rp
->atomname
[j
]);
96 for (int k
= i0
; k
< i
; k
++)
98 bFound
= (bFound
|| (gmx_strcasecmp(*(at
->atomname
[k
]), name
) == 0));
105 "atom %s is missing in residue %s %d in the pdb file\n",
106 name
, *(at
->resinfo
[resind
].name
), at
->resinfo
[resind
].nr
);
107 if (name
[0] == 'H' || name
[0] == 'h')
110 " You might need to add atom %s to the hydrogen database of "
111 "building block %s\n"
112 " in the file %s.hdb (see the manual)\n",
113 name
, *(at
->resinfo
[resind
].rtp
), rp
->filebase
.c_str());
115 fprintf(stderr
, "\n");
122 bool is_int(double x
)
124 const double tol
= 1e-4;
131 ix
= gmx::roundToInt(x
);
133 return (fabs(x
- ix
) < tol
);
136 static void choose_ff_impl(const char* ffsel
, char* forcefield
, int ff_maxlen
, char* ffdir
, int ffdir_maxlen
)
138 std::vector
<gmx::DataFileInfo
> ffdirs
= fflib_enumerate_forcefields();
139 const int nff
= ssize(ffdirs
);
141 /* Replace with unix path separators */
142 #if DIR_SEPARATOR != '/'
143 for (int i
= 0; i
< nff
; ++i
)
145 std::replace(ffdirs
[i
].dir
.begin(), ffdirs
[i
].dir
.end(), DIR_SEPARATOR
, '/');
149 /* Store the force field names in ffs */
150 std::vector
<std::string
> ffs
;
151 ffs
.reserve(ffdirs
.size());
152 for (int i
= 0; i
< nff
; ++i
)
154 ffs
.push_back(gmx::stripSuffixIfPresent(ffdirs
[i
].name
, fflib_forcefield_dir_ext()));
158 if (ffsel
!= nullptr)
163 for (int i
= 0; i
< nff
; ++i
)
167 /* Matching ff name */
171 if (ffdirs
[i
].dir
== ".")
188 "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
189 "current directory. Use interactive selection (not the -ff option) if\n"
190 "you would prefer a different one.\n",
195 std::string message
= gmx::formatString(
196 "Force field '%s' occurs in %d places, but not in "
197 "the current directory.\n"
198 "Run without the -ff switch and select the force "
199 "field interactively.",
201 GMX_THROW(gmx::InconsistentInputError(message
));
204 else if (nfound
== 0)
206 std::string message
= gmx::formatString(
207 "Could not find force field '%s' in current directory, "
208 "install tree or GMXLIB path.",
210 GMX_THROW(gmx::InconsistentInputError(message
));
215 std::vector
<std::string
> desc
;
216 desc
.reserve(ffdirs
.size());
217 for (int i
= 0; i
< nff
; ++i
)
219 std::string
docFileName(gmx::Path::join(ffdirs
[i
].dir
, ffdirs
[i
].name
, fflib_forcefield_doc()));
220 // TODO: Just try to open the file with a method that does not
221 // throw/bail out with a fatal error instead of multiple checks.
222 if (gmx::File::exists(docFileName
, gmx::File::returnFalseOnError
))
224 // TODO: Use a C++ API without such an intermediate/fixed-length buffer.
226 /* We don't use fflib_open, because we don't want printf's */
227 FILE* fp
= gmx_ffopen(docFileName
, "r");
228 get_a_line(fp
, buf
, STRLEN
);
230 desc
.emplace_back(buf
);
234 desc
.push_back(ffs
[i
]);
237 /* Order force fields from the same dir alphabetically
238 * and put deprecated force fields at the end.
240 for (int i
= 0; i
< nff
; ++i
)
242 for (int j
= i
+ 1; j
< nff
; ++j
)
244 if (ffdirs
[i
].dir
== ffdirs
[j
].dir
245 && ((desc
[i
][0] == '[' && desc
[j
][0] != '[')
246 || ((desc
[i
][0] == '[' || desc
[j
][0] != '[')
247 && gmx_strcasecmp(desc
[i
].c_str(), desc
[j
].c_str()) > 0)))
249 std::swap(ffdirs
[i
].name
, ffdirs
[j
].name
);
250 std::swap(ffs
[i
], ffs
[j
]);
251 std::swap(desc
[i
], desc
[j
]);
256 printf("\nSelect the Force Field:\n");
257 for (int i
= 0; i
< nff
; ++i
)
259 if (i
== 0 || ffdirs
[i
- 1].dir
!= ffdirs
[i
].dir
)
261 if (ffdirs
[i
].dir
== ".")
263 printf("From current directory:\n");
267 printf("From '%s':\n", ffdirs
[i
].dir
.c_str());
270 printf("%2d: %s\n", i
+ 1, desc
[i
].c_str());
274 // TODO: Add a C++ API for this.
279 pret
= fgets(buf
, STRLEN
, stdin
);
283 sel
= strtol(buf
, nullptr, 10);
286 } while (pret
== nullptr || (sel
< 0) || (sel
>= nff
));
288 /* Check for a current limitation of the fflib code.
289 * It will always read from the first ff directory in the list.
290 * This check assumes that the order of ffs matches the order
291 * in which fflib_open searches ff library files.
293 for (int i
= 0; i
< sel
; i
++)
295 if (ffs
[i
] == ffs
[sel
])
297 std::string message
= gmx::formatString(
298 "Can only select the first of multiple force "
299 "field entries with directory name '%s%s' in "
300 "the list. If you want to use the next entry, "
301 "run pdb2gmx in a different directory, set GMXLIB "
302 "to point to the desired force field first, and/or "
303 "rename or move the force field directory present "
304 "in the current working directory.",
305 ffs
[sel
].c_str(), fflib_forcefield_dir_ext());
306 GMX_THROW(gmx::NotImplementedError(message
));
315 if (ffs
[sel
].length() >= static_cast<size_t>(ff_maxlen
))
317 std::string message
= gmx::formatString("Length of force field name (%d) >= maxlen (%d)",
318 static_cast<int>(ffs
[sel
].length()), ff_maxlen
);
319 GMX_THROW(gmx::InvalidInputError(message
));
321 strcpy(forcefield
, ffs
[sel
].c_str());
324 if (ffdirs
[sel
].bFromDefaultDir
)
326 ffpath
= ffdirs
[sel
].name
;
330 ffpath
= gmx::Path::join(ffdirs
[sel
].dir
, ffdirs
[sel
].name
);
332 if (ffpath
.length() >= static_cast<size_t>(ffdir_maxlen
))
334 std::string message
= gmx::formatString("Length of force field dir (%d) >= maxlen (%d)",
335 static_cast<int>(ffpath
.length()), ffdir_maxlen
);
336 GMX_THROW(gmx::InvalidInputError(message
));
338 strcpy(ffdir
, ffpath
.c_str());
341 void choose_ff(const char* ffsel
, char* forcefield
, int ff_maxlen
, char* ffdir
, int ffdir_maxlen
)
345 choose_ff_impl(ffsel
, forcefield
, ff_maxlen
, ffdir
, ffdir_maxlen
);
347 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
350 void choose_watermodel(const char* wmsel
, const char* ffdir
, char** watermodel
)
352 const char* fn_watermodels
= "watermodels.dat";
359 if (strcmp(wmsel
, "none") == 0)
361 *watermodel
= nullptr;
365 else if (strcmp(wmsel
, "select") != 0)
367 *watermodel
= gmx_strdup(wmsel
);
372 std::string filename
= gmx::Path::join(ffdir
, fn_watermodels
);
373 if (!fflib_fexist(filename
))
375 fprintf(stderr
, "No file '%s' found, will not include a water model\n", fn_watermodels
);
376 *watermodel
= nullptr;
381 fp
= fflib_open(filename
);
382 printf("\nSelect the Water Model:\n");
385 while (get_a_line(fp
, buf
, STRLEN
))
387 srenew(model
, nwm
+ 1);
388 snew(model
[nwm
], STRLEN
);
389 sscanf(buf
, "%s%n", model
[nwm
], &i
);
393 fprintf(stderr
, "%2d: %s\n", nwm
+ 1, buf
+ i
);
402 fprintf(stderr
, "%2d: %s\n", nwm
+ 1, "None");
407 pret
= fgets(buf
, STRLEN
, stdin
);
411 sel
= strtol(buf
, nullptr, 10);
414 } while (pret
== nullptr || sel
< 0 || sel
> nwm
);
418 *watermodel
= nullptr;
422 *watermodel
= gmx_strdup(model
[sel
]);
425 for (i
= 0; i
< nwm
; i
++)
432 static int name2type(t_atoms
* at
, int** cgnr
, gmx::ArrayRef
<const PreprocessResidue
> usedPpResidues
, ResidueType
* rt
)
434 int i
, j
, prevresind
, i0
, prevcg
, cg
, curcg
;
450 for (i
= 0; (i
< at
->nr
); i
++)
453 if (at
->atom
[i
].resind
!= resind
)
455 resind
= at
->atom
[i
].resind
;
456 bool bProt
= rt
->namedResidueHasType(*(at
->resinfo
[resind
].name
), "Protein");
457 bNterm
= bProt
&& (resind
== 0);
460 nmissat
+= missing_atoms(&usedPpResidues
[prevresind
], prevresind
, at
, i0
, i
);
464 if (at
->atom
[i
].m
== 0)
468 name
= *(at
->atomname
[i
]);
469 j
= search_jtype(usedPpResidues
[resind
], name
, bNterm
);
470 at
->atom
[i
].type
= usedPpResidues
[resind
].atom
[j
].type
;
471 at
->atom
[i
].q
= usedPpResidues
[resind
].atom
[j
].q
;
472 at
->atom
[i
].m
= usedPpResidues
[resind
].atom
[j
].m
;
473 cg
= usedPpResidues
[resind
].cgnr
[j
];
474 /* A charge group number -1 signals a separate charge group
477 if ((cg
== -1) || (cg
!= prevcg
) || (resind
!= prevresind
))
493 at
->atom
[i
].typeB
= at
->atom
[i
].type
;
494 at
->atom
[i
].qB
= at
->atom
[i
].q
;
495 at
->atom
[i
].mB
= at
->atom
[i
].m
;
497 nmissat
+= missing_atoms(&usedPpResidues
[resind
], resind
, at
, i0
, i
);
502 static void print_top_heavy_H(FILE* out
, real mHmult
)
506 fprintf(out
, "; Using deuterium instead of hydrogen\n\n");
508 else if (mHmult
== 4.0)
510 fprintf(out
, "#define HEAVY_H\n\n");
512 else if (mHmult
!= 1.0)
515 "WARNING: unsupported proton mass multiplier (%g) "
521 void print_top_comment(FILE* out
, const char* filename
, const char* ffdir
, bool bITP
)
523 char ffdir_parent
[STRLEN
];
528 gmx::TextWriter
writer(out
);
529 gmx::niceHeader(&writer
, filename
, ';');
530 writer
.writeLine(gmx::formatString(";\tThis is a %s topology file", bITP
? "include" : "standalone"));
531 writer
.writeLine(";");
533 gmx::BinaryInformationSettings settings
;
534 settings
.generatedByHeader(true);
535 settings
.linePrefix(";\t");
536 gmx::printBinaryInformation(&writer
, gmx::getProgramContext(), settings
);
538 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
540 if (strchr(ffdir
, '/') == nullptr)
542 fprintf(out
, ";\tForce field was read from the standard GROMACS share directory.\n;\n\n");
544 else if (ffdir
[0] == '.')
547 ";\tForce field was read from current directory or a relative path - path "
552 strncpy(ffdir_parent
, ffdir
, STRLEN
- 1);
553 ffdir_parent
[STRLEN
- 1] = '\0'; /*make sure it is 0-terminated even for long string*/
554 p
= strrchr(ffdir_parent
, '/');
559 ";\tForce field data was read from:\n"
563 ";\tThis might be a non-standard force field location. When you use this topology, "
565 ";\tforce field must either be present in the current directory, or the location\n"
566 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file "
572 void print_top_header(FILE* out
, const char* filename
, bool bITP
, const char* ffdir
, real mHmult
)
576 print_top_comment(out
, filename
, ffdir
, bITP
);
578 print_top_heavy_H(out
, mHmult
);
579 fprintf(out
, "; Include forcefield parameters\n");
581 p
= strrchr(ffdir
, '/');
582 p
= (ffdir
[0] == '.' || p
== nullptr) ? ffdir
: p
+ 1;
584 fprintf(out
, "#include \"%s/%s\"\n\n", p
, fflib_forcefield_itp());
587 static void print_top_posre(FILE* out
, const char* pr
)
589 fprintf(out
, "; Include Position restraint file\n");
590 fprintf(out
, "#ifdef POSRES\n");
591 fprintf(out
, "#include \"%s\"\n", pr
);
592 fprintf(out
, "#endif\n\n");
595 static void print_top_water(FILE* out
, const char* ffdir
, const char* water
)
600 fprintf(out
, "; Include water topology\n");
602 p
= strrchr(ffdir
, '/');
603 p
= (ffdir
[0] == '.' || p
== nullptr) ? ffdir
: p
+ 1;
604 fprintf(out
, "#include \"%s/%s.itp\"\n", p
, water
);
607 fprintf(out
, "#ifdef POSRES_WATER\n");
608 fprintf(out
, "; Position restraint for each water oxygen\n");
609 fprintf(out
, "[ position_restraints ]\n");
610 fprintf(out
, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
611 fprintf(out
, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
612 fprintf(out
, "#endif\n");
615 sprintf(buf
, "%s/ions.itp", p
);
617 if (fflib_fexist(buf
))
619 fprintf(out
, "; Include topology for ions\n");
620 fprintf(out
, "#include \"%s\"\n", buf
);
625 static void print_top_system(FILE* out
, const char* title
)
627 fprintf(out
, "[ %s ]\n", dir2str(Directive::d_system
));
628 fprintf(out
, "; Name\n");
629 fprintf(out
, "%s\n\n", title
[0] ? title
: "Protein");
632 void print_top_mols(FILE* out
,
636 gmx::ArrayRef
<const std::string
> incls
,
637 gmx::ArrayRef
<const t_mols
> mols
)
642 fprintf(out
, "; Include chain topologies\n");
643 for (const auto& incl
: incls
)
645 fprintf(out
, "#include \"%s\"\n", gmx::Path::getFilename(incl
).c_str());
652 print_top_water(out
, ffdir
, water
);
654 print_top_system(out
, title
);
658 fprintf(out
, "[ %s ]\n", dir2str(Directive::d_molecules
));
659 fprintf(out
, "; %-15s %5s\n", "Compound", "#mols");
660 for (const auto& mol
: mols
)
662 fprintf(out
, "%-15s %5d\n", mol
.name
.c_str(), mol
.nr
);
667 void write_top(FILE* out
,
673 gmx::ArrayRef
<const InteractionsOfType
> plist
,
675 PreprocessingAtomTypes
* atype
,
678 /* NOTE: nrexcl is not the size of *excl! */
680 if (at
&& atype
&& cgnr
)
682 fprintf(out
, "[ %s ]\n", dir2str(Directive::d_moleculetype
));
683 fprintf(out
, "; %-15s %5s\n", "Name", "nrexcl");
684 fprintf(out
, "%-15s %5d\n\n", molname
? molname
: "Protein", nrexcl
);
686 print_atoms(out
, atype
, at
, cgnr
, bRTPresname
);
687 print_bondeds(out
, at
->nr
, Directive::d_bonds
, F_BONDS
, bts
[ebtsBONDS
], plist
);
688 print_bondeds(out
, at
->nr
, Directive::d_constraints
, F_CONSTR
, 0, plist
);
689 print_bondeds(out
, at
->nr
, Directive::d_constraints
, F_CONSTRNC
, 0, plist
);
690 print_bondeds(out
, at
->nr
, Directive::d_pairs
, F_LJ14
, 0, plist
);
691 print_excl(out
, at
->nr
, excls
);
692 print_bondeds(out
, at
->nr
, Directive::d_angles
, F_ANGLES
, bts
[ebtsANGLES
], plist
);
693 print_bondeds(out
, at
->nr
, Directive::d_dihedrals
, F_PDIHS
, bts
[ebtsPDIHS
], plist
);
694 print_bondeds(out
, at
->nr
, Directive::d_dihedrals
, F_IDIHS
, bts
[ebtsIDIHS
], plist
);
695 print_bondeds(out
, at
->nr
, Directive::d_cmap
, F_CMAP
, bts
[ebtsCMAP
], plist
);
696 print_bondeds(out
, at
->nr
, Directive::d_polarization
, F_POLARIZATION
, 0, plist
);
697 print_bondeds(out
, at
->nr
, Directive::d_thole_polarization
, F_THOLE_POL
, 0, plist
);
698 print_bondeds(out
, at
->nr
, Directive::d_vsites2
, F_VSITE2
, 0, plist
);
699 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3
, 0, plist
);
700 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3FD
, 0, plist
);
701 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3FAD
, 0, plist
);
702 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3OUT
, 0, plist
);
703 print_bondeds(out
, at
->nr
, Directive::d_vsites4
, F_VSITE4FD
, 0, plist
);
704 print_bondeds(out
, at
->nr
, Directive::d_vsites4
, F_VSITE4FDN
, 0, plist
);
708 print_top_posre(out
, pr
);
714 static void do_ssbonds(InteractionsOfType
* ps
,
716 gmx::ArrayRef
<const DisulfideBond
> ssbonds
,
719 for (const auto& bond
: ssbonds
)
721 int ri
= bond
.firstResidue
;
722 int rj
= bond
.secondResidue
;
723 int ai
= search_res_atom(bond
.firstAtom
.c_str(), ri
, atoms
, "special bond", bAllowMissing
);
724 int aj
= search_res_atom(bond
.secondAtom
.c_str(), rj
, atoms
, "special bond", bAllowMissing
);
725 if ((ai
== -1) || (aj
== -1))
727 gmx_fatal(FARGS
, "Trying to make impossible special bond (%s-%s)!",
728 bond
.firstAtom
.c_str(), bond
.secondAtom
.c_str());
730 add_param(ps
, ai
, aj
, {}, nullptr);
734 static void at2bonds(InteractionsOfType
* psb
,
735 gmx::ArrayRef
<MoleculePatchDatabase
> globalPatches
,
737 gmx::ArrayRef
<const gmx::RVec
> x
,
739 real short_bond_dist
)
741 real long_bond_dist2
, short_bond_dist2
;
744 long_bond_dist2
= gmx::square(long_bond_dist
);
745 short_bond_dist2
= gmx::square(short_bond_dist
);
756 fprintf(stderr
, "Making bonds...\n");
758 for (int resind
= 0; (resind
< atoms
->nres
) && (i
< atoms
->nr
); resind
++)
760 /* add bonds from list of bonded interactions */
761 for (const auto& patch
: globalPatches
[resind
].rb
[ebtsBONDS
].b
)
763 /* Unfortunately we can not issue errors or warnings
764 * for missing atoms in bonds, as the hydrogens and terminal atoms
765 * have not been added yet.
767 int ai
= search_atom(patch
.ai().c_str(), i
, atoms
, ptr
, TRUE
);
768 int aj
= search_atom(patch
.aj().c_str(), i
, atoms
, ptr
, TRUE
);
769 if (ai
!= -1 && aj
!= -1)
771 real dist2
= distance2(x
[ai
], x
[aj
]);
772 if (dist2
> long_bond_dist2
)
775 fprintf(stderr
, "Warning: Long Bond (%d-%d = %g nm)\n", ai
+ 1, aj
+ 1,
778 else if (dist2
< short_bond_dist2
)
780 fprintf(stderr
, "Warning: Short Bond (%d-%d = %g nm)\n", ai
+ 1, aj
+ 1,
783 add_param(psb
, ai
, aj
, {}, patch
.s
.c_str());
786 /* add bonds from list of hacks (each added atom gets a bond) */
787 while ((i
< atoms
->nr
) && (atoms
->atom
[i
].resind
== resind
))
789 for (const auto& patch
: globalPatches
[resind
].hack
)
791 if ((patch
.tp
> 0 || patch
.type() == MoleculePatchType::Add
)
792 && patch
.a
[0] == *(atoms
->atomname
[i
]))
796 case 9: /* COOH terminus */
797 add_param(psb
, i
, i
+ 1, {}, nullptr); /* C-O */
798 add_param(psb
, i
, i
+ 2, {}, nullptr); /* C-OA */
799 add_param(psb
, i
+ 2, i
+ 3, {}, nullptr); /* OA-H */
802 for (int k
= 0; (k
< patch
.nr
); k
++)
804 add_param(psb
, i
, i
+ k
+ 1, {}, nullptr);
811 /* we're now at the start of the next residue */
815 static bool pcompar(const InteractionOfType
& a
, const InteractionOfType
& b
)
819 if ((d
= a
.ai() - b
.ai()) != 0)
823 else if ((d
= a
.aj() - b
.aj()) != 0)
829 return a
.interactionTypeName().length() > b
.interactionTypeName().length();
833 static void clean_bonds(InteractionsOfType
* ps
)
838 for (auto& bond
: ps
->interactionTypes
)
842 std::sort(ps
->interactionTypes
.begin(), ps
->interactionTypes
.end(), pcompar
);
844 /* remove doubles, keep the first one always. */
845 int oldNumber
= ps
->size();
846 for (auto parm
= ps
->interactionTypes
.begin() + 1; parm
!= ps
->interactionTypes
.end();)
848 auto prev
= parm
- 1;
849 if (parm
->ai() == prev
->ai() && parm
->aj() == prev
->aj())
851 parm
= ps
->interactionTypes
.erase(parm
);
858 fprintf(stderr
, "Number of bonds was %d, now %zu\n", oldNumber
, ps
->size());
862 fprintf(stderr
, "No bonds\n");
866 void print_sums(const t_atoms
* atoms
, bool bSystem
)
874 where
= " in system";
883 for (i
= 0; (i
< atoms
->nr
); i
++)
885 m
+= atoms
->atom
[i
].m
;
886 qtot
+= atoms
->atom
[i
].q
;
889 fprintf(stderr
, "Total mass%s %.3f a.m.u.\n", where
, m
);
890 fprintf(stderr
, "Total charge%s %.3f e\n", where
, qtot
);
893 static void check_restp_type(const char* name
, int t1
, int t2
)
897 gmx_fatal(FARGS
, "Residues in one molecule have a different '%s' type: %d and %d", name
, t1
, t2
);
901 static void check_restp_types(const PreprocessResidue
& r0
, const PreprocessResidue
& r1
)
903 check_restp_type("all dihedrals", static_cast<int>(r0
.bKeepAllGeneratedDihedrals
),
904 static_cast<int>(r1
.bKeepAllGeneratedDihedrals
));
905 check_restp_type("nrexcl", r0
.nrexcl
, r1
.nrexcl
);
906 check_restp_type("HH14", static_cast<int>(r0
.bGenerateHH14Interactions
),
907 static_cast<int>(r1
.bGenerateHH14Interactions
));
908 check_restp_type("remove dihedrals", static_cast<int>(r0
.bRemoveDihedralIfWithImproper
),
909 static_cast<int>(r1
.bRemoveDihedralIfWithImproper
));
911 for (int i
= 0; i
< ebtsNR
; i
++)
913 check_restp_type(btsNames
[i
], r0
.rb
[i
].type
, r1
.rb
[i
].type
);
917 static void add_atom_to_restp(PreprocessResidue
* usedPpResidues
,
920 const MoleculePatch
* patch
)
923 for (int k
= 0; k
< patch
->nr
; k
++)
925 /* set counter in atomname */
926 std::string buf
= patch
->nname
;
929 buf
.append(gmx::formatString("%d", k
+ 1));
931 usedPpResidues
->atomname
.insert(usedPpResidues
->atomname
.begin() + at_start
+ 1 + k
,
932 put_symtab(symtab
, buf
.c_str()));
933 usedPpResidues
->atom
.insert(usedPpResidues
->atom
.begin() + at_start
+ 1 + k
, patch
->atom
.back());
934 if (patch
->cgnr
!= NOTSET
)
936 usedPpResidues
->cgnr
.insert(usedPpResidues
->cgnr
.begin() + at_start
+ 1 + k
, patch
->cgnr
);
940 usedPpResidues
->cgnr
.insert(usedPpResidues
->cgnr
.begin() + at_start
+ 1 + k
,
941 usedPpResidues
->cgnr
[at_start
]);
946 void get_hackblocks_rtp(std::vector
<MoleculePatchDatabase
>* globalPatches
,
947 std::vector
<PreprocessResidue
>* usedPpResidues
,
948 gmx::ArrayRef
<const PreprocessResidue
> rtpFFDB
,
953 gmx::ArrayRef
<MoleculePatchDatabase
*> ntdb
,
954 gmx::ArrayRef
<MoleculePatchDatabase
*> ctdb
,
955 gmx::ArrayRef
<const int> rn
,
956 gmx::ArrayRef
<const int> rc
,
962 globalPatches
->resize(nres
);
963 usedPpResidues
->clear();
964 /* first the termini */
965 for (int i
= 0; i
< nterpairs
; i
++)
967 if (rn
[i
] >= 0 && ntdb
[i
] != nullptr)
969 copyModificationBlocks(*ntdb
[i
], &globalPatches
->at(rn
[i
]));
971 if (rc
[i
] >= 0 && ctdb
[i
] != nullptr)
973 mergeAtomAndBondModifications(*ctdb
[i
], &globalPatches
->at(rc
[i
]));
977 /* then the whole rtp */
978 for (int i
= 0; i
< nres
; i
++)
980 /* Here we allow a mismatch of one character when looking for the rtp entry.
981 * For such a mismatch there should be only one mismatching name.
982 * This is mainly useful for small molecules such as ions.
983 * Note that this will usually not work for protein, DNA and RNA,
984 * since there the residue names should be listed in residuetypes.dat
985 * and an error will have been generated earlier in the process.
987 key
= *resinfo
[i
].rtp
;
989 resinfo
[i
].rtp
= put_symtab(symtab
, searchResidueDatabase(key
, rtpFFDB
).c_str());
990 auto res
= getDatabaseEntry(*resinfo
[i
].rtp
, rtpFFDB
);
991 usedPpResidues
->push_back(PreprocessResidue());
992 PreprocessResidue
* newentry
= &usedPpResidues
->back();
993 copyPreprocessResidues(*res
, newentry
, symtab
);
995 /* Check that we do not have different bonded types in one molecule */
996 check_restp_types(usedPpResidues
->front(), *newentry
);
999 for (int j
= 0; j
< nterpairs
&& tern
== -1; j
++)
1007 for (int j
= 0; j
< nterpairs
&& terc
== -1; j
++)
1014 bRM
= mergeBondedInteractionList(res
->rb
, globalPatches
->at(i
).rb
, tern
>= 0, terc
>= 0);
1016 if (bRM
&& ((tern
>= 0 && ntdb
[tern
] == nullptr) || (terc
>= 0 && ctdb
[terc
] == nullptr)))
1018 const char* errString
=
1019 "There is a dangling bond at at least one of the terminal ends and the force "
1020 "field does not provide terminal entries or files. Fix your terminal "
1021 "residues so that they match the residue database (.rtp) entries, or provide "
1022 "terminal database entries (.tdb).";
1025 fprintf(stderr
, "%s\n", errString
);
1029 gmx_fatal(FARGS
, "%s", errString
);
1032 else if (bRM
&& ((tern
>= 0 && ntdb
[tern
]->nhack() == 0) || (terc
>= 0 && ctdb
[terc
]->nhack() == 0)))
1034 const char* errString
=
1035 "There is a dangling bond at at least one of the terminal ends. Fix your "
1036 "coordinate file, add a new terminal database entry (.tdb), or select the "
1037 "proper existing terminal entry.";
1040 fprintf(stderr
, "%s\n", errString
);
1044 gmx_fatal(FARGS
, "%s", errString
);
1049 /* Apply patchs to t_restp entries
1050 i.e. add's and deletes from termini database will be
1051 added to/removed from residue topology
1052 we'll do this on one big dirty loop, so it won't make easy reading! */
1053 for (auto modifiedResidue
= globalPatches
->begin(); modifiedResidue
!= globalPatches
->end();
1056 const int pos
= std::distance(globalPatches
->begin(), modifiedResidue
);
1057 PreprocessResidue
* posres
= &usedPpResidues
->at(pos
);
1058 for (auto patch
= modifiedResidue
->hack
.begin(); patch
!= modifiedResidue
->hack
.end(); patch
++)
1062 /* find atom in restp */
1063 auto found
= std::find_if(posres
->atomname
.begin(), posres
->atomname
.end(),
1064 [&patch
](char** name
) {
1065 return (patch
->oname
.empty() && patch
->a
[0] == *name
)
1066 || (patch
->oname
== *name
);
1069 if (found
== posres
->atomname
.end())
1071 /* If we are doing an atom rename only, we don't need
1072 * to generate a fatal error if the old name is not found
1075 /* Deleting can happen also only on the input atoms,
1076 * not necessarily always on the rtp entry.
1078 if (patch
->type() == MoleculePatchType::Add
)
1081 "atom %s not found in buiding block %d%s "
1082 "while combining tdb and rtp",
1083 patch
->oname
.empty() ? patch
->a
[0].c_str() : patch
->oname
.c_str(),
1084 pos
+ 1, *resinfo
[pos
].rtp
);
1089 int l
= std::distance(posres
->atomname
.begin(), found
);
1090 switch (patch
->type())
1092 case MoleculePatchType::Add
:
1095 add_atom_to_restp(posres
, symtab
, l
, &(*patch
));
1098 case MoleculePatchType::Delete
:
1099 { /* we're deleting */
1100 posres
->atom
.erase(posres
->atom
.begin() + l
);
1101 posres
->atomname
.erase(posres
->atomname
.begin() + l
);
1102 posres
->cgnr
.erase(posres
->cgnr
.begin() + l
);
1105 case MoleculePatchType::Replace
:
1107 /* we're replacing */
1108 posres
->atom
[l
] = patch
->atom
.back();
1109 posres
->atomname
[l
] = put_symtab(symtab
, patch
->nname
.c_str());
1110 if (patch
->cgnr
!= NOTSET
)
1112 posres
->cgnr
[l
] = patch
->cgnr
;
1123 static bool atomname_cmp_nr(const char* anm
, const MoleculePatch
* patch
, int* nr
)
1130 return (gmx::equalCaseInsensitive(anm
, patch
->nname
));
1134 if (isdigit(anm
[strlen(anm
) - 1]))
1136 *nr
= anm
[strlen(anm
) - 1] - '0';
1142 if (*nr
<= 0 || *nr
> patch
->nr
)
1148 return (strlen(anm
) == patch
->nname
.length() + 1
1149 && gmx_strncasecmp(anm
, patch
->nname
.c_str(), patch
->nname
.length()) == 0);
1154 static bool match_atomnames_with_rtp_atom(t_atoms
* pdba
,
1155 gmx::ArrayRef
<gmx::RVec
> x
,
1158 PreprocessResidue
* localPpResidue
,
1159 const MoleculePatchDatabase
& singlePatch
,
1167 oldnm
= *pdba
->atomname
[atind
];
1168 resnr
= pdba
->resinfo
[pdba
->atom
[atind
].resind
].nr
;
1171 for (auto patch
= singlePatch
.hack
.begin(); patch
!= singlePatch
.hack
.end(); patch
++)
1173 if (patch
->type() == MoleculePatchType::Replace
&& gmx::equalCaseInsensitive(oldnm
, patch
->oname
))
1175 /* This is a replace entry. */
1176 /* Check if we are not replacing a replaced atom. */
1177 bool bReplaceReplace
= false;
1178 for (auto selfPatch
= singlePatch
.hack
.begin(); selfPatch
!= singlePatch
.hack
.end(); selfPatch
++)
1180 if (patch
!= selfPatch
&& selfPatch
->type() == MoleculePatchType::Replace
1181 && gmx::equalCaseInsensitive(selfPatch
->nname
, patch
->oname
))
1183 /* The replace in patch replaces an atom that
1184 * was already replaced in selfPatch, we do not want
1185 * second or higher level replaces at this stage.
1187 bReplaceReplace
= true;
1190 if (bReplaceReplace
)
1192 /* Skip this replace. */
1196 /* This atom still has the old name, rename it */
1197 std::string newnm
= patch
->nname
;
1198 auto found
= std::find_if(
1199 localPpResidue
->atomname
.begin(), localPpResidue
->atomname
.end(),
1200 [&newnm
](char** name
) { return gmx::equalCaseInsensitive(newnm
, *name
); });
1201 if (found
== localPpResidue
->atomname
.end())
1203 /* The new name is not present in the rtp.
1204 * We need to apply the replace also to the rtp entry.
1207 /* We need to find the add hack that can add this atom
1208 * to find out after which atom it should be added.
1210 bool bFoundInAdd
= false;
1211 for (auto rtpModification
= singlePatch
.hack
.begin();
1212 rtpModification
!= singlePatch
.hack
.end(); rtpModification
++)
1214 int k
= std::distance(localPpResidue
->atomname
.begin(), found
);
1215 std::string start_at
;
1216 if (rtpModification
->type() == MoleculePatchType::Add
1217 && atomname_cmp_nr(newnm
.c_str(), &(*rtpModification
), &anmnr
))
1221 start_at
= singlePatch
.hack
[k
].a
[0];
1225 start_at
= gmx::formatString("%s%d", singlePatch
.hack
[k
].nname
.c_str(),
1229 std::find_if(localPpResidue
->atomname
.begin(),
1230 localPpResidue
->atomname
.end(), [&start_at
](char** name
) {
1231 return gmx::equalCaseInsensitive(start_at
, *name
);
1233 if (found2
== localPpResidue
->atomname
.end())
1236 "Could not find atom '%s' in residue building block '%s' to "
1238 start_at
.c_str(), localPpResidue
->resname
.c_str(), newnm
.c_str());
1240 /* We can add the atom after atom start_nr */
1241 add_atom_to_restp(localPpResidue
, symtab
,
1242 std::distance(localPpResidue
->atomname
.begin(), found2
),
1252 "Could not find an 'add' entry for atom named '%s' corresponding to "
1253 "the 'replace' entry from atom name '%s' to '%s' for tdb or hdb "
1254 "database of residue type '%s'",
1255 newnm
.c_str(), patch
->oname
.c_str(), patch
->nname
.c_str(),
1256 localPpResidue
->resname
.c_str());
1262 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n", oldnm
,
1263 localPpResidue
->resname
.c_str(), resnr
, newnm
.c_str());
1265 /* Rename the atom in pdba */
1266 pdba
->atomname
[atind
] = put_symtab(symtab
, newnm
.c_str());
1268 else if (patch
->type() == MoleculePatchType::Delete
1269 && gmx::equalCaseInsensitive(oldnm
, patch
->oname
))
1271 /* This is a delete entry, check if this atom is present
1272 * in the rtp entry of this residue.
1274 auto found3
= std::find_if(
1275 localPpResidue
->atomname
.begin(), localPpResidue
->atomname
.end(),
1276 [&oldnm
](char** name
) { return gmx::equalCaseInsensitive(oldnm
, *name
); });
1277 if (found3
== localPpResidue
->atomname
.end())
1279 /* This atom is not present in the rtp entry,
1280 * delete is now from the input pdba.
1284 printf("Deleting atom '%s' in residue '%s' %d\n", oldnm
,
1285 localPpResidue
->resname
.c_str(), resnr
);
1287 /* We should free the atom name,
1288 * but it might be used multiple times in the symtab.
1289 * sfree(pdba->atomname[atind]);
1291 for (int k
= atind
+ 1; k
< pdba
->nr
; k
++)
1293 pdba
->atom
[k
- 1] = pdba
->atom
[k
];
1294 pdba
->atomname
[k
- 1] = pdba
->atomname
[k
];
1295 copy_rvec(x
[k
], x
[k
- 1]);
1306 void match_atomnames_with_rtp(gmx::ArrayRef
<PreprocessResidue
> usedPpResidues
,
1307 gmx::ArrayRef
<MoleculePatchDatabase
> globalPatches
,
1310 gmx::ArrayRef
<gmx::RVec
> x
,
1313 for (int i
= 0; i
< pdba
->nr
; i
++)
1315 const char* oldnm
= *pdba
->atomname
[i
];
1316 PreprocessResidue
* localPpResidue
= &usedPpResidues
[pdba
->atom
[i
].resind
];
1317 auto found
= std::find_if(
1318 localPpResidue
->atomname
.begin(), localPpResidue
->atomname
.end(),
1319 [&oldnm
](char** name
) { return gmx::equalCaseInsensitive(oldnm
, *name
); });
1320 if (found
== localPpResidue
->atomname
.end())
1322 /* Not found yet, check if we have to rename this atom */
1323 if (match_atomnames_with_rtp_atom(pdba
, x
, symtab
, i
, localPpResidue
,
1324 globalPatches
[pdba
->atom
[i
].resind
], bVerbose
))
1326 /* We deleted this atom, decrease the atom counter by 1. */
1333 #define NUM_CMAP_ATOMS 5
1334 static void gen_cmap(InteractionsOfType
* psb
, gmx::ArrayRef
<const PreprocessResidue
> usedPpResidues
, t_atoms
* atoms
)
1338 t_resinfo
* resinfo
= atoms
->resinfo
;
1339 int nres
= atoms
->nres
;
1340 int cmap_atomid
[NUM_CMAP_ATOMS
];
1341 int cmap_chainnum
= -1;
1352 fprintf(stderr
, "Making cmap torsions...\n");
1354 /* Most cmap entries use the N atom from the next residue, so the last
1355 * residue should not have its CMAP entry in that case, but for things like
1356 * dipeptides we sometimes define a complete CMAP entry inside a residue,
1357 * and in this case we need to process everything through the last residue.
1359 for (residx
= 0; residx
< nres
; residx
++)
1361 /* Add CMAP terms from the list of CMAP interactions */
1362 for (const auto& b
: usedPpResidues
[residx
].rb
[ebtsCMAP
].b
)
1364 bool bAddCMAP
= true;
1365 /* Loop over atoms in a candidate CMAP interaction and
1366 * check that they exist, are from the same chain and are
1367 * from residues labelled as protein. */
1368 for (int k
= 0; k
< NUM_CMAP_ATOMS
&& bAddCMAP
; k
++)
1370 /* Assign the pointer to the name of the next reference atom.
1371 * This can use -/+ labels to refer to previous/next residue.
1373 const char* pname
= b
.a
[k
].c_str();
1374 /* Skip this CMAP entry if it refers to residues before the
1375 * first or after the last residue.
1377 if (((strchr(pname
, '-') != nullptr) && (residx
== 0))
1378 || ((strchr(pname
, '+') != nullptr) && (residx
== nres
- 1)))
1384 cmap_atomid
[k
] = search_atom(pname
, i
, atoms
, ptr
, TRUE
);
1385 bAddCMAP
= bAddCMAP
&& (cmap_atomid
[k
] != -1);
1388 /* This break is necessary, because cmap_atomid[k]
1389 * == -1 cannot be safely used as an index
1390 * into the atom array. */
1393 int this_residue_index
= atoms
->atom
[cmap_atomid
[k
]].resind
;
1396 cmap_chainnum
= resinfo
[this_residue_index
].chainnum
;
1400 /* Does the residue for this atom have the same
1401 * chain number as the residues for previous
1403 bAddCMAP
= bAddCMAP
&& cmap_chainnum
== resinfo
[this_residue_index
].chainnum
;
1405 /* Here we used to check that the residuetype was protein and
1406 * disable bAddCMAP if that was not the case. However, some
1407 * special residues (say, alanine dipeptides) might not adhere
1408 * to standard naming, and if we start calling them normal
1409 * protein residues the user will be bugged to select termini.
1411 * Instead, I believe that the right course of action is to
1412 * keep the CMAP interaction if it is present in the RTP file
1413 * and we correctly identified all atoms (which is the case
1420 add_cmap_param(psb
, cmap_atomid
[0], cmap_atomid
[1], cmap_atomid
[2], cmap_atomid
[3],
1421 cmap_atomid
[4], b
.s
.c_str());
1425 if (residx
< nres
- 1)
1427 while (atoms
->atom
[i
].resind
< residx
+ 1)
1433 /* Start the next residue */
1436 static void scrub_charge_groups(int* cgnr
, int natoms
)
1440 for (i
= 0; i
< natoms
; i
++)
1447 void pdb2top(FILE* top_file
,
1448 const char* posre_fn
,
1449 const char* molname
,
1451 std::vector
<gmx::RVec
>* x
,
1452 PreprocessingAtomTypes
* atype
,
1454 gmx::ArrayRef
<const PreprocessResidue
> rtpFFDB
,
1455 gmx::ArrayRef
<PreprocessResidue
> usedPpResidues
,
1456 gmx::ArrayRef
<MoleculePatchDatabase
> globalPatches
,
1459 bool bVsiteAromatics
,
1462 gmx::ArrayRef
<const DisulfideBond
> ssbonds
,
1463 real long_bond_dist
,
1464 real short_bond_dist
,
1471 std::array
<InteractionsOfType
, F_NRE
> plist
;
1481 at2bonds(&(plist
[F_BONDS
]), globalPatches
, atoms
, *x
, long_bond_dist
, short_bond_dist
);
1483 /* specbonds: disulphide bonds & heme-his */
1484 do_ssbonds(&(plist
[F_BONDS
]), atoms
, ssbonds
, bAllowMissing
);
1486 nmissat
= name2type(atoms
, &cgnr
, usedPpResidues
, &rt
);
1491 fprintf(stderr
, "There were %d missing atoms in molecule %s\n", nmissat
, molname
);
1496 "There were %d missing atoms in molecule %s, if you want to use this "
1497 "incomplete topology anyhow, use the option -missing",
1502 /* Cleanup bonds (sort and rm doubles) */
1503 clean_bonds(&(plist
[F_BONDS
]));
1505 snew(vsite_type
, atoms
->nr
);
1506 for (i
= 0; i
< atoms
->nr
; i
++)
1508 vsite_type
[i
] = NOTSET
;
1512 if (bVsiteAromatics
)
1515 "The conversion of aromatic rings into virtual sites is deprecated "
1516 "and may be removed in a future version of GROMACS");
1518 /* determine which atoms will be vsites and add dummy masses
1519 also renumber atom numbers in plist[0..F_NRE]! */
1520 do_vsites(rtpFFDB
, atype
, atoms
, tab
, x
, plist
, &vsite_type
, &cgnr
, mHmult
, bVsiteAromatics
, ffdir
);
1523 /* Make Angles and Dihedrals */
1524 fprintf(stderr
, "Generating angles, dihedrals and pairs...\n");
1525 snew(excls
, atoms
->nr
);
1526 gen_pad(atoms
, usedPpResidues
, plist
, excls
, globalPatches
, bAllowMissing
);
1531 gen_cmap(&(plist
[F_CMAP
]), usedPpResidues
, atoms
);
1532 if (plist
[F_CMAP
].size() > 0)
1534 fprintf(stderr
, "There are %4zu cmap torsion pairs\n", plist
[F_CMAP
].size());
1538 /* set mass of all remaining hydrogen atoms */
1541 do_h_mass(&(plist
[F_BONDS
]), vsite_type
, atoms
, mHmult
, bDeuterate
);
1545 /* Cleanup bonds (sort and rm doubles) */
1546 /* clean_bonds(&(plist[F_BONDS]));*/
1549 "There are %4zu dihedrals, %4zu impropers, %4zu angles\n"
1550 " %4zu pairs, %4zu bonds and %4zu virtual sites\n",
1551 plist
[F_PDIHS
].size(), plist
[F_IDIHS
].size(), plist
[F_ANGLES
].size(),
1552 plist
[F_LJ14
].size(), plist
[F_BONDS
].size(),
1553 plist
[F_VSITE2
].size() + plist
[F_VSITE3
].size() + plist
[F_VSITE3FD
].size()
1554 + plist
[F_VSITE3FAD
].size() + plist
[F_VSITE3OUT
].size()
1555 + plist
[F_VSITE4FD
].size() + plist
[F_VSITE4FDN
].size());
1557 print_sums(atoms
, FALSE
);
1561 scrub_charge_groups(cgnr
, atoms
->nr
);
1566 for (i
= 0; i
< atoms
->nres
; i
++)
1568 atoms
->resinfo
[i
].nr
= i
+ 1;
1569 atoms
->resinfo
[i
].ic
= ' ';
1575 fprintf(stderr
, "Writing topology\n");
1576 /* We can copy the bonded types from the first restp,
1577 * since the types have to be identical for all residues in one molecule.
1579 for (i
= 0; i
< ebtsNR
; i
++)
1581 bts
[i
] = usedPpResidues
[0].rb
[i
].type
;
1583 write_top(top_file
, posre_fn
, molname
, atoms
, bRTPresname
, bts
, plist
, excls
, atype
, cgnr
,
1584 usedPpResidues
[0].nrexcl
);
1588 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1590 for (i
= 0; i
< atoms
->nr
; i
++)