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/logger.h"
75 #include "gromacs/utility/niceheader.h"
76 #include "gromacs/utility/path.h"
77 #include "gromacs/utility/programcontext.h"
78 #include "gromacs/utility/smalloc.h"
79 #include "gromacs/utility/strconvert.h"
80 #include "gromacs/utility/strdb.h"
81 #include "gromacs/utility/stringutil.h"
82 #include "gromacs/utility/textwriter.h"
84 #include "hackblock.h"
87 /* this must correspond to enum in pdb2top.h */
88 const char* hh
[ehisNR
] = { "HISD", "HISE", "HISH", "HIS1" };
90 static int missing_atoms(const PreprocessResidue
* rp
, int resind
, t_atoms
* at
, int i0
, int i
, const gmx::MDLogger
& logger
)
93 for (int j
= 0; j
< rp
->natom(); j
++)
95 const char* name
= *(rp
->atomname
[j
]);
97 for (int k
= i0
; k
< i
; k
++)
99 bFound
= (bFound
|| (gmx_strcasecmp(*(at
->atomname
[k
]), name
) == 0));
104 GMX_LOG(logger
.warning
)
106 .appendTextFormatted("atom %s is missing in residue %s %d in the pdb file",
107 name
, *(at
->resinfo
[resind
].name
), at
->resinfo
[resind
].nr
);
108 if (name
[0] == 'H' || name
[0] == 'h')
110 GMX_LOG(logger
.warning
)
112 .appendTextFormatted(
113 "You might need to add atom %s to the hydrogen database of "
115 "in the file %s.hdb (see the manual)",
116 name
, *(at
->resinfo
[resind
].rtp
), rp
->filebase
.c_str());
124 bool is_int(double x
)
126 const double tol
= 1e-4;
133 ix
= gmx::roundToInt(x
);
135 return (fabs(x
- ix
) < tol
);
138 static void choose_ff_impl(const char* ffsel
,
143 const gmx::MDLogger
& logger
)
145 std::vector
<gmx::DataFileInfo
> ffdirs
= fflib_enumerate_forcefields();
146 const int nff
= ssize(ffdirs
);
148 /* Replace with unix path separators */
149 #if DIR_SEPARATOR != '/'
150 for (int i
= 0; i
< nff
; ++i
)
152 std::replace(ffdirs
[i
].dir
.begin(), ffdirs
[i
].dir
.end(), DIR_SEPARATOR
, '/');
156 /* Store the force field names in ffs */
157 std::vector
<std::string
> ffs
;
158 ffs
.reserve(ffdirs
.size());
159 for (int i
= 0; i
< nff
; ++i
)
161 ffs
.push_back(gmx::stripSuffixIfPresent(ffdirs
[i
].name
, fflib_forcefield_dir_ext()));
165 if (ffsel
!= nullptr)
170 for (int i
= 0; i
< nff
; ++i
)
174 /* Matching ff name */
178 if (ffdirs
[i
].dir
== ".")
194 GMX_LOG(logger
.warning
)
196 .appendTextFormatted(
197 "Force field '%s' occurs in %d places. pdb2gmx is using the one in "
198 "the current directory. Use interactive selection "
199 "(not the -ff option) if you would prefer a different one.",
204 std::string message
= gmx::formatString(
205 "Force field '%s' occurs in %d places, but not in "
206 "the current directory.\n"
207 "Run without the -ff switch and select the force "
208 "field interactively.",
210 GMX_THROW(gmx::InconsistentInputError(message
));
213 else if (nfound
== 0)
215 std::string message
= gmx::formatString(
216 "Could not find force field '%s' in current directory, "
217 "install tree or GMXLIB path.",
219 GMX_THROW(gmx::InconsistentInputError(message
));
224 std::vector
<std::string
> desc
;
225 desc
.reserve(ffdirs
.size());
226 for (int i
= 0; i
< nff
; ++i
)
228 std::string
docFileName(gmx::Path::join(ffdirs
[i
].dir
, ffdirs
[i
].name
, fflib_forcefield_doc()));
229 // TODO: Just try to open the file with a method that does not
230 // throw/bail out with a fatal error instead of multiple checks.
231 if (gmx::File::exists(docFileName
, gmx::File::returnFalseOnError
))
233 // TODO: Use a C++ API without such an intermediate/fixed-length buffer.
235 /* We don't use fflib_open, because we don't want printf's */
236 FILE* fp
= gmx_ffopen(docFileName
, "r");
237 get_a_line(fp
, buf
, STRLEN
);
239 desc
.emplace_back(buf
);
243 desc
.push_back(ffs
[i
]);
246 /* Order force fields from the same dir alphabetically
247 * and put deprecated force fields at the end.
249 for (int i
= 0; i
< nff
; ++i
)
251 for (int j
= i
+ 1; j
< nff
; ++j
)
253 if (ffdirs
[i
].dir
== ffdirs
[j
].dir
254 && ((desc
[i
][0] == '[' && desc
[j
][0] != '[')
255 || ((desc
[i
][0] == '[' || desc
[j
][0] != '[')
256 && gmx_strcasecmp(desc
[i
].c_str(), desc
[j
].c_str()) > 0)))
258 std::swap(ffdirs
[i
].name
, ffdirs
[j
].name
);
259 std::swap(ffs
[i
], ffs
[j
]);
260 std::swap(desc
[i
], desc
[j
]);
265 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("Select the Force Field:");
266 for (int i
= 0; i
< nff
; ++i
)
268 if (i
== 0 || ffdirs
[i
- 1].dir
!= ffdirs
[i
].dir
)
270 if (ffdirs
[i
].dir
== ".")
274 .appendTextFormatted("From current directory:");
280 .appendTextFormatted("From '%s':", ffdirs
[i
].dir
.c_str());
283 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("%2d: %s", i
+ 1, desc
[i
].c_str());
287 // TODO: Add a C++ API for this.
292 pret
= fgets(buf
, STRLEN
, stdin
);
296 sel
= strtol(buf
, nullptr, 10);
299 } while (pret
== nullptr || (sel
< 0) || (sel
>= nff
));
301 /* Check for a current limitation of the fflib code.
302 * It will always read from the first ff directory in the list.
303 * This check assumes that the order of ffs matches the order
304 * in which fflib_open searches ff library files.
306 for (int i
= 0; i
< sel
; i
++)
308 if (ffs
[i
] == ffs
[sel
])
310 std::string message
= gmx::formatString(
311 "Can only select the first of multiple force "
312 "field entries with directory name '%s%s' in "
313 "the list. If you want to use the next entry, "
314 "run pdb2gmx in a different directory, set GMXLIB "
315 "to point to the desired force field first, and/or "
316 "rename or move the force field directory present "
317 "in the current working directory.",
318 ffs
[sel
].c_str(), fflib_forcefield_dir_ext());
319 GMX_THROW(gmx::NotImplementedError(message
));
328 if (ffs
[sel
].length() >= static_cast<size_t>(ff_maxlen
))
330 std::string message
= gmx::formatString("Length of force field name (%d) >= maxlen (%d)",
331 static_cast<int>(ffs
[sel
].length()), ff_maxlen
);
332 GMX_THROW(gmx::InvalidInputError(message
));
334 strcpy(forcefield
, ffs
[sel
].c_str());
337 if (ffdirs
[sel
].bFromDefaultDir
)
339 ffpath
= ffdirs
[sel
].name
;
343 ffpath
= gmx::Path::join(ffdirs
[sel
].dir
, ffdirs
[sel
].name
);
345 if (ffpath
.length() >= static_cast<size_t>(ffdir_maxlen
))
347 std::string message
= gmx::formatString("Length of force field dir (%d) >= maxlen (%d)",
348 static_cast<int>(ffpath
.length()), ffdir_maxlen
);
349 GMX_THROW(gmx::InvalidInputError(message
));
351 strcpy(ffdir
, ffpath
.c_str());
354 void choose_ff(const char* ffsel
, char* forcefield
, int ff_maxlen
, char* ffdir
, int ffdir_maxlen
, const gmx::MDLogger
& logger
)
358 choose_ff_impl(ffsel
, forcefield
, ff_maxlen
, ffdir
, ffdir_maxlen
, logger
);
360 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
363 void choose_watermodel(const char* wmsel
, const char* ffdir
, char** watermodel
, const gmx::MDLogger
& logger
)
365 const char* fn_watermodels
= "watermodels.dat";
372 if (strcmp(wmsel
, "none") == 0)
374 *watermodel
= nullptr;
378 else if (strcmp(wmsel
, "select") != 0)
380 *watermodel
= gmx_strdup(wmsel
);
385 std::string filename
= gmx::Path::join(ffdir
, fn_watermodels
);
386 if (!fflib_fexist(filename
))
388 GMX_LOG(logger
.warning
)
390 .appendTextFormatted("No file '%s' found, will not include a water model", fn_watermodels
);
391 *watermodel
= nullptr;
396 fp
= fflib_open(filename
);
397 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("Select the Water Model:");
400 while (get_a_line(fp
, buf
, STRLEN
))
402 srenew(model
, nwm
+ 1);
403 snew(model
[nwm
], STRLEN
);
404 sscanf(buf
, "%s%n", model
[nwm
], &i
);
408 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("%2d: %s", nwm
+ 1, buf
+ i
);
417 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("%2d: %s", nwm
+ 1, "None");
422 pret
= fgets(buf
, STRLEN
, stdin
);
426 sel
= strtol(buf
, nullptr, 10);
429 } while (pret
== nullptr || sel
< 0 || sel
> nwm
);
433 *watermodel
= nullptr;
437 *watermodel
= gmx_strdup(model
[sel
]);
440 for (i
= 0; i
< nwm
; i
++)
447 static int name2type(t_atoms
* at
,
449 gmx::ArrayRef
<const PreprocessResidue
> usedPpResidues
,
451 const gmx::MDLogger
& logger
)
453 int i
, j
, prevresind
, i0
, prevcg
, cg
, curcg
;
469 for (i
= 0; (i
< at
->nr
); i
++)
472 if (at
->atom
[i
].resind
!= resind
)
474 resind
= at
->atom
[i
].resind
;
475 bool bProt
= rt
->namedResidueHasType(*(at
->resinfo
[resind
].name
), "Protein");
476 bNterm
= bProt
&& (resind
== 0);
479 nmissat
+= missing_atoms(&usedPpResidues
[prevresind
], prevresind
, at
, i0
, i
, logger
);
483 if (at
->atom
[i
].m
== 0)
487 name
= *(at
->atomname
[i
]);
488 j
= search_jtype(usedPpResidues
[resind
], name
, bNterm
);
489 at
->atom
[i
].type
= usedPpResidues
[resind
].atom
[j
].type
;
490 at
->atom
[i
].q
= usedPpResidues
[resind
].atom
[j
].q
;
491 at
->atom
[i
].m
= usedPpResidues
[resind
].atom
[j
].m
;
492 cg
= usedPpResidues
[resind
].cgnr
[j
];
493 /* A charge group number -1 signals a separate charge group
496 if ((cg
== -1) || (cg
!= prevcg
) || (resind
!= prevresind
))
512 at
->atom
[i
].typeB
= at
->atom
[i
].type
;
513 at
->atom
[i
].qB
= at
->atom
[i
].q
;
514 at
->atom
[i
].mB
= at
->atom
[i
].m
;
516 nmissat
+= missing_atoms(&usedPpResidues
[resind
], resind
, at
, i0
, i
, logger
);
521 static void print_top_heavy_H(FILE* out
, real mHmult
)
525 fprintf(out
, "; Using deuterium instead of hydrogen\n\n");
527 else if (mHmult
== 4.0)
529 fprintf(out
, "#define HEAVY_H\n\n");
531 else if (mHmult
!= 1.0)
534 "WARNING: unsupported proton mass multiplier (%g) "
540 void print_top_comment(FILE* out
, const char* filename
, const char* ffdir
, bool bITP
)
542 char ffdir_parent
[STRLEN
];
547 gmx::TextWriter
writer(out
);
548 gmx::niceHeader(&writer
, filename
, ';');
549 writer
.writeLine(gmx::formatString(";\tThis is a %s topology file", bITP
? "include" : "standalone"));
550 writer
.writeLine(";");
552 gmx::BinaryInformationSettings settings
;
553 settings
.generatedByHeader(true);
554 settings
.linePrefix(";\t");
555 gmx::printBinaryInformation(&writer
, gmx::getProgramContext(), settings
);
557 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
559 if (strchr(ffdir
, '/') == nullptr)
561 fprintf(out
, ";\tForce field was read from the standard GROMACS share directory.\n;\n\n");
563 else if (ffdir
[0] == '.')
566 ";\tForce field was read from current directory or a relative path - path "
571 strncpy(ffdir_parent
, ffdir
, STRLEN
- 1);
572 ffdir_parent
[STRLEN
- 1] = '\0'; /*make sure it is 0-terminated even for long string*/
573 p
= strrchr(ffdir_parent
, '/');
578 ";\tForce field data was read from:\n"
582 ";\tThis might be a non-standard force field location. When you use this topology, "
584 ";\tforce field must either be present in the current directory, or the location\n"
585 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file "
591 void print_top_header(FILE* out
, const char* filename
, bool bITP
, const char* ffdir
, real mHmult
)
595 print_top_comment(out
, filename
, ffdir
, bITP
);
597 print_top_heavy_H(out
, mHmult
);
598 fprintf(out
, "; Include forcefield parameters\n");
600 p
= strrchr(ffdir
, '/');
601 p
= (ffdir
[0] == '.' || p
== nullptr) ? ffdir
: p
+ 1;
603 fprintf(out
, "#include \"%s/%s\"\n\n", p
, fflib_forcefield_itp());
606 static void print_top_posre(FILE* out
, const char* pr
)
608 fprintf(out
, "; Include Position restraint file\n");
609 fprintf(out
, "#ifdef POSRES\n");
610 fprintf(out
, "#include \"%s\"\n", pr
);
611 fprintf(out
, "#endif\n\n");
614 static void print_top_water(FILE* out
, const char* ffdir
, const char* water
)
619 fprintf(out
, "; Include water topology\n");
621 p
= strrchr(ffdir
, '/');
622 p
= (ffdir
[0] == '.' || p
== nullptr) ? ffdir
: p
+ 1;
623 fprintf(out
, "#include \"%s/%s.itp\"\n", p
, water
);
626 fprintf(out
, "#ifdef POSRES_WATER\n");
627 fprintf(out
, "; Position restraint for each water oxygen\n");
628 fprintf(out
, "[ position_restraints ]\n");
629 fprintf(out
, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
630 fprintf(out
, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
631 fprintf(out
, "#endif\n");
634 sprintf(buf
, "%s/ions.itp", p
);
636 if (fflib_fexist(buf
))
638 fprintf(out
, "; Include topology for ions\n");
639 fprintf(out
, "#include \"%s\"\n", buf
);
644 static void print_top_system(FILE* out
, const char* title
)
646 fprintf(out
, "[ %s ]\n", dir2str(Directive::d_system
));
647 fprintf(out
, "; Name\n");
648 fprintf(out
, "%s\n\n", title
[0] ? title
: "Protein");
651 void print_top_mols(FILE* out
,
655 gmx::ArrayRef
<const std::string
> incls
,
656 gmx::ArrayRef
<const t_mols
> mols
)
661 fprintf(out
, "; Include chain topologies\n");
662 for (const auto& incl
: incls
)
664 fprintf(out
, "#include \"%s\"\n", gmx::Path::getFilename(incl
).c_str());
671 print_top_water(out
, ffdir
, water
);
673 print_top_system(out
, title
);
677 fprintf(out
, "[ %s ]\n", dir2str(Directive::d_molecules
));
678 fprintf(out
, "; %-15s %5s\n", "Compound", "#mols");
679 for (const auto& mol
: mols
)
681 fprintf(out
, "%-15s %5d\n", mol
.name
.c_str(), mol
.nr
);
686 void write_top(FILE* out
,
692 gmx::ArrayRef
<const InteractionsOfType
> plist
,
694 PreprocessingAtomTypes
* atype
,
697 /* NOTE: nrexcl is not the size of *excl! */
699 if (at
&& atype
&& cgnr
)
701 fprintf(out
, "[ %s ]\n", dir2str(Directive::d_moleculetype
));
702 fprintf(out
, "; %-15s %5s\n", "Name", "nrexcl");
703 fprintf(out
, "%-15s %5d\n\n", molname
? molname
: "Protein", nrexcl
);
705 print_atoms(out
, atype
, at
, cgnr
, bRTPresname
);
706 print_bondeds(out
, at
->nr
, Directive::d_bonds
, F_BONDS
, bts
[ebtsBONDS
], plist
);
707 print_bondeds(out
, at
->nr
, Directive::d_constraints
, F_CONSTR
, 0, plist
);
708 print_bondeds(out
, at
->nr
, Directive::d_constraints
, F_CONSTRNC
, 0, plist
);
709 print_bondeds(out
, at
->nr
, Directive::d_pairs
, F_LJ14
, 0, plist
);
710 print_excl(out
, at
->nr
, excls
);
711 print_bondeds(out
, at
->nr
, Directive::d_angles
, F_ANGLES
, bts
[ebtsANGLES
], plist
);
712 print_bondeds(out
, at
->nr
, Directive::d_dihedrals
, F_PDIHS
, bts
[ebtsPDIHS
], plist
);
713 print_bondeds(out
, at
->nr
, Directive::d_dihedrals
, F_IDIHS
, bts
[ebtsIDIHS
], plist
);
714 print_bondeds(out
, at
->nr
, Directive::d_cmap
, F_CMAP
, bts
[ebtsCMAP
], plist
);
715 print_bondeds(out
, at
->nr
, Directive::d_polarization
, F_POLARIZATION
, 0, plist
);
716 print_bondeds(out
, at
->nr
, Directive::d_thole_polarization
, F_THOLE_POL
, 0, plist
);
717 print_bondeds(out
, at
->nr
, Directive::d_vsites2
, F_VSITE2
, 0, plist
);
718 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3
, 0, plist
);
719 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3FD
, 0, plist
);
720 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3FAD
, 0, plist
);
721 print_bondeds(out
, at
->nr
, Directive::d_vsites3
, F_VSITE3OUT
, 0, plist
);
722 print_bondeds(out
, at
->nr
, Directive::d_vsites4
, F_VSITE4FD
, 0, plist
);
723 print_bondeds(out
, at
->nr
, Directive::d_vsites4
, F_VSITE4FDN
, 0, plist
);
727 print_top_posre(out
, pr
);
733 static void do_ssbonds(InteractionsOfType
* ps
,
735 gmx::ArrayRef
<const DisulfideBond
> ssbonds
,
738 for (const auto& bond
: ssbonds
)
740 int ri
= bond
.firstResidue
;
741 int rj
= bond
.secondResidue
;
742 int ai
= search_res_atom(bond
.firstAtom
.c_str(), ri
, atoms
, "special bond", bAllowMissing
);
743 int aj
= search_res_atom(bond
.secondAtom
.c_str(), rj
, atoms
, "special bond", bAllowMissing
);
744 if ((ai
== -1) || (aj
== -1))
746 gmx_fatal(FARGS
, "Trying to make impossible special bond (%s-%s)!",
747 bond
.firstAtom
.c_str(), bond
.secondAtom
.c_str());
749 add_param(ps
, ai
, aj
, {}, nullptr);
753 static void at2bonds(InteractionsOfType
* psb
,
754 gmx::ArrayRef
<MoleculePatchDatabase
> globalPatches
,
756 gmx::ArrayRef
<const gmx::RVec
> x
,
758 real short_bond_dist
,
759 gmx::ArrayRef
<const int> cyclicBondsIndex
,
760 const gmx::MDLogger
& logger
)
762 real long_bond_dist2
, short_bond_dist2
;
765 long_bond_dist2
= gmx::square(long_bond_dist
);
766 short_bond_dist2
= gmx::square(short_bond_dist
);
777 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("Making bonds...");
779 for (int resind
= 0; (resind
< atoms
->nres
) && (i
< atoms
->nr
); resind
++)
781 /* add bonds from list of bonded interactions */
782 for (const auto& patch
: globalPatches
[resind
].rb
[ebtsBONDS
].b
)
784 /* Unfortunately we can not issue errors or warnings
785 * for missing atoms in bonds, as the hydrogens and terminal atoms
786 * have not been added yet.
788 int ai
= search_atom(patch
.ai().c_str(), i
, atoms
, ptr
, TRUE
, cyclicBondsIndex
);
789 int aj
= search_atom(patch
.aj().c_str(), i
, atoms
, ptr
, TRUE
, cyclicBondsIndex
);
790 if (ai
!= -1 && aj
!= -1)
792 real dist2
= distance2(x
[ai
], x
[aj
]);
793 if (dist2
> long_bond_dist2
)
796 GMX_LOG(logger
.warning
)
798 .appendTextFormatted("Long Bond (%d-%d = %g nm)", ai
+ 1, aj
+ 1,
801 else if (dist2
< short_bond_dist2
)
803 GMX_LOG(logger
.warning
)
805 .appendTextFormatted("Short Bond (%d-%d = %g nm)", ai
+ 1, aj
+ 1,
808 add_param(psb
, ai
, aj
, {}, patch
.s
.c_str());
811 /* add bonds from list of hacks (each added atom gets a bond) */
812 while ((i
< atoms
->nr
) && (atoms
->atom
[i
].resind
== resind
))
814 for (const auto& patch
: globalPatches
[resind
].hack
)
816 if ((patch
.tp
> 0 || patch
.type() == MoleculePatchType::Add
)
817 && patch
.a
[0] == *(atoms
->atomname
[i
]))
821 case 9: /* COOH terminus */
822 add_param(psb
, i
, i
+ 1, {}, nullptr); /* C-O */
823 add_param(psb
, i
, i
+ 2, {}, nullptr); /* C-OA */
824 add_param(psb
, i
+ 2, i
+ 3, {}, nullptr); /* OA-H */
827 for (int k
= 0; (k
< patch
.nr
); k
++)
829 add_param(psb
, i
, i
+ k
+ 1, {}, nullptr);
836 /* we're now at the start of the next residue */
840 static bool pcompar(const InteractionOfType
& a
, const InteractionOfType
& b
)
844 if (((d
= a
.ai() - b
.ai()) != 0) || ((d
= a
.aj() - b
.aj()) != 0))
850 return a
.interactionTypeName().length() > b
.interactionTypeName().length();
854 static void clean_bonds(InteractionsOfType
* ps
, const gmx::MDLogger
& logger
)
859 for (auto& bond
: ps
->interactionTypes
)
863 std::sort(ps
->interactionTypes
.begin(), ps
->interactionTypes
.end(), pcompar
);
865 /* remove doubles, keep the first one always. */
866 int oldNumber
= ps
->size();
867 for (auto parm
= ps
->interactionTypes
.begin() + 1; parm
!= ps
->interactionTypes
.end();)
869 auto prev
= parm
- 1;
870 if (parm
->ai() == prev
->ai() && parm
->aj() == prev
->aj())
872 parm
= ps
->interactionTypes
.erase(parm
);
881 .appendTextFormatted("Number of bonds was %d, now %zu", oldNumber
, ps
->size());
885 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("No bonds");
889 void print_sums(const t_atoms
* atoms
, bool bSystem
, const gmx::MDLogger
& logger
)
897 where
= " in system";
906 for (i
= 0; (i
< atoms
->nr
); i
++)
908 m
+= atoms
->atom
[i
].m
;
909 qtot
+= atoms
->atom
[i
].q
;
912 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("Total mass%s %.3f a.m.u.", where
, m
);
913 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("Total charge%s %.3f e", where
, qtot
);
916 static void check_restp_type(const char* name
, int t1
, int t2
)
920 gmx_fatal(FARGS
, "Residues in one molecule have a different '%s' type: %d and %d", name
, t1
, t2
);
924 static void check_restp_types(const PreprocessResidue
& r0
, const PreprocessResidue
& r1
)
926 check_restp_type("all dihedrals", static_cast<int>(r0
.bKeepAllGeneratedDihedrals
),
927 static_cast<int>(r1
.bKeepAllGeneratedDihedrals
));
928 check_restp_type("nrexcl", r0
.nrexcl
, r1
.nrexcl
);
929 check_restp_type("HH14", static_cast<int>(r0
.bGenerateHH14Interactions
),
930 static_cast<int>(r1
.bGenerateHH14Interactions
));
931 check_restp_type("remove dihedrals", static_cast<int>(r0
.bRemoveDihedralIfWithImproper
),
932 static_cast<int>(r1
.bRemoveDihedralIfWithImproper
));
934 for (int i
= 0; i
< ebtsNR
; i
++)
936 check_restp_type(btsNames
[i
], r0
.rb
[i
].type
, r1
.rb
[i
].type
);
940 static void add_atom_to_restp(PreprocessResidue
* usedPpResidues
,
943 const MoleculePatch
* patch
)
946 for (int k
= 0; k
< patch
->nr
; k
++)
948 /* set counter in atomname */
949 std::string buf
= patch
->nname
;
952 buf
.append(gmx::formatString("%d", k
+ 1));
954 usedPpResidues
->atomname
.insert(usedPpResidues
->atomname
.begin() + at_start
+ 1 + k
,
955 put_symtab(symtab
, buf
.c_str()));
956 usedPpResidues
->atom
.insert(usedPpResidues
->atom
.begin() + at_start
+ 1 + k
, patch
->atom
.back());
957 if (patch
->cgnr
!= NOTSET
)
959 usedPpResidues
->cgnr
.insert(usedPpResidues
->cgnr
.begin() + at_start
+ 1 + k
, patch
->cgnr
);
963 usedPpResidues
->cgnr
.insert(usedPpResidues
->cgnr
.begin() + at_start
+ 1 + k
,
964 usedPpResidues
->cgnr
[at_start
]);
969 void get_hackblocks_rtp(std::vector
<MoleculePatchDatabase
>* globalPatches
,
970 std::vector
<PreprocessResidue
>* usedPpResidues
,
971 gmx::ArrayRef
<const PreprocessResidue
> rtpFFDB
,
976 gmx::ArrayRef
<MoleculePatchDatabase
*> ntdb
,
977 gmx::ArrayRef
<MoleculePatchDatabase
*> ctdb
,
978 gmx::ArrayRef
<const int> rn
,
979 gmx::ArrayRef
<const int> rc
,
981 const gmx::MDLogger
& logger
)
986 globalPatches
->resize(nres
);
987 usedPpResidues
->clear();
988 /* first the termini */
989 for (int i
= 0; i
< nterpairs
; i
++)
991 if (rn
[i
] >= 0 && ntdb
[i
] != nullptr)
993 copyModificationBlocks(*ntdb
[i
], &globalPatches
->at(rn
[i
]));
995 if (rc
[i
] >= 0 && ctdb
[i
] != nullptr)
997 mergeAtomAndBondModifications(*ctdb
[i
], &globalPatches
->at(rc
[i
]));
1001 /* then the whole rtp */
1002 for (int i
= 0; i
< nres
; i
++)
1004 /* Here we allow a mismatch of one character when looking for the rtp entry.
1005 * For such a mismatch there should be only one mismatching name.
1006 * This is mainly useful for small molecules such as ions.
1007 * Note that this will usually not work for protein, DNA and RNA,
1008 * since there the residue names should be listed in residuetypes.dat
1009 * and an error will have been generated earlier in the process.
1011 key
= *resinfo
[i
].rtp
;
1013 resinfo
[i
].rtp
= put_symtab(symtab
, searchResidueDatabase(key
, rtpFFDB
, logger
).c_str());
1014 auto res
= getDatabaseEntry(*resinfo
[i
].rtp
, rtpFFDB
);
1015 usedPpResidues
->push_back(PreprocessResidue());
1016 PreprocessResidue
* newentry
= &usedPpResidues
->back();
1017 copyPreprocessResidues(*res
, newentry
, symtab
);
1019 /* Check that we do not have different bonded types in one molecule */
1020 check_restp_types(usedPpResidues
->front(), *newentry
);
1023 for (int j
= 0; j
< nterpairs
&& tern
== -1; j
++)
1031 for (int j
= 0; j
< nterpairs
&& terc
== -1; j
++)
1038 bRM
= mergeBondedInteractionList(res
->rb
, globalPatches
->at(i
).rb
, tern
>= 0, terc
>= 0);
1040 if (bRM
&& ((tern
>= 0 && ntdb
[tern
] == nullptr) || (terc
>= 0 && ctdb
[terc
] == nullptr)))
1042 const char* errString
=
1043 "There is a dangling bond at at least one of the terminal ends and the force "
1044 "field does not provide terminal entries or files. Fix your terminal "
1045 "residues so that they match the residue database (.rtp) entries, or provide "
1046 "terminal database entries (.tdb).";
1049 GMX_LOG(logger
.warning
).asParagraph().appendTextFormatted("%s", errString
);
1053 gmx_fatal(FARGS
, "%s", errString
);
1056 else if (bRM
&& ((tern
>= 0 && ntdb
[tern
]->nhack() == 0) || (terc
>= 0 && ctdb
[terc
]->nhack() == 0)))
1058 const char* errString
=
1059 "There is a dangling bond at at least one of the terminal ends. Fix your "
1060 "coordinate file, add a new terminal database entry (.tdb), or select the "
1061 "proper existing terminal entry.";
1064 GMX_LOG(logger
.warning
).asParagraph().appendTextFormatted("%s", errString
);
1068 gmx_fatal(FARGS
, "%s", errString
);
1073 /* Apply patchs to t_restp entries
1074 i.e. add's and deletes from termini database will be
1075 added to/removed from residue topology
1076 we'll do this on one big dirty loop, so it won't make easy reading! */
1077 for (auto modifiedResidue
= globalPatches
->begin(); modifiedResidue
!= globalPatches
->end();
1080 const int pos
= std::distance(globalPatches
->begin(), modifiedResidue
);
1081 PreprocessResidue
* posres
= &usedPpResidues
->at(pos
);
1082 for (auto patch
= modifiedResidue
->hack
.begin(); patch
!= modifiedResidue
->hack
.end(); patch
++)
1086 /* find atom in restp */
1087 auto found
= std::find_if(posres
->atomname
.begin(), posres
->atomname
.end(),
1088 [&patch
](char** name
) {
1089 return (patch
->oname
.empty() && patch
->a
[0] == *name
)
1090 || (patch
->oname
== *name
);
1093 if (found
== posres
->atomname
.end())
1095 /* If we are doing an atom rename only, we don't need
1096 * to generate a fatal error if the old name is not found
1099 /* Deleting can happen also only on the input atoms,
1100 * not necessarily always on the rtp entry.
1102 if (patch
->type() == MoleculePatchType::Add
)
1105 "atom %s not found in buiding block %d%s "
1106 "while combining tdb and rtp",
1107 patch
->oname
.empty() ? patch
->a
[0].c_str() : patch
->oname
.c_str(),
1108 pos
+ 1, *resinfo
[pos
].rtp
);
1113 int l
= std::distance(posres
->atomname
.begin(), found
);
1114 switch (patch
->type())
1116 case MoleculePatchType::Add
:
1119 add_atom_to_restp(posres
, symtab
, l
, &(*patch
));
1122 case MoleculePatchType::Delete
:
1123 { /* we're deleting */
1124 posres
->atom
.erase(posres
->atom
.begin() + l
);
1125 posres
->atomname
.erase(posres
->atomname
.begin() + l
);
1126 posres
->cgnr
.erase(posres
->cgnr
.begin() + l
);
1129 case MoleculePatchType::Replace
:
1131 /* we're replacing */
1132 posres
->atom
[l
] = patch
->atom
.back();
1133 posres
->atomname
[l
] = put_symtab(symtab
, patch
->nname
.c_str());
1134 if (patch
->cgnr
!= NOTSET
)
1136 posres
->cgnr
[l
] = patch
->cgnr
;
1147 static bool atomname_cmp_nr(const char* anm
, const MoleculePatch
* patch
, int* nr
)
1154 return (gmx::equalCaseInsensitive(anm
, patch
->nname
));
1158 if (isdigit(anm
[strlen(anm
) - 1]))
1160 *nr
= anm
[strlen(anm
) - 1] - '0';
1166 if (*nr
<= 0 || *nr
> patch
->nr
)
1172 return (strlen(anm
) == patch
->nname
.length() + 1
1173 && gmx_strncasecmp(anm
, patch
->nname
.c_str(), patch
->nname
.length()) == 0);
1178 static bool match_atomnames_with_rtp_atom(t_atoms
* pdba
,
1179 gmx::ArrayRef
<gmx::RVec
> x
,
1182 PreprocessResidue
* localPpResidue
,
1183 const MoleculePatchDatabase
& singlePatch
,
1185 const gmx::MDLogger
& logger
)
1192 oldnm
= *pdba
->atomname
[atind
];
1193 resnr
= pdba
->resinfo
[pdba
->atom
[atind
].resind
].nr
;
1196 for (auto patch
= singlePatch
.hack
.begin(); patch
!= singlePatch
.hack
.end(); patch
++)
1198 if (patch
->type() == MoleculePatchType::Replace
&& gmx::equalCaseInsensitive(oldnm
, patch
->oname
))
1200 /* This is a replace entry. */
1201 /* Check if we are not replacing a replaced atom. */
1202 bool bReplaceReplace
= false;
1203 for (auto selfPatch
= singlePatch
.hack
.begin(); selfPatch
!= singlePatch
.hack
.end(); selfPatch
++)
1205 if (patch
!= selfPatch
&& selfPatch
->type() == MoleculePatchType::Replace
1206 && gmx::equalCaseInsensitive(selfPatch
->nname
, patch
->oname
))
1208 /* The replace in patch replaces an atom that
1209 * was already replaced in selfPatch, we do not want
1210 * second or higher level replaces at this stage.
1212 bReplaceReplace
= true;
1215 if (bReplaceReplace
)
1217 /* Skip this replace. */
1221 /* This atom still has the old name, rename it */
1222 std::string newnm
= patch
->nname
;
1223 auto found
= std::find_if(
1224 localPpResidue
->atomname
.begin(), localPpResidue
->atomname
.end(),
1225 [&newnm
](char** name
) { return gmx::equalCaseInsensitive(newnm
, *name
); });
1226 if (found
== localPpResidue
->atomname
.end())
1228 /* The new name is not present in the rtp.
1229 * We need to apply the replace also to the rtp entry.
1232 /* We need to find the add hack that can add this atom
1233 * to find out after which atom it should be added.
1235 bool bFoundInAdd
= false;
1236 for (auto rtpModification
= singlePatch
.hack
.begin();
1237 rtpModification
!= singlePatch
.hack
.end(); rtpModification
++)
1239 int k
= std::distance(localPpResidue
->atomname
.begin(), found
);
1240 std::string start_at
;
1241 if (rtpModification
->type() == MoleculePatchType::Add
1242 && atomname_cmp_nr(newnm
.c_str(), &(*rtpModification
), &anmnr
))
1246 start_at
= singlePatch
.hack
[k
].a
[0];
1250 start_at
= gmx::formatString("%s%d", singlePatch
.hack
[k
].nname
.c_str(),
1254 std::find_if(localPpResidue
->atomname
.begin(),
1255 localPpResidue
->atomname
.end(), [&start_at
](char** name
) {
1256 return gmx::equalCaseInsensitive(start_at
, *name
);
1258 if (found2
== localPpResidue
->atomname
.end())
1261 "Could not find atom '%s' in residue building block '%s' to "
1263 start_at
.c_str(), localPpResidue
->resname
.c_str(), newnm
.c_str());
1265 /* We can add the atom after atom start_nr */
1266 add_atom_to_restp(localPpResidue
, symtab
,
1267 std::distance(localPpResidue
->atomname
.begin(), found2
),
1277 "Could not find an 'add' entry for atom named '%s' corresponding to "
1278 "the 'replace' entry from atom name '%s' to '%s' for tdb or hdb "
1279 "database of residue type '%s'",
1280 newnm
.c_str(), patch
->oname
.c_str(), patch
->nname
.c_str(),
1281 localPpResidue
->resname
.c_str());
1287 GMX_LOG(logger
.info
)
1289 .appendTextFormatted("Renaming atom '%s' in residue '%s' %d to '%s'", oldnm
,
1290 localPpResidue
->resname
.c_str(), resnr
, newnm
.c_str());
1292 /* Rename the atom in pdba */
1293 pdba
->atomname
[atind
] = put_symtab(symtab
, newnm
.c_str());
1295 else if (patch
->type() == MoleculePatchType::Delete
1296 && gmx::equalCaseInsensitive(oldnm
, patch
->oname
))
1298 /* This is a delete entry, check if this atom is present
1299 * in the rtp entry of this residue.
1301 auto found3
= std::find_if(
1302 localPpResidue
->atomname
.begin(), localPpResidue
->atomname
.end(),
1303 [&oldnm
](char** name
) { return gmx::equalCaseInsensitive(oldnm
, *name
); });
1304 if (found3
== localPpResidue
->atomname
.end())
1306 /* This atom is not present in the rtp entry,
1307 * delete is now from the input pdba.
1311 GMX_LOG(logger
.info
)
1313 .appendTextFormatted("Deleting atom '%s' in residue '%s' %d", oldnm
,
1314 localPpResidue
->resname
.c_str(), resnr
);
1316 /* We should free the atom name,
1317 * but it might be used multiple times in the symtab.
1318 * sfree(pdba->atomname[atind]);
1320 for (int k
= atind
+ 1; k
< pdba
->nr
; k
++)
1322 pdba
->atom
[k
- 1] = pdba
->atom
[k
];
1323 pdba
->atomname
[k
- 1] = pdba
->atomname
[k
];
1324 copy_rvec(x
[k
], x
[k
- 1]);
1335 void match_atomnames_with_rtp(gmx::ArrayRef
<PreprocessResidue
> usedPpResidues
,
1336 gmx::ArrayRef
<MoleculePatchDatabase
> globalPatches
,
1339 gmx::ArrayRef
<gmx::RVec
> x
,
1341 const gmx::MDLogger
& logger
)
1343 for (int i
= 0; i
< pdba
->nr
; i
++)
1345 const char* oldnm
= *pdba
->atomname
[i
];
1346 PreprocessResidue
* localPpResidue
= &usedPpResidues
[pdba
->atom
[i
].resind
];
1347 auto found
= std::find_if(
1348 localPpResidue
->atomname
.begin(), localPpResidue
->atomname
.end(),
1349 [&oldnm
](char** name
) { return gmx::equalCaseInsensitive(oldnm
, *name
); });
1350 if (found
== localPpResidue
->atomname
.end())
1352 /* Not found yet, check if we have to rename this atom */
1353 if (match_atomnames_with_rtp_atom(pdba
, x
, symtab
, i
, localPpResidue
,
1354 globalPatches
[pdba
->atom
[i
].resind
], bVerbose
, logger
))
1356 /* We deleted this atom, decrease the atom counter by 1. */
1363 #define NUM_CMAP_ATOMS 5
1364 static void gen_cmap(InteractionsOfType
* psb
,
1365 gmx::ArrayRef
<const PreprocessResidue
> usedPpResidues
,
1367 gmx::ArrayRef
<const int> cyclicBondsIndex
,
1368 const gmx::MDLogger
& logger
)
1372 t_resinfo
* resinfo
= atoms
->resinfo
;
1373 int nres
= atoms
->nres
;
1374 int cmap_atomid
[NUM_CMAP_ATOMS
];
1375 int cmap_chainnum
= -1;
1386 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("Making cmap torsions...");
1388 /* Most cmap entries use the N atom from the next residue, so the last
1389 * residue should not have its CMAP entry in that case, but for things like
1390 * dipeptides we sometimes define a complete CMAP entry inside a residue,
1391 * and in this case we need to process everything through the last residue.
1393 for (residx
= 0; residx
< nres
; residx
++)
1395 /* Add CMAP terms from the list of CMAP interactions */
1396 for (const auto& b
: usedPpResidues
[residx
].rb
[ebtsCMAP
].b
)
1398 bool bAddCMAP
= true;
1399 /* Loop over atoms in a candidate CMAP interaction and
1400 * check that they exist, are from the same chain and are
1401 * from residues labelled as protein. */
1402 for (int k
= 0; k
< NUM_CMAP_ATOMS
&& bAddCMAP
; k
++)
1404 /* Assign the pointer to the name of the next reference atom.
1405 * This can use -/+ labels to refer to previous/next residue.
1407 const char* pname
= b
.a
[k
].c_str();
1408 /* Skip this CMAP entry if it refers to residues before the
1409 * first or after the last residue.
1411 if ((cyclicBondsIndex
.empty() && ((strchr(pname
, '-') != nullptr) && (residx
== 0)))
1412 || ((strchr(pname
, '+') != nullptr) && (residx
== nres
- 1)))
1418 cmap_atomid
[k
] = search_atom(pname
, i
, atoms
, ptr
, TRUE
, cyclicBondsIndex
);
1419 bAddCMAP
= bAddCMAP
&& (cmap_atomid
[k
] != -1);
1422 /* This break is necessary, because cmap_atomid[k]
1423 * == -1 cannot be safely used as an index
1424 * into the atom array. */
1427 int this_residue_index
= atoms
->atom
[cmap_atomid
[k
]].resind
;
1430 cmap_chainnum
= resinfo
[this_residue_index
].chainnum
;
1434 /* Does the residue for this atom have the same
1435 * chain number as the residues for previous
1437 bAddCMAP
= bAddCMAP
&& cmap_chainnum
== resinfo
[this_residue_index
].chainnum
;
1439 /* Here we used to check that the residuetype was protein and
1440 * disable bAddCMAP if that was not the case. However, some
1441 * special residues (say, alanine dipeptides) might not adhere
1442 * to standard naming, and if we start calling them normal
1443 * protein residues the user will be bugged to select termini.
1445 * Instead, I believe that the right course of action is to
1446 * keep the CMAP interaction if it is present in the RTP file
1447 * and we correctly identified all atoms (which is the case
1454 add_cmap_param(psb
, cmap_atomid
[0], cmap_atomid
[1], cmap_atomid
[2], cmap_atomid
[3],
1455 cmap_atomid
[4], b
.s
.c_str());
1459 if (residx
< nres
- 1)
1461 while (atoms
->atom
[i
].resind
< residx
+ 1)
1467 /* Start the next residue */
1470 static void scrub_charge_groups(int* cgnr
, int natoms
)
1474 for (i
= 0; i
< natoms
; i
++)
1481 void pdb2top(FILE* top_file
,
1482 const char* posre_fn
,
1483 const char* molname
,
1485 std::vector
<gmx::RVec
>* x
,
1486 PreprocessingAtomTypes
* atype
,
1488 gmx::ArrayRef
<const PreprocessResidue
> rtpFFDB
,
1489 gmx::ArrayRef
<PreprocessResidue
> usedPpResidues
,
1490 gmx::ArrayRef
<MoleculePatchDatabase
> globalPatches
,
1493 bool bVsiteAromatics
,
1496 gmx::ArrayRef
<const DisulfideBond
> ssbonds
,
1497 real long_bond_dist
,
1498 real short_bond_dist
,
1504 gmx::ArrayRef
<const int> cyclicBondsIndex
,
1505 const gmx::MDLogger
& logger
)
1507 std::array
<InteractionsOfType
, F_NRE
> plist
;
1517 at2bonds(&(plist
[F_BONDS
]), globalPatches
, atoms
, *x
, long_bond_dist
, short_bond_dist
,
1518 cyclicBondsIndex
, logger
);
1520 /* specbonds: disulphide bonds & heme-his */
1521 do_ssbonds(&(plist
[F_BONDS
]), atoms
, ssbonds
, bAllowMissing
);
1523 nmissat
= name2type(atoms
, &cgnr
, usedPpResidues
, &rt
, logger
);
1528 GMX_LOG(logger
.warning
)
1530 .appendTextFormatted("There were %d missing atoms in molecule %s", nmissat
, molname
);
1535 "There were %d missing atoms in molecule %s, if you want to use this "
1536 "incomplete topology anyhow, use the option -missing",
1541 /* Cleanup bonds (sort and rm doubles) */
1542 clean_bonds(&(plist
[F_BONDS
]), logger
);
1544 snew(vsite_type
, atoms
->nr
);
1545 for (i
= 0; i
< atoms
->nr
; i
++)
1547 vsite_type
[i
] = NOTSET
;
1551 if (bVsiteAromatics
)
1553 GMX_LOG(logger
.info
)
1555 .appendTextFormatted(
1556 "The conversion of aromatic rings into virtual sites is deprecated "
1557 "and may be removed in a future version of GROMACS");
1559 /* determine which atoms will be vsites and add dummy masses
1560 also renumber atom numbers in plist[0..F_NRE]! */
1561 do_vsites(rtpFFDB
, atype
, atoms
, tab
, x
, plist
, &vsite_type
, &cgnr
, mHmult
, bVsiteAromatics
, ffdir
);
1564 /* Make Angles and Dihedrals */
1565 GMX_LOG(logger
.info
)
1567 .appendTextFormatted("Generating angles, dihedrals and pairs...");
1568 snew(excls
, atoms
->nr
);
1569 gen_pad(atoms
, usedPpResidues
, plist
, excls
, globalPatches
, bAllowMissing
, cyclicBondsIndex
);
1574 gen_cmap(&(plist
[F_CMAP
]), usedPpResidues
, atoms
, cyclicBondsIndex
, logger
);
1575 if (plist
[F_CMAP
].size() > 0)
1577 GMX_LOG(logger
.info
)
1579 .appendTextFormatted("There are %4zu cmap torsion pairs", plist
[F_CMAP
].size());
1583 /* set mass of all remaining hydrogen atoms */
1586 do_h_mass(&(plist
[F_BONDS
]), vsite_type
, atoms
, mHmult
, bDeuterate
);
1590 /* Cleanup bonds (sort and rm doubles) */
1591 /* clean_bonds(&(plist[F_BONDS]));*/
1593 GMX_LOG(logger
.info
)
1595 .appendTextFormatted(
1596 "There are %4zu dihedrals, %4zu impropers, %4zu angles\n"
1597 " %4zu pairs, %4zu bonds and %4zu virtual sites",
1598 plist
[F_PDIHS
].size(), plist
[F_IDIHS
].size(), plist
[F_ANGLES
].size(),
1599 plist
[F_LJ14
].size(), plist
[F_BONDS
].size(),
1600 plist
[F_VSITE2
].size() + plist
[F_VSITE3
].size() + plist
[F_VSITE3FD
].size()
1601 + plist
[F_VSITE3FAD
].size() + plist
[F_VSITE3OUT
].size()
1602 + plist
[F_VSITE4FD
].size() + plist
[F_VSITE4FDN
].size());
1604 print_sums(atoms
, FALSE
, logger
);
1608 scrub_charge_groups(cgnr
, atoms
->nr
);
1613 for (i
= 0; i
< atoms
->nres
; i
++)
1615 atoms
->resinfo
[i
].nr
= i
+ 1;
1616 atoms
->resinfo
[i
].ic
= ' ';
1622 GMX_LOG(logger
.info
).asParagraph().appendTextFormatted("Writing topology");
1623 /* We can copy the bonded types from the first restp,
1624 * since the types have to be identical for all residues in one molecule.
1626 for (i
= 0; i
< ebtsNR
; i
++)
1628 bts
[i
] = usedPpResidues
[0].rb
[i
].type
;
1630 write_top(top_file
, posre_fn
, molname
, atoms
, bRTPresname
, bts
, plist
, excls
, atype
, cgnr
,
1631 usedPpResidues
[0].nrexcl
);
1635 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1637 for (i
= 0; i
< atoms
->nr
; i
++)