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,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * 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/gpp_nextnb.h"
57 #include "gromacs/gmxpreprocess/h_db.h"
58 #include "gromacs/gmxpreprocess/notset.h"
59 #include "gromacs/gmxpreprocess/pgutil.h"
60 #include "gromacs/gmxpreprocess/resall.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/strdb.h"
79 #include "gromacs/utility/stringutil.h"
80 #include "gromacs/utility/textwriter.h"
82 /* this must correspond to enum in pdb2top.h */
83 const char *hh[ehisNR] = { "HISD", "HISE", "HISH", "HIS1" };
85 static int missing_atoms(t_restp *rp, int resind, t_atoms *at, int i0, int i)
92 for (j = 0; j < rp->natom; j++)
94 name = *(rp->atomname[j]);
96 for (k = i0; k < i; k++)
98 bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]), name));
103 fprintf(stderr, "\nWARNING: "
104 "atom %s is missing in residue %s %d in the pdb file\n",
105 name, *(at->resinfo[resind].name), at->resinfo[resind].nr);
106 if (name[0] == 'H' || name[0] == 'h')
108 fprintf(stderr, " You might need to add atom %s to the hydrogen database of building block %s\n"
109 " in the file %s.hdb (see the manual)\n",
110 name, *(at->resinfo[resind].rtp), rp->filebase);
112 fprintf(stderr, "\n");
119 gmx_bool is_int(double x)
121 const double tol = 1e-4;
130 return (fabs(x-ix) < tol);
134 choose_ff_impl(const char *ffsel,
135 char *forcefield, int ff_maxlen,
136 char *ffdir, int ffdir_maxlen)
138 std::vector<gmx::DataFileInfo> ffdirs = fflib_enumerate_forcefields();
139 const int nff = static_cast<int>(ffdirs.size());
141 /* Replace with unix path separators */
142 if (DIR_SEPARATOR != '/')
144 for (int i = 0; i < nff; ++i)
146 std::replace(ffdirs[i].dir.begin(), ffdirs[i].dir.end(), DIR_SEPARATOR, '/');
150 /* Store the force field names in ffs */
151 std::vector<std::string> ffs;
152 ffs.reserve(ffdirs.size());
153 for (int i = 0; i < nff; ++i)
155 ffs.push_back(gmx::stripSuffixIfPresent(ffdirs[i].name,
156 fflib_forcefield_dir_ext()));
160 if (ffsel != nullptr)
165 for (int i = 0; i < nff; ++i)
169 /* Matching ff name */
173 if (ffdirs[i].dir == ".")
190 "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
191 "current directory. Use interactive selection (not the -ff option) if\n"
192 "you would prefer a different one.\n", ffsel, nfound);
196 std::string message = gmx::formatString(
197 "Force field '%s' occurs in %d places, but not in "
198 "the current directory.\n"
199 "Run without the -ff switch and select the force "
200 "field interactively.", ffsel, nfound);
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.", ffsel);
209 GMX_THROW(gmx::InconsistentInputError(message));
214 std::vector<std::string> desc;
215 desc.reserve(ffdirs.size());
216 for (int i = 0; i < nff; ++i)
218 std::string docFileName(
219 gmx::Path::join(ffdirs[i].dir, ffdirs[i].name,
220 fflib_forcefield_doc()));
221 // TODO: Just try to open the file with a method that does not
222 // throw/bail out with a fatal error instead of multiple checks.
223 if (gmx::File::exists(docFileName, gmx::File::returnFalseOnError))
225 // TODO: Use a C++ API without such an intermediate/fixed-length buffer.
227 /* We don't use fflib_open, because we don't want printf's */
228 FILE *fp = gmx_ffopen(docFileName.c_str(), "r");
229 get_a_line(fp, buf, STRLEN);
231 desc.emplace_back(buf);
235 desc.push_back(ffs[i]);
238 /* Order force fields from the same dir alphabetically
239 * and put deprecated force fields at the end.
241 for (int i = 0; i < nff; ++i)
243 for (int j = i + 1; j < nff; ++j)
245 if (ffdirs[i].dir == ffdirs[j].dir &&
246 ((desc[i][0] == '[' && desc[j][0] != '[') ||
247 ((desc[i][0] == '[' || desc[j][0] != '[') &&
248 gmx_strcasecmp(desc[i].c_str(), desc[j].c_str()) > 0)))
250 std::swap(ffdirs[i].name, ffdirs[j].name);
251 std::swap(ffs[i], ffs[j]);
252 std::swap(desc[i], desc[j]);
257 printf("\nSelect the Force Field:\n");
258 for (int i = 0; i < nff; ++i)
260 if (i == 0 || ffdirs[i-1].dir != ffdirs[i].dir)
262 if (ffdirs[i].dir == ".")
264 printf("From current directory:\n");
268 printf("From '%s':\n", ffdirs[i].dir.c_str());
271 printf("%2d: %s\n", i+1, desc[i].c_str());
275 // TODO: Add a C++ API for this.
280 pret = fgets(buf, STRLEN, stdin);
284 sel = strtol(buf, nullptr, 10);
288 while (pret == nullptr || (sel < 0) || (sel >= nff));
290 /* Check for a current limitation of the fflib code.
291 * It will always read from the first ff directory in the list.
292 * This check assumes that the order of ffs matches the order
293 * in which fflib_open searches ff library files.
295 for (int i = 0; i < sel; i++)
297 if (ffs[i] == ffs[sel])
299 std::string message = gmx::formatString(
300 "Can only select the first of multiple force "
301 "field entries with directory name '%s%s' in "
302 "the list. If you want to use the next entry, "
303 "run pdb2gmx in a different directory, set GMXLIB "
304 "to point to the desired force field first, and/or "
305 "rename or move the force field directory present "
306 "in the current working directory.",
307 ffs[sel].c_str(), fflib_forcefield_dir_ext());
308 GMX_THROW(gmx::NotImplementedError(message));
317 if (ffs[sel].length() >= static_cast<size_t>(ff_maxlen))
319 std::string message = gmx::formatString(
320 "Length of force field name (%d) >= maxlen (%d)",
321 static_cast<int>(ffs[sel].length()), ff_maxlen);
322 GMX_THROW(gmx::InvalidInputError(message));
324 strcpy(forcefield, ffs[sel].c_str());
327 if (ffdirs[sel].bFromDefaultDir)
329 ffpath = ffdirs[sel].name;
333 ffpath = gmx::Path::join(ffdirs[sel].dir, ffdirs[sel].name);
335 if (ffpath.length() >= static_cast<size_t>(ffdir_maxlen))
337 std::string message = gmx::formatString(
338 "Length of force field dir (%d) >= maxlen (%d)",
339 static_cast<int>(ffpath.length()), ffdir_maxlen);
340 GMX_THROW(gmx::InvalidInputError(message));
342 strcpy(ffdir, ffpath.c_str());
346 choose_ff(const char *ffsel,
347 char *forcefield, int ff_maxlen,
348 char *ffdir, int ffdir_maxlen)
352 choose_ff_impl(ffsel, forcefield, ff_maxlen, ffdir, ffdir_maxlen);
354 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
357 void choose_watermodel(const char *wmsel, const char *ffdir,
360 const char *fn_watermodels = "watermodels.dat";
361 char fn_list[STRLEN];
368 if (strcmp(wmsel, "none") == 0)
370 *watermodel = nullptr;
374 else if (strcmp(wmsel, "select") != 0)
376 *watermodel = gmx_strdup(wmsel);
381 sprintf(fn_list, "%s%c%s", ffdir, DIR_SEPARATOR, fn_watermodels);
383 if (!fflib_fexist(fn_list))
385 fprintf(stderr, "No file '%s' found, will not include a water model\n",
387 *watermodel = nullptr;
392 fp = fflib_open(fn_list);
393 printf("\nSelect the Water Model:\n");
396 while (get_a_line(fp, buf, STRLEN))
398 srenew(model, nwm+1);
399 snew(model[nwm], STRLEN);
400 sscanf(buf, "%s%n", model[nwm], &i);
404 fprintf(stderr, "%2d: %s\n", nwm+1, buf+i);
413 fprintf(stderr, "%2d: %s\n", nwm+1, "None");
418 pret = fgets(buf, STRLEN, stdin);
422 sel = strtol(buf, nullptr, 10);
426 while (pret == nullptr || sel < 0 || sel > nwm);
430 *watermodel = nullptr;
434 *watermodel = gmx_strdup(model[sel]);
437 for (i = 0; i < nwm; i++)
444 static int name2type(t_atoms *at, int **cgnr,
445 t_restp restp[], gmx_residuetype_t *rt)
447 int i, j, prevresind, resind, i0, prevcg, cg, curcg;
465 for (i = 0; (i < at->nr); i++)
468 if (at->atom[i].resind != resind)
471 resind = at->atom[i].resind;
472 bProt = gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name));
473 bNterm = bProt && (resind == 0);
476 nmissat += missing_atoms(&restp[prevresind], prevresind, at, i0, i);
480 if (at->atom[i].m == 0)
484 fprintf(debug, "atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
485 i+1, *(at->atomname[i]), curcg, prevcg,
486 j == NOTSET ? NOTSET : restp[resind].cgnr[j]);
490 name = *(at->atomname[i]);
491 j = search_jtype(&restp[resind], name, bNterm);
492 at->atom[i].type = restp[resind].atom[j].type;
493 at->atom[i].q = restp[resind].atom[j].q;
494 at->atom[i].m = restp[resind].atom[j].m;
495 cg = restp[resind].cgnr[j];
496 /* A charge group number -1 signals a separate charge group
499 if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) )
508 fprintf(debug, "atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
509 i+1, *(at->atomname[i]), curcg, qt, is_int(qt));
520 at->atom[i].typeB = at->atom[i].type;
521 at->atom[i].qB = at->atom[i].q;
522 at->atom[i].mB = at->atom[i].m;
524 nmissat += missing_atoms(&restp[resind], resind, at, i0, i);
529 static void print_top_heavy_H(FILE *out, real mHmult)
533 fprintf(out, "; Using deuterium instead of hydrogen\n\n");
535 else if (mHmult == 4.0)
537 fprintf(out, "#define HEAVY_H\n\n");
539 else if (mHmult != 1.0)
541 fprintf(stderr, "WARNING: unsupported proton mass multiplier (%g) "
542 "in pdb2top\n", mHmult);
546 void print_top_comment(FILE *out,
547 const char *filename,
551 char ffdir_parent[STRLEN];
556 gmx::TextWriter writer(out);
557 gmx::niceHeader(&writer, filename, ';');
558 writer.writeLine(gmx::formatString(";\tThis is a %s topology file", bITP ? "include" : "standalone"));
559 writer.writeLine(";");
561 gmx::BinaryInformationSettings settings;
562 settings.generatedByHeader(true);
563 settings.linePrefix(";\t");
564 gmx::printBinaryInformation(&writer, gmx::getProgramContext(), settings);
566 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
568 if (strchr(ffdir, '/') == nullptr)
570 fprintf(out, ";\tForce field was read from the standard GROMACS share directory.\n;\n\n");
572 else if (ffdir[0] == '.')
574 fprintf(out, ";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
578 strncpy(ffdir_parent, ffdir, STRLEN-1);
579 ffdir_parent[STRLEN-1] = '\0'; /*make sure it is 0-terminated even for long string*/
580 p = strrchr(ffdir_parent, '/');
585 ";\tForce field data was read from:\n"
589 ";\tThis might be a non-standard force field location. When you use this topology, the\n"
590 ";\tforce field must either be present in the current directory, or the location\n"
591 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
596 void print_top_header(FILE *out, const char *filename,
597 gmx_bool bITP, const char *ffdir, real mHmult)
601 print_top_comment(out, filename, ffdir, bITP);
603 print_top_heavy_H(out, mHmult);
604 fprintf(out, "; Include forcefield parameters\n");
606 p = strrchr(ffdir, '/');
607 p = (ffdir[0] == '.' || p == nullptr) ? ffdir : p+1;
609 fprintf(out, "#include \"%s/%s\"\n\n", p, fflib_forcefield_itp());
612 static void print_top_posre(FILE *out, const char *pr)
614 fprintf(out, "; Include Position restraint file\n");
615 fprintf(out, "#ifdef POSRES\n");
616 fprintf(out, "#include \"%s\"\n", pr);
617 fprintf(out, "#endif\n\n");
620 static void print_top_water(FILE *out, const char *ffdir, const char *water)
625 fprintf(out, "; Include water topology\n");
627 p = strrchr(ffdir, '/');
628 p = (ffdir[0] == '.' || p == nullptr) ? ffdir : p+1;
629 fprintf(out, "#include \"%s/%s.itp\"\n", p, water);
632 fprintf(out, "#ifdef POSRES_WATER\n");
633 fprintf(out, "; Position restraint for each water oxygen\n");
634 fprintf(out, "[ position_restraints ]\n");
635 fprintf(out, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
636 fprintf(out, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
637 fprintf(out, "#endif\n");
640 sprintf(buf, "%s/ions.itp", p);
642 if (fflib_fexist(buf))
644 fprintf(out, "; Include topology for ions\n");
645 fprintf(out, "#include \"%s\"\n", buf);
650 static void print_top_system(FILE *out, const char *title)
652 fprintf(out, "[ %s ]\n", dir2str(d_system));
653 fprintf(out, "; Name\n");
654 fprintf(out, "%s\n\n", title[0] ? title : "Protein");
657 void print_top_mols(FILE *out,
658 const char *title, const char *ffdir, const char *water,
659 int nincl, char **incls, int nmol, t_mols *mols)
666 fprintf(out, "; Include chain topologies\n");
667 for (i = 0; (i < nincl); i++)
669 incl = strrchr(incls[i], DIR_SEPARATOR);
676 /* Remove the path from the include name */
679 fprintf(out, "#include \"%s\"\n", incl);
686 print_top_water(out, ffdir, water);
688 print_top_system(out, title);
692 fprintf(out, "[ %s ]\n", dir2str(d_molecules));
693 fprintf(out, "; %-15s %5s\n", "Compound", "#mols");
694 for (i = 0; (i < nmol); i++)
696 fprintf(out, "%-15s %5d\n", mols[i].name, mols[i].nr);
701 void write_top(FILE *out, char *pr, const char *molname,
702 t_atoms *at, gmx_bool bRTPresname,
703 int bts[], t_params plist[], t_excls excls[],
704 gpp_atomtype_t atype, int *cgnr, int nrexcl)
705 /* NOTE: nrexcl is not the size of *excl! */
707 if (at && atype && cgnr)
709 fprintf(out, "[ %s ]\n", dir2str(d_moleculetype));
710 fprintf(out, "; %-15s %5s\n", "Name", "nrexcl");
711 fprintf(out, "%-15s %5d\n\n", molname ? molname : "Protein", nrexcl);
713 print_atoms(out, atype, at, cgnr, bRTPresname);
714 print_bondeds(out, at->nr, d_bonds, F_BONDS, bts[ebtsBONDS], plist);
715 print_bondeds(out, at->nr, d_constraints, F_CONSTR, 0, plist);
716 print_bondeds(out, at->nr, d_constraints, F_CONSTRNC, 0, plist);
717 print_bondeds(out, at->nr, d_pairs, F_LJ14, 0, plist);
718 print_excl(out, at->nr, excls);
719 print_bondeds(out, at->nr, d_angles, F_ANGLES, bts[ebtsANGLES], plist);
720 print_bondeds(out, at->nr, d_dihedrals, F_PDIHS, bts[ebtsPDIHS], plist);
721 print_bondeds(out, at->nr, d_dihedrals, F_IDIHS, bts[ebtsIDIHS], plist);
722 print_bondeds(out, at->nr, d_cmap, F_CMAP, bts[ebtsCMAP], plist);
723 print_bondeds(out, at->nr, d_polarization, F_POLARIZATION, 0, plist);
724 print_bondeds(out, at->nr, d_thole_polarization, F_THOLE_POL, 0, plist);
725 print_bondeds(out, at->nr, d_vsites2, F_VSITE2, 0, plist);
726 print_bondeds(out, at->nr, d_vsites3, F_VSITE3, 0, plist);
727 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FD, 0, plist);
728 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FAD, 0, plist);
729 print_bondeds(out, at->nr, d_vsites3, F_VSITE3OUT, 0, plist);
730 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FD, 0, plist);
731 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FDN, 0, plist);
735 print_top_posre(out, pr);
742 static void do_ssbonds(t_params *ps, t_atoms *atoms,
743 int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
748 for (i = 0; (i < nssbonds); i++)
750 ri = ssbonds[i].res1;
751 rj = ssbonds[i].res2;
752 ai = search_res_atom(ssbonds[i].a1, ri, atoms,
753 "special bond", bAllowMissing);
754 aj = search_res_atom(ssbonds[i].a2, rj, atoms,
755 "special bond", bAllowMissing);
756 if ((ai == -1) || (aj == -1))
758 gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
759 ssbonds[i].a1, ssbonds[i].a2);
761 add_param(ps, ai, aj, nullptr, nullptr);
765 static void at2bonds(t_params *psb, t_hackblock *hb,
768 real long_bond_dist, real short_bond_dist)
772 real dist2, long_bond_dist2, short_bond_dist2;
775 long_bond_dist2 = gmx::square(long_bond_dist);
776 short_bond_dist2 = gmx::square(short_bond_dist);
787 fprintf(stderr, "Making bonds...\n");
789 for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
791 /* add bonds from list of bonded interactions */
792 for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
794 /* Unfortunately we can not issue errors or warnings
795 * for missing atoms in bonds, as the hydrogens and terminal atoms
796 * have not been added yet.
798 ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[0], i, atoms,
800 aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[1], i, atoms,
802 if (ai != -1 && aj != -1)
804 dist2 = distance2(x[ai], x[aj]);
805 if (dist2 > long_bond_dist2)
808 fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
809 ai+1, aj+1, std::sqrt(dist2));
811 else if (dist2 < short_bond_dist2)
813 fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
814 ai+1, aj+1, std::sqrt(dist2));
816 add_param(psb, ai, aj, nullptr, hb[resind].rb[ebtsBONDS].b[j].s);
819 /* add bonds from list of hacks (each added atom gets a bond) */
820 while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
822 for (j = 0; j < hb[resind].nhack; j++)
824 if ( ( hb[resind].hack[j].tp > 0 ||
825 hb[resind].hack[j].oname == nullptr ) &&
826 strcmp(hb[resind].hack[j].a[0], *(atoms->atomname[i])) == 0)
828 switch (hb[resind].hack[j].tp)
830 case 9: /* COOH terminus */
831 add_param(psb, i, i+1, nullptr, nullptr); /* C-O */
832 add_param(psb, i, i+2, nullptr, nullptr); /* C-OA */
833 add_param(psb, i+2, i+3, nullptr, nullptr); /* OA-H */
836 for (k = 0; (k < hb[resind].hack[j].nr); k++)
838 add_param(psb, i, i+k+1, nullptr, nullptr);
845 /* we're now at the start of the next residue */
849 static int pcompar(const void *a, const void *b)
851 const t_param *pa, *pb;
853 pa = static_cast<const t_param *>(a);
854 pb = static_cast<const t_param *>(b);
856 d = pa->a[0] - pb->a[0];
859 d = pa->a[1] - pb->a[1];
863 return strlen(pb->s) - strlen(pa->s);
871 static void clean_bonds(t_params *ps)
878 /* swap atomnumbers in bond if first larger than second: */
879 for (i = 0; (i < ps->nr); i++)
881 if (ps->param[i].a[1] < ps->param[i].a[0])
883 a = ps->param[i].a[0];
884 ps->param[i].a[0] = ps->param[i].a[1];
885 ps->param[i].a[1] = a;
890 qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
892 /* remove doubles, keep the first one always. */
894 for (i = 1; (i < ps->nr); i++)
896 if ((ps->param[i].a[0] != ps->param[j-1].a[0]) ||
897 (ps->param[i].a[1] != ps->param[j-1].a[1]) )
901 cp_param(&(ps->param[j]), &(ps->param[i]));
906 fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
911 fprintf(stderr, "No bonds\n");
915 void print_sums(t_atoms *atoms, gmx_bool bSystem)
923 where = " in system";
932 for (i = 0; (i < atoms->nr); i++)
934 m += atoms->atom[i].m;
935 qtot += atoms->atom[i].q;
938 fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
939 fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
942 static void check_restp_type(const char *name, int t1, int t2)
946 gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
950 static void check_restp_types(t_restp *r0, t_restp *r1)
954 check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
955 check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
956 check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
957 check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
959 for (i = 0; i < ebtsNR; i++)
961 check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
965 static void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
969 const char *Hnum = "123456";
973 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
975 * restp->atomname[at_start], resnr, restp->resname);
977 strcpy(buf, hack->nname);
978 buf[strlen(buf)+1] = '\0';
981 buf[strlen(buf)] = '-';
984 restp->natom += hack->nr;
985 srenew(restp->atom, restp->natom);
986 srenew(restp->atomname, restp->natom);
987 srenew(restp->cgnr, restp->natom);
989 for (k = restp->natom-1; k > at_start+hack->nr; k--)
992 restp->atom [k - hack->nr];
994 restp->atomname[k - hack->nr];
996 restp->cgnr [k - hack->nr];
999 for (k = 0; k < hack->nr; k++)
1001 /* set counter in atomname */
1004 buf[strlen(buf)-1] = Hnum[k];
1006 snew( restp->atomname[at_start+1+k], 1);
1007 restp->atom [at_start+1+k] = *hack->atom;
1008 *restp->atomname[at_start+1+k] = gmx_strdup(buf);
1009 if (hack->cgnr != NOTSET)
1011 restp->cgnr [at_start+1+k] = hack->cgnr;
1015 restp->cgnr [at_start+1+k] = restp->cgnr[at_start];
1020 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1021 int nrtp, t_restp rtp[],
1022 int nres, t_resinfo *resinfo,
1024 t_hackblock **ntdb, t_hackblock **ctdb,
1025 const int *rn, const int *rc,
1026 gmx_bool bAllowMissing)
1036 /* first the termini */
1037 for (i = 0; i < nterpairs; i++)
1039 if (rn[i] >= 0 && ntdb[i] != nullptr)
1041 copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1043 if (rc[i] >= 0 && ctdb[i] != nullptr)
1045 merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1049 /* then the whole rtp */
1050 for (i = 0; i < nres; i++)
1052 /* Here we allow a mismatch of one character when looking for the rtp entry.
1053 * For such a mismatch there should be only one mismatching name.
1054 * This is mainly useful for small molecules such as ions.
1055 * Note that this will usually not work for protein, DNA and RNA,
1056 * since there the residue names should be listed in residuetypes.dat
1057 * and an error will have been generated earlier in the process.
1059 key = *resinfo[i].rtp;
1060 snew(resinfo[i].rtp, 1);
1061 *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1062 res = get_restp(*resinfo[i].rtp, nrtp, rtp);
1063 copy_t_restp(res, &(*restp)[i]);
1065 /* Check that we do not have different bonded types in one molecule */
1066 check_restp_types(&(*restp)[0], &(*restp)[i]);
1069 for (j = 0; j < nterpairs && tern == -1; j++)
1077 for (j = 0; j < nterpairs && terc == -1; j++)
1084 bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1086 if (bRM && ((tern >= 0 && ntdb[tern] == nullptr) ||
1087 (terc >= 0 && ctdb[terc] == nullptr)))
1089 const char *errString = "There is a dangling bond at at least one of the terminal ends and the force field does not provide terminal entries or files. Fix your terminal residues so that they match the residue database (.rtp) entries, or provide terminal database entries (.tdb).";
1092 fprintf(stderr, "%s\n", errString);
1096 gmx_fatal(FARGS, errString);
1099 else if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1100 (terc >= 0 && ctdb[terc]->nhack == 0)))
1102 const char *errString = "There is a dangling bond at at least one of the terminal ends. Fix your coordinate file, add a new terminal database entry (.tdb), or select the proper existing terminal entry.";
1105 fprintf(stderr, "%s\n", errString);
1109 gmx_fatal(FARGS, errString);
1114 /* now perform t_hack's on t_restp's,
1115 i.e. add's and deletes from termini database will be
1116 added to/removed from residue topology
1117 we'll do this on one big dirty loop, so it won't make easy reading! */
1118 for (i = 0; i < nres; i++)
1120 for (j = 0; j < (*hb)[i].nhack; j++)
1122 if ( (*hb)[i].hack[j].nr)
1124 /* find atom in restp */
1125 for (l = 0; l < (*restp)[i].natom; l++)
1127 if ( ( (*hb)[i].hack[j].oname == nullptr &&
1128 strcmp((*hb)[i].hack[j].a[0], *(*restp)[i].atomname[l]) == 0 ) ||
1129 ( (*hb)[i].hack[j].oname != nullptr &&
1130 strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1135 if (l == (*restp)[i].natom)
1137 /* If we are doing an atom rename only, we don't need
1138 * to generate a fatal error if the old name is not found
1141 /* Deleting can happen also only on the input atoms,
1142 * not necessarily always on the rtp entry.
1144 if (!((*hb)[i].hack[j].oname != nullptr &&
1145 (*hb)[i].hack[j].nname != nullptr) &&
1146 !((*hb)[i].hack[j].oname != nullptr &&
1147 (*hb)[i].hack[j].nname == nullptr))
1150 "atom %s not found in buiding block %d%s "
1151 "while combining tdb and rtp",
1152 (*hb)[i].hack[j].oname != nullptr ?
1153 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].a[0],
1154 i+1, *resinfo[i].rtp);
1159 if ( (*hb)[i].hack[j].oname == nullptr)
1162 add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1167 if ( (*hb)[i].hack[j].nname == nullptr)
1169 /* we're deleting */
1172 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1173 *(*restp)[i].atomname[l],
1174 i+1, (*restp)[i].resname);
1176 /* shift the rest */
1177 (*restp)[i].natom--;
1178 for (k = l; k < (*restp)[i].natom; k++)
1180 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
1181 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1182 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
1184 /* give back space */
1185 srenew((*restp)[i].atom, (*restp)[i].natom);
1186 srenew((*restp)[i].atomname, (*restp)[i].natom);
1187 srenew((*restp)[i].cgnr, (*restp)[i].natom);
1189 else /* nname != NULL */
1190 { /* we're replacing */
1193 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1194 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1195 i+1, (*restp)[i].resname);
1197 snew( (*restp)[i].atomname[l], 1);
1198 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
1199 *(*restp)[i].atomname[l] = gmx_strdup((*hb)[i].hack[j].nname);
1200 if ( (*hb)[i].hack[j].cgnr != NOTSET)
1202 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
1212 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1219 return (gmx_strcasecmp(anm, hack->nname) == 0);
1223 if (isdigit(anm[strlen(anm)-1]))
1225 *nr = anm[strlen(anm)-1] - '0';
1231 if (*nr <= 0 || *nr > hack->nr)
1237 return (strlen(anm) == strlen(hack->nname) + 1 &&
1238 gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1243 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1244 t_restp *rptr, t_hackblock *hbr,
1249 char *oldnm, *newnm;
1251 char *start_at, buf[STRLEN];
1253 gmx_bool bReplaceReplace, bFoundInAdd;
1256 oldnm = *pdba->atomname[atind];
1257 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1260 for (j = 0; j < hbr->nhack; j++)
1262 if (hbr->hack[j].oname != nullptr && hbr->hack[j].nname != nullptr &&
1263 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1265 /* This is a replace entry. */
1266 /* Check if we are not replacing a replaced atom. */
1267 bReplaceReplace = FALSE;
1268 for (k = 0; k < hbr->nhack; k++)
1271 hbr->hack[k].oname != nullptr && hbr->hack[k].nname != nullptr &&
1272 gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1274 /* The replace in hack[j] replaces an atom that
1275 * was already replaced in hack[k], we do not want
1276 * second or higher level replaces at this stage.
1278 bReplaceReplace = TRUE;
1281 if (bReplaceReplace)
1283 /* Skip this replace. */
1287 /* This atom still has the old name, rename it */
1288 newnm = hbr->hack[j].nname;
1289 for (k = 0; k < rptr->natom; k++)
1291 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1296 if (k == rptr->natom)
1298 /* The new name is not present in the rtp.
1299 * We need to apply the replace also to the rtp entry.
1302 /* We need to find the add hack that can add this atom
1303 * to find out after which atom it should be added.
1305 bFoundInAdd = FALSE;
1306 for (k = 0; k < hbr->nhack; k++)
1308 if (hbr->hack[k].oname == nullptr &&
1309 hbr->hack[k].nname != nullptr &&
1310 atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1314 start_at = hbr->hack[k].a[0];
1318 sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1321 for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1323 if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1328 if (start_nr == rptr->natom)
1330 gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1331 start_at, rptr->resname, newnm);
1333 /* We can add the atom after atom start_nr */
1334 add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1342 gmx_fatal(FARGS, "Could not find an 'add' entry for atom named '%s' corresponding to the 'replace' entry from atom name '%s' to '%s' for tdb or hdb database of residue type '%s'",
1344 hbr->hack[j].oname, hbr->hack[j].nname,
1351 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1352 oldnm, rptr->resname, resnr, newnm);
1354 /* Rename the atom in pdba */
1355 snew(pdba->atomname[atind], 1);
1356 *pdba->atomname[atind] = gmx_strdup(newnm);
1358 else if (hbr->hack[j].oname != nullptr && hbr->hack[j].nname == nullptr &&
1359 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1361 /* This is a delete entry, check if this atom is present
1362 * in the rtp entry of this residue.
1364 for (k = 0; k < rptr->natom; k++)
1366 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1371 if (k == rptr->natom)
1373 /* This atom is not present in the rtp entry,
1374 * delete is now from the input pdba.
1378 printf("Deleting atom '%s' in residue '%s' %d\n",
1379 oldnm, rptr->resname, resnr);
1381 /* We should free the atom name,
1382 * but it might be used multiple times in the symtab.
1383 * sfree(pdba->atomname[atind]);
1385 for (k = atind+1; k < pdba->nr; k++)
1387 pdba->atom[k-1] = pdba->atom[k];
1388 pdba->atomname[k-1] = pdba->atomname[k];
1389 copy_rvec(x[k], x[k-1]);
1400 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1401 t_atoms *pdba, rvec *x,
1408 for (i = 0; i < pdba->nr; i++)
1410 oldnm = *pdba->atomname[i];
1411 rptr = &restp[pdba->atom[i].resind];
1412 for (j = 0; (j < rptr->natom); j++)
1414 if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1419 if (j == rptr->natom)
1421 /* Not found yet, check if we have to rename this atom */
1422 if (match_atomnames_with_rtp_atom(pdba, x, i,
1423 rptr, &(hb[pdba->atom[i].resind]),
1426 /* We deleted this atom, decrease the atom counter by 1. */
1433 #define NUM_CMAP_ATOMS 5
1434 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms)
1436 int residx, i, j, k;
1439 t_resinfo *resinfo = atoms->resinfo;
1440 int nres = atoms->nres;
1442 int cmap_atomid[NUM_CMAP_ATOMS];
1443 int cmap_chainnum = -1, this_residue_index;
1454 fprintf(stderr, "Making cmap torsions...\n");
1456 /* Most cmap entries use the N atom from the next residue, so the last
1457 * residue should not have its CMAP entry in that case, but for things like
1458 * dipeptides we sometimes define a complete CMAP entry inside a residue,
1459 * and in this case we need to process everything through the last residue.
1461 for (residx = 0; residx < nres; residx++)
1463 /* Add CMAP terms from the list of CMAP interactions */
1464 for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1467 /* Loop over atoms in a candidate CMAP interaction and
1468 * check that they exist, are from the same chain and are
1469 * from residues labelled as protein. */
1470 for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1472 /* Assign the pointer to the name of the next reference atom.
1473 * This can use -/+ labels to refer to previous/next residue.
1475 pname = restp[residx].rb[ebtsCMAP].b[j].a[k];
1476 /* Skip this CMAP entry if it refers to residues before the
1477 * first or after the last residue.
1479 if (((strchr(pname, '-') != nullptr) && (residx == 0)) ||
1480 ((strchr(pname, '+') != nullptr) && (residx == nres-1)))
1486 cmap_atomid[k] = search_atom(pname,
1487 i, atoms, ptr, TRUE);
1488 bAddCMAP = bAddCMAP && (cmap_atomid[k] != -1);
1491 /* This break is necessary, because cmap_atomid[k]
1492 * == -1 cannot be safely used as an index
1493 * into the atom array. */
1496 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1499 cmap_chainnum = resinfo[this_residue_index].chainnum;
1503 /* Does the residue for this atom have the same
1504 * chain number as the residues for previous
1506 bAddCMAP = bAddCMAP &&
1507 cmap_chainnum == resinfo[this_residue_index].chainnum;
1509 /* Here we used to check that the residuetype was protein and
1510 * disable bAddCMAP if that was not the case. However, some
1511 * special residues (say, alanine dipeptides) might not adhere
1512 * to standard naming, and if we start calling them normal
1513 * protein residues the user will be bugged to select termini.
1515 * Instead, I believe that the right course of action is to
1516 * keep the CMAP interaction if it is present in the RTP file
1517 * and we correctly identified all atoms (which is the case
1524 add_cmap_param(psb, cmap_atomid[0], cmap_atomid[1], cmap_atomid[2], cmap_atomid[3], cmap_atomid[4], restp[residx].rb[ebtsCMAP].b[j].s);
1528 if (residx < nres-1)
1530 while (atoms->atom[i].resind < residx+1)
1536 /* Start the next residue */
1540 scrub_charge_groups(int *cgnr, int natoms)
1544 for (i = 0; i < natoms; i++)
1551 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1552 t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1553 int nrtp, t_restp rtp[],
1554 t_restp *restp, t_hackblock *hb,
1555 gmx_bool bAllowMissing,
1556 gmx_bool bVsites, gmx_bool bVsiteAromatics,
1559 int nssbonds, t_ssbond *ssbonds,
1560 real long_bond_dist, real short_bond_dist,
1561 gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1562 gmx_bool bRenumRes, gmx_bool bRTPresname)
1568 t_params plist[F_NRE];
1575 gmx_residuetype_t*rt;
1578 gmx_residuetype_init(&rt);
1582 print_resall(debug, atoms->nres, restp, atype);
1583 dump_hb(debug, atoms->nres, hb);
1587 at2bonds(&(plist[F_BONDS]), hb,
1589 long_bond_dist, short_bond_dist);
1591 /* specbonds: disulphide bonds & heme-his */
1592 do_ssbonds(&(plist[F_BONDS]),
1593 atoms, nssbonds, ssbonds,
1596 nmissat = name2type(atoms, &cgnr, restp, rt);
1601 fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1606 gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1611 /* Cleanup bonds (sort and rm doubles) */
1612 clean_bonds(&(plist[F_BONDS]));
1614 snew(vsite_type, atoms->nr);
1615 for (i = 0; i < atoms->nr; i++)
1617 vsite_type[i] = NOTSET;
1621 /* determine which atoms will be vsites and add dummy masses
1622 also renumber atom numbers in plist[0..F_NRE]! */
1623 do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1624 &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1627 /* Make Angles and Dihedrals */
1628 fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1629 snew(excls, atoms->nr);
1630 init_nnb(&nnb, atoms->nr, 4);
1631 gen_nnb(&nnb, plist);
1632 print_nnb(&nnb, "NNB");
1633 gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1639 gen_cmap(&(plist[F_CMAP]), restp, atoms);
1640 if (plist[F_CMAP].nr > 0)
1642 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1647 /* set mass of all remaining hydrogen atoms */
1650 do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1654 /* Cleanup bonds (sort and rm doubles) */
1655 /* clean_bonds(&(plist[F_BONDS]));*/
1658 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1659 " %4d pairs, %4d bonds and %4d virtual sites\n",
1660 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1661 plist[F_LJ14].nr, plist[F_BONDS].nr,
1662 plist[F_VSITE2].nr +
1663 plist[F_VSITE3].nr +
1664 plist[F_VSITE3FD].nr +
1665 plist[F_VSITE3FAD].nr +
1666 plist[F_VSITE3OUT].nr +
1667 plist[F_VSITE4FD].nr +
1668 plist[F_VSITE4FDN].nr );
1670 print_sums(atoms, FALSE);
1672 if (FALSE == bChargeGroups)
1674 scrub_charge_groups(cgnr, atoms->nr);
1679 for (i = 0; i < atoms->nres; i++)
1681 atoms->resinfo[i].nr = i + 1;
1682 atoms->resinfo[i].ic = ' ';
1688 fprintf(stderr, "Writing topology\n");
1689 /* We can copy the bonded types from the first restp,
1690 * since the types have to be identical for all residues in one molecule.
1692 for (i = 0; i < ebtsNR; i++)
1694 bts[i] = restp[0].rb[i].type;
1696 write_top(top_file, posre_fn, molname,
1698 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1702 free_t_hackblock(atoms->nres, &hb);
1703 free_t_restp(atoms->nres, &restp);
1704 gmx_residuetype_destroy(rt);
1706 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1708 for (i = 0; i < F_NRE; i++)
1710 sfree(plist[i].param);
1712 for (i = 0; i < atoms->nr; i++)