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, 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.
45 #include "gromacs/math/vec.h"
47 #include "gromacs/utility/futil.h"
49 #include "gpp_nextnb.h"
56 #include "gromacs/utility/cstringutil.h"
57 #include "gromacs/fileio/pdbio.h"
59 #include "gromacs/fileio/filenm.h"
60 #include "gen_vsite.h"
63 #include "fflibutil.h"
66 #include "gromacs/fileio/strdb.h"
67 #include "gromacs/topology/residuetypes.h"
68 #include "gromacs/topology/symtab.h"
69 #include "gromacs/utility/exceptions.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/programcontext.h"
72 #include "gromacs/utility/smalloc.h"
74 /* this must correspond to enum in pdb2top.h */
75 const char *hh[ehisNR] = { "HISD", "HISE", "HISH", "HIS1" };
77 static int missing_atoms(t_restp *rp, int resind, t_atoms *at, int i0, int i)
84 for (j = 0; j < rp->natom; j++)
86 name = *(rp->atomname[j]);
88 for (k = i0; k < i; k++)
90 bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]), name));
95 fprintf(stderr, "\nWARNING: "
96 "atom %s is missing in residue %s %d in the pdb file\n",
97 name, *(at->resinfo[resind].name), at->resinfo[resind].nr);
98 if (name[0] == 'H' || name[0] == 'h')
100 fprintf(stderr, " You might need to add atom %s to the hydrogen database of building block %s\n"
101 " in the file %s.hdb (see the manual)\n",
102 name, *(at->resinfo[resind].rtp), rp->filebase);
104 fprintf(stderr, "\n");
111 gmx_bool is_int(double x)
113 const double tol = 1e-4;
122 return (fabs(x-ix) < tol);
125 static void swap_strings(char **s, int i, int j)
135 choose_ff(const char *ffsel,
136 char *forcefield, int ff_maxlen,
137 char *ffdir, int ffdir_maxlen)
140 char **ffdirs, **ffs, **ffs_dir, *ptr;
141 int i, j, sel, cwdsel, nfound;
142 char buf[STRLEN], **desc;
146 nff = fflib_search_file_in_dirend(fflib_forcefield_itp(),
147 fflib_forcefield_dir_ext(),
152 gmx_fatal(FARGS, "No force fields found (files with name '%s' in subdirectories ending on '%s')",
153 fflib_forcefield_itp(), fflib_forcefield_dir_ext());
156 /* Replace with unix path separators */
157 if (DIR_SEPARATOR != '/')
159 for (i = 0; i < nff; i++)
161 while ( (ptr = strchr(ffdirs[i], DIR_SEPARATOR)) != NULL)
168 /* Store the force field names in ffs */
171 for (i = 0; i < nff; i++)
173 /* Remove the path from the ffdir name - use our unix standard here! */
174 ptr = strrchr(ffdirs[i], '/');
177 ffs[i] = strdup(ffdirs[i]);
178 ffs_dir[i] = low_gmxlibfn(ffdirs[i], FALSE, FALSE);
179 if (ffs_dir[i] == NULL)
181 gmx_fatal(FARGS, "Can no longer find file '%s'", ffdirs[i]);
186 ffs[i] = strdup(ptr+1);
187 ffs_dir[i] = strdup(ffdirs[i]);
189 ffs_dir[i][strlen(ffs_dir[i])-strlen(ffs[i])-1] = '\0';
190 /* Remove the extension from the ffdir name */
191 ffs[i][strlen(ffs[i])-strlen(fflib_forcefield_dir_ext())] = '\0';
199 for (i = 0; i < nff; i++)
201 if (strcmp(ffs[i], ffsel) == 0)
203 /* Matching ff name */
207 if (strncmp(ffs_dir[i], ".", 1) == 0)
224 "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
225 "current directory. Use interactive selection (not the -ff option) if\n"
226 "you would prefer a different one.\n", ffsel, nfound);
231 "Force field '%s' occurs in %d places, but not in the current directory.\n"
232 "Run without the -ff switch and select the force field interactively.", ffsel, nfound);
235 else if (nfound == 0)
237 gmx_fatal(FARGS, "Could not find force field '%s' in current directory, install tree or GMXDATA path.", ffsel);
243 for (i = 0; (i < nff); i++)
245 sprintf(buf, "%s%c%s%s%c%s",
246 ffs_dir[i], DIR_SEPARATOR,
247 ffs[i], fflib_forcefield_dir_ext(), DIR_SEPARATOR,
248 fflib_forcefield_doc());
251 /* We don't use fflib_open, because we don't want printf's */
252 fp = gmx_ffopen(buf, "r");
253 snew(desc[i], STRLEN);
254 get_a_line(fp, desc[i], STRLEN);
259 desc[i] = strdup(ffs[i]);
262 /* Order force fields from the same dir alphabetically
263 * and put deprecated force fields at the end.
265 for (i = 0; (i < nff); i++)
267 for (j = i+1; (j < nff); j++)
269 if (strcmp(ffs_dir[i], ffs_dir[j]) == 0 &&
270 ((desc[i][0] == '[' && desc[j][0] != '[') ||
271 ((desc[i][0] == '[' || desc[j][0] != '[') &&
272 gmx_strcasecmp(desc[i], desc[j]) > 0)))
274 swap_strings(ffdirs, i, j);
275 swap_strings(ffs, i, j);
276 swap_strings(desc, i, j);
281 printf("\nSelect the Force Field:\n");
282 for (i = 0; (i < nff); i++)
284 if (i == 0 || strcmp(ffs_dir[i-1], ffs_dir[i]) != 0)
286 if (strcmp(ffs_dir[i], ".") == 0)
288 printf("From current directory:\n");
292 printf("From '%s':\n", ffs_dir[i]);
295 printf("%2d: %s\n", i+1, desc[i]);
303 pret = fgets(buf, STRLEN, stdin);
307 sel = strtol(buf, NULL, 10);
311 while (pret == NULL || (sel < 0) || (sel >= nff));
313 /* Check for a current limitation of the fflib code.
314 * It will always read from the first ff directory in the list.
315 * This check assumes that the order of ffs matches the order
316 * in which fflib_open searches ff library files.
318 for (i = 0; i < sel; i++)
320 if (strcmp(ffs[i], ffs[sel]) == 0)
322 gmx_fatal(FARGS, "Can only select the first of multiple force field entries with directory name '%s%s' in the list. If you want to use the next entry, run pdb2gmx in a different directory or rename or move the force field directory present in the current working directory.",
323 ffs[sel], fflib_forcefield_dir_ext());
332 if (strlen(ffs[sel]) >= (size_t)ff_maxlen)
334 gmx_fatal(FARGS, "Length of force field name (%d) >= maxlen (%d)",
335 strlen(ffs[sel]), ff_maxlen);
337 strcpy(forcefield, ffs[sel]);
339 if (strlen(ffdirs[sel]) >= (size_t)ffdir_maxlen)
341 gmx_fatal(FARGS, "Length of force field dir (%d) >= maxlen (%d)",
342 strlen(ffdirs[sel]), ffdir_maxlen);
344 strcpy(ffdir, ffdirs[sel]);
346 for (i = 0; (i < nff); i++)
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)
374 else if (strcmp(wmsel, "select") != 0)
376 *watermodel = 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",
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, NULL, 10);
426 while (pret == NULL || sel < 0 || sel > nwm);
434 *watermodel = strdup(model[sel]);
437 for (i = 0; i < nwm; i++)
444 static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
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 = get_atomtype_massA(restp[resind].atom[j].type,
496 cg = restp[resind].cgnr[j];
497 /* A charge group number -1 signals a separate charge group
500 if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) )
509 fprintf(debug, "atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
510 i+1, *(at->atomname[i]), curcg, qt, is_int(qt));
521 at->atom[i].typeB = at->atom[i].type;
522 at->atom[i].qB = at->atom[i].q;
523 at->atom[i].mB = at->atom[i].m;
525 nmissat += missing_atoms(&restp[resind], resind, at, i0, i);
530 static void print_top_heavy_H(FILE *out, real mHmult)
534 fprintf(out, "; Using deuterium instead of hydrogen\n\n");
536 else if (mHmult == 4.0)
538 fprintf(out, "#define HEAVY_H\n\n");
540 else if (mHmult != 1.0)
542 fprintf(stderr, "WARNING: unsupported proton mass multiplier (%g) "
543 "in pdb2top\n", mHmult);
547 void print_top_comment(FILE *out,
548 const char *filename,
552 char ffdir_parent[STRLEN];
555 nice_header(out, filename);
556 fprintf(out, ";\tThis is a %s topology file\n;\n", bITP ? "include" : "standalone");
559 gmx::BinaryInformationSettings settings;
560 settings.generatedByHeader(true);
561 settings.linePrefix(";\t");
562 gmx::printBinaryInformation(out, gmx::getProgramContext(), settings);
564 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
566 if (strchr(ffdir, '/') == NULL)
568 fprintf(out, ";\tForce field was read from the standard Gromacs share directory.\n;\n\n");
570 else if (ffdir[0] == '.')
572 fprintf(out, ";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
576 strncpy(ffdir_parent, ffdir, STRLEN-1);
577 ffdir_parent[STRLEN-1] = '\0'; /*make sure it is 0-terminated even for long string*/
578 p = strrchr(ffdir_parent, '/');
583 ";\tForce field data was read from:\n"
587 ";\tThis might be a non-standard force field location. When you use this topology, the\n"
588 ";\tforce field must either be present in the current directory, or the location\n"
589 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
594 void print_top_header(FILE *out, const char *filename,
595 gmx_bool bITP, const char *ffdir, real mHmult)
599 print_top_comment(out, filename, ffdir, bITP);
601 print_top_heavy_H(out, mHmult);
602 fprintf(out, "; Include forcefield parameters\n");
604 p = strrchr(ffdir, '/');
605 p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
607 fprintf(out, "#include \"%s/%s\"\n\n", p, fflib_forcefield_itp());
610 static void print_top_posre(FILE *out, const char *pr)
612 fprintf(out, "; Include Position restraint file\n");
613 fprintf(out, "#ifdef POSRES\n");
614 fprintf(out, "#include \"%s\"\n", pr);
615 fprintf(out, "#endif\n\n");
618 static void print_top_water(FILE *out, const char *ffdir, const char *water)
623 fprintf(out, "; Include water topology\n");
625 p = strrchr(ffdir, '/');
626 p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
627 fprintf(out, "#include \"%s/%s.itp\"\n", p, water);
630 fprintf(out, "#ifdef POSRES_WATER\n");
631 fprintf(out, "; Position restraint for each water oxygen\n");
632 fprintf(out, "[ position_restraints ]\n");
633 fprintf(out, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
634 fprintf(out, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
635 fprintf(out, "#endif\n");
638 sprintf(buf, "%s/ions.itp", p);
640 if (fflib_fexist(buf))
642 fprintf(out, "; Include topology for ions\n");
643 fprintf(out, "#include \"%s\"\n", buf);
648 static void print_top_system(FILE *out, const char *title)
650 fprintf(out, "[ %s ]\n", dir2str(d_system));
651 fprintf(out, "; Name\n");
652 fprintf(out, "%s\n\n", title[0] ? title : "Protein");
655 void print_top_mols(FILE *out,
656 const char *title, const char *ffdir, const char *water,
657 int nincl, char **incls, int nmol, t_mols *mols)
664 fprintf(out, "; Include chain topologies\n");
665 for (i = 0; (i < nincl); i++)
667 incl = strrchr(incls[i], DIR_SEPARATOR);
674 /* Remove the path from the include name */
677 fprintf(out, "#include \"%s\"\n", incl);
684 print_top_water(out, ffdir, water);
686 print_top_system(out, title);
690 fprintf(out, "[ %s ]\n", dir2str(d_molecules));
691 fprintf(out, "; %-15s %5s\n", "Compound", "#mols");
692 for (i = 0; (i < nmol); i++)
694 fprintf(out, "%-15s %5d\n", mols[i].name, mols[i].nr);
699 void write_top(FILE *out, char *pr, char *molname,
700 t_atoms *at, gmx_bool bRTPresname,
701 int bts[], t_params plist[], t_excls excls[],
702 gpp_atomtype_t atype, int *cgnr, int nrexcl)
703 /* NOTE: nrexcl is not the size of *excl! */
705 if (at && atype && cgnr)
707 fprintf(out, "[ %s ]\n", dir2str(d_moleculetype));
708 fprintf(out, "; %-15s %5s\n", "Name", "nrexcl");
709 fprintf(out, "%-15s %5d\n\n", molname ? molname : "Protein", nrexcl);
711 print_atoms(out, atype, at, cgnr, bRTPresname);
712 print_bondeds(out, at->nr, d_bonds, F_BONDS, bts[ebtsBONDS], plist);
713 print_bondeds(out, at->nr, d_constraints, F_CONSTR, 0, plist);
714 print_bondeds(out, at->nr, d_constraints, F_CONSTRNC, 0, plist);
715 print_bondeds(out, at->nr, d_pairs, F_LJ14, 0, plist);
716 print_excl(out, at->nr, excls);
717 print_bondeds(out, at->nr, d_angles, F_ANGLES, bts[ebtsANGLES], plist);
718 print_bondeds(out, at->nr, d_dihedrals, F_PDIHS, bts[ebtsPDIHS], plist);
719 print_bondeds(out, at->nr, d_dihedrals, F_IDIHS, bts[ebtsIDIHS], plist);
720 print_bondeds(out, at->nr, d_cmap, F_CMAP, bts[ebtsCMAP], plist);
721 print_bondeds(out, at->nr, d_polarization, F_POLARIZATION, 0, plist);
722 print_bondeds(out, at->nr, d_thole_polarization, F_THOLE_POL, 0, plist);
723 print_bondeds(out, at->nr, d_vsites2, F_VSITE2, 0, plist);
724 print_bondeds(out, at->nr, d_vsites3, F_VSITE3, 0, plist);
725 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FD, 0, plist);
726 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FAD, 0, plist);
727 print_bondeds(out, at->nr, d_vsites3, F_VSITE3OUT, 0, plist);
728 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FD, 0, plist);
729 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FDN, 0, plist);
733 print_top_posre(out, pr);
738 static atom_id search_res_atom(const char *type, int resind,
740 const char *bondtype, gmx_bool bAllowMissing)
744 for (i = 0; (i < atoms->nr); i++)
746 if (atoms->atom[i].resind == resind)
748 return search_atom(type, i, atoms, bondtype, bAllowMissing);
755 static void do_ssbonds(t_params *ps, t_atoms *atoms,
756 int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
761 for (i = 0; (i < nssbonds); i++)
763 ri = ssbonds[i].res1;
764 rj = ssbonds[i].res2;
765 ai = search_res_atom(ssbonds[i].a1, ri, atoms,
766 "special bond", bAllowMissing);
767 aj = search_res_atom(ssbonds[i].a2, rj, atoms,
768 "special bond", bAllowMissing);
769 if ((ai == NO_ATID) || (aj == NO_ATID))
771 gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
772 ssbonds[i].a1, ssbonds[i].a2);
774 add_param(ps, ai, aj, NULL, NULL);
778 static void at2bonds(t_params *psb, t_hackblock *hb,
781 real long_bond_dist, real short_bond_dist)
785 real dist2, long_bond_dist2, short_bond_dist2;
788 long_bond_dist2 = sqr(long_bond_dist);
789 short_bond_dist2 = sqr(short_bond_dist);
800 fprintf(stderr, "Making bonds...\n");
802 for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
804 /* add bonds from list of bonded interactions */
805 for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
807 /* Unfortunately we can not issue errors or warnings
808 * for missing atoms in bonds, as the hydrogens and terminal atoms
809 * have not been added yet.
811 ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[0], i, atoms,
813 aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[1], i, atoms,
815 if (ai != NO_ATID && aj != NO_ATID)
817 dist2 = distance2(x[ai], x[aj]);
818 if (dist2 > long_bond_dist2)
820 fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
821 ai+1, aj+1, sqrt(dist2));
823 else if (dist2 < short_bond_dist2)
825 fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
826 ai+1, aj+1, sqrt(dist2));
828 add_param(psb, ai, aj, NULL, hb[resind].rb[ebtsBONDS].b[j].s);
831 /* add bonds from list of hacks (each added atom gets a bond) */
832 while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
834 for (j = 0; j < hb[resind].nhack; j++)
836 if ( ( hb[resind].hack[j].tp > 0 ||
837 hb[resind].hack[j].oname == NULL ) &&
838 strcmp(hb[resind].hack[j].a[0], *(atoms->atomname[i])) == 0)
840 switch (hb[resind].hack[j].tp)
842 case 9: /* COOH terminus */
843 add_param(psb, i, i+1, NULL, NULL); /* C-O */
844 add_param(psb, i, i+2, NULL, NULL); /* C-OA */
845 add_param(psb, i+2, i+3, NULL, NULL); /* OA-H */
848 for (k = 0; (k < hb[resind].hack[j].nr); k++)
850 add_param(psb, i, i+k+1, NULL, NULL);
857 /* we're now at the start of the next residue */
861 static int pcompar(const void *a, const void *b)
868 d = pa->a[0] - pb->a[0];
871 d = pa->a[1] - pb->a[1];
875 return strlen(pb->s) - strlen(pa->s);
883 static void clean_bonds(t_params *ps)
890 /* swap atomnumbers in bond if first larger than second: */
891 for (i = 0; (i < ps->nr); i++)
893 if (ps->param[i].a[1] < ps->param[i].a[0])
895 a = ps->param[i].a[0];
896 ps->param[i].a[0] = ps->param[i].a[1];
897 ps->param[i].a[1] = a;
902 qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
904 /* remove doubles, keep the first one always. */
906 for (i = 1; (i < ps->nr); i++)
908 if ((ps->param[i].a[0] != ps->param[j-1].a[0]) ||
909 (ps->param[i].a[1] != ps->param[j-1].a[1]) )
913 cp_param(&(ps->param[j]), &(ps->param[i]));
918 fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
923 fprintf(stderr, "No bonds\n");
927 void print_sums(t_atoms *atoms, gmx_bool bSystem)
935 where = " in system";
944 for (i = 0; (i < atoms->nr); i++)
946 m += atoms->atom[i].m;
947 qtot += atoms->atom[i].q;
950 fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
951 fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
954 static void check_restp_type(const char *name, int t1, int t2)
958 gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
962 static void check_restp_types(t_restp *r0, t_restp *r1)
966 check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
967 check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
968 check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
969 check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
971 for (i = 0; i < ebtsNR; i++)
973 check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
977 void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
981 const char *Hnum = "123456";
985 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
987 * restp->atomname[at_start], resnr, restp->resname);
989 strcpy(buf, hack->nname);
990 buf[strlen(buf)+1] = '\0';
993 buf[strlen(buf)] = '-';
996 restp->natom += hack->nr;
997 srenew(restp->atom, restp->natom);
998 srenew(restp->atomname, restp->natom);
999 srenew(restp->cgnr, restp->natom);
1000 /* shift the rest */
1001 for (k = restp->natom-1; k > at_start+hack->nr; k--)
1004 restp->atom [k - hack->nr];
1005 restp->atomname[k] =
1006 restp->atomname[k - hack->nr];
1008 restp->cgnr [k - hack->nr];
1011 for (k = 0; k < hack->nr; k++)
1013 /* set counter in atomname */
1016 buf[strlen(buf)-1] = Hnum[k];
1018 snew( restp->atomname[at_start+1+k], 1);
1019 restp->atom [at_start+1+k] = *hack->atom;
1020 *restp->atomname[at_start+1+k] = strdup(buf);
1021 if (hack->cgnr != NOTSET)
1023 restp->cgnr [at_start+1+k] = hack->cgnr;
1027 restp->cgnr [at_start+1+k] = restp->cgnr[at_start];
1032 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1033 int nrtp, t_restp rtp[],
1034 int nres, t_resinfo *resinfo,
1036 t_hackblock **ntdb, t_hackblock **ctdb,
1047 /* first the termini */
1048 for (i = 0; i < nterpairs; i++)
1050 if (rn[i] >= 0 && ntdb[i] != NULL)
1052 copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1054 if (rc[i] >= 0 && ctdb[i] != NULL)
1056 merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1060 /* then the whole rtp */
1061 for (i = 0; i < nres; i++)
1063 /* Here we allow a mismatch of one character when looking for the rtp entry.
1064 * For such a mismatch there should be only one mismatching name.
1065 * This is mainly useful for small molecules such as ions.
1066 * Note that this will usually not work for protein, DNA and RNA,
1067 * since there the residue names should be listed in residuetypes.dat
1068 * and an error will have been generated earlier in the process.
1070 key = *resinfo[i].rtp;
1071 snew(resinfo[i].rtp, 1);
1072 *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1073 res = get_restp(*resinfo[i].rtp, nrtp, rtp);
1074 copy_t_restp(res, &(*restp)[i]);
1076 /* Check that we do not have different bonded types in one molecule */
1077 check_restp_types(&(*restp)[0], &(*restp)[i]);
1080 for (j = 0; j < nterpairs && tern == -1; j++)
1088 for (j = 0; j < nterpairs && terc == -1; j++)
1095 bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1097 if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1098 (terc >= 0 && ctdb[terc] == NULL)))
1100 gmx_fatal(FARGS, "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).");
1102 if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1103 (terc >= 0 && ctdb[terc]->nhack == 0)))
1105 gmx_fatal(FARGS, "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.");
1109 /* now perform t_hack's on t_restp's,
1110 i.e. add's and deletes from termini database will be
1111 added to/removed from residue topology
1112 we'll do this on one big dirty loop, so it won't make easy reading! */
1113 for (i = 0; i < nres; i++)
1115 for (j = 0; j < (*hb)[i].nhack; j++)
1117 if ( (*hb)[i].hack[j].nr)
1119 /* find atom in restp */
1120 for (l = 0; l < (*restp)[i].natom; l++)
1122 if ( ( (*hb)[i].hack[j].oname == NULL &&
1123 strcmp((*hb)[i].hack[j].a[0], *(*restp)[i].atomname[l]) == 0 ) ||
1124 ( (*hb)[i].hack[j].oname != NULL &&
1125 strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1130 if (l == (*restp)[i].natom)
1132 /* If we are doing an atom rename only, we don't need
1133 * to generate a fatal error if the old name is not found
1136 /* Deleting can happen also only on the input atoms,
1137 * not necessarily always on the rtp entry.
1139 if (!((*hb)[i].hack[j].oname != NULL &&
1140 (*hb)[i].hack[j].nname != NULL) &&
1141 !((*hb)[i].hack[j].oname != NULL &&
1142 (*hb)[i].hack[j].nname == NULL))
1145 "atom %s not found in buiding block %d%s "
1146 "while combining tdb and rtp",
1147 (*hb)[i].hack[j].oname != NULL ?
1148 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].a[0],
1149 i+1, *resinfo[i].rtp);
1154 if ( (*hb)[i].hack[j].oname == NULL)
1157 add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1162 if ( (*hb)[i].hack[j].nname == NULL)
1164 /* we're deleting */
1167 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1168 *(*restp)[i].atomname[l],
1169 i+1, (*restp)[i].resname);
1171 /* shift the rest */
1172 (*restp)[i].natom--;
1173 for (k = l; k < (*restp)[i].natom; k++)
1175 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
1176 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1177 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
1179 /* give back space */
1180 srenew((*restp)[i].atom, (*restp)[i].natom);
1181 srenew((*restp)[i].atomname, (*restp)[i].natom);
1182 srenew((*restp)[i].cgnr, (*restp)[i].natom);
1184 else /* nname != NULL */
1185 { /* we're replacing */
1188 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1189 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1190 i+1, (*restp)[i].resname);
1192 snew( (*restp)[i].atomname[l], 1);
1193 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
1194 *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1195 if ( (*hb)[i].hack[j].cgnr != NOTSET)
1197 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
1207 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1214 return (gmx_strcasecmp(anm, hack->nname) == 0);
1218 if (isdigit(anm[strlen(anm)-1]))
1220 *nr = anm[strlen(anm)-1] - '0';
1226 if (*nr <= 0 || *nr > hack->nr)
1232 return (strlen(anm) == strlen(hack->nname) + 1 &&
1233 gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1238 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1239 t_restp *rptr, t_hackblock *hbr,
1244 char *oldnm, *newnm;
1246 char *start_at, buf[STRLEN];
1248 gmx_bool bReplaceReplace, bFoundInAdd;
1251 oldnm = *pdba->atomname[atind];
1252 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1255 for (j = 0; j < hbr->nhack; j++)
1257 if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1258 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1260 /* This is a replace entry. */
1261 /* Check if we are not replacing a replaced atom. */
1262 bReplaceReplace = FALSE;
1263 for (k = 0; k < hbr->nhack; k++)
1266 hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1267 gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1269 /* The replace in hack[j] replaces an atom that
1270 * was already replaced in hack[k], we do not want
1271 * second or higher level replaces at this stage.
1273 bReplaceReplace = TRUE;
1276 if (bReplaceReplace)
1278 /* Skip this replace. */
1282 /* This atom still has the old name, rename it */
1283 newnm = hbr->hack[j].nname;
1284 for (k = 0; k < rptr->natom; k++)
1286 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1291 if (k == rptr->natom)
1293 /* The new name is not present in the rtp.
1294 * We need to apply the replace also to the rtp entry.
1297 /* We need to find the add hack that can add this atom
1298 * to find out after which atom it should be added.
1300 bFoundInAdd = FALSE;
1301 for (k = 0; k < hbr->nhack; k++)
1303 if (hbr->hack[k].oname == NULL &&
1304 hbr->hack[k].nname != NULL &&
1305 atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1309 start_at = hbr->hack[k].a[0];
1313 sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1316 for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1318 if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1323 if (start_nr == rptr->natom)
1325 gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1326 start_at, rptr->resname, newnm);
1328 /* We can add the atom after atom start_nr */
1329 add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1337 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'",
1339 hbr->hack[j].oname, hbr->hack[j].nname,
1346 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1347 oldnm, rptr->resname, resnr, newnm);
1349 /* Rename the atom in pdba */
1350 snew(pdba->atomname[atind], 1);
1351 *pdba->atomname[atind] = strdup(newnm);
1353 else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1354 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1356 /* This is a delete entry, check if this atom is present
1357 * in the rtp entry of this residue.
1359 for (k = 0; k < rptr->natom; k++)
1361 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1366 if (k == rptr->natom)
1368 /* This atom is not present in the rtp entry,
1369 * delete is now from the input pdba.
1373 printf("Deleting atom '%s' in residue '%s' %d\n",
1374 oldnm, rptr->resname, resnr);
1376 /* We should free the atom name,
1377 * but it might be used multiple times in the symtab.
1378 * sfree(pdba->atomname[atind]);
1380 for (k = atind+1; k < pdba->nr; k++)
1382 pdba->atom[k-1] = pdba->atom[k];
1383 pdba->atomname[k-1] = pdba->atomname[k];
1384 copy_rvec(x[k], x[k-1]);
1395 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1396 t_atoms *pdba, rvec *x,
1403 for (i = 0; i < pdba->nr; i++)
1405 oldnm = *pdba->atomname[i];
1406 rptr = &restp[pdba->atom[i].resind];
1407 for (j = 0; (j < rptr->natom); j++)
1409 if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1414 if (j == rptr->natom)
1416 /* Not found yet, check if we have to rename this atom */
1417 if (match_atomnames_with_rtp_atom(pdba, x, i,
1418 rptr, &(hb[pdba->atom[i].resind]),
1421 /* We deleted this atom, decrease the atom counter by 1. */
1428 #define NUM_CMAP_ATOMS 5
1429 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t *rt)
1431 int residx, i, j, k;
1433 t_resinfo *resinfo = atoms->resinfo;
1434 int nres = atoms->nres;
1436 atom_id cmap_atomid[NUM_CMAP_ATOMS];
1437 int cmap_chainnum = -1, this_residue_index;
1448 fprintf(stderr, "Making cmap torsions...");
1450 /* End loop at nres-1, since the very last residue does not have a +N atom, and
1451 * therefore we get a valgrind invalid 4 byte read error with atom am */
1452 for (residx = 0; residx < nres-1; residx++)
1454 /* Add CMAP terms from the list of CMAP interactions */
1455 for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1458 /* Loop over atoms in a candidate CMAP interaction and
1459 * check that they exist, are from the same chain and are
1460 * from residues labelled as protein. */
1461 for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1463 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
1464 i, atoms, ptr, TRUE);
1465 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1468 /* This break is necessary, because cmap_atomid[k]
1469 * == NO_ATID cannot be safely used as an index
1470 * into the atom array. */
1473 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1476 cmap_chainnum = resinfo[this_residue_index].chainnum;
1480 /* Does the residue for this atom have the same
1481 * chain number as the residues for previous
1483 bAddCMAP = bAddCMAP &&
1484 cmap_chainnum == resinfo[this_residue_index].chainnum;
1486 bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt, *(resinfo[this_residue_index].name));
1491 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);
1495 if (residx < nres-1)
1497 while (atoms->atom[i].resind < residx+1)
1504 /* Start the next residue */
1508 scrub_charge_groups(int *cgnr, int natoms)
1512 for (i = 0; i < natoms; i++)
1519 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1520 t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1521 int nrtp, t_restp rtp[],
1522 t_restp *restp, t_hackblock *hb,
1523 gmx_bool bAllowMissing,
1524 gmx_bool bVsites, gmx_bool bVsiteAromatics,
1527 int nssbonds, t_ssbond *ssbonds,
1528 real long_bond_dist, real short_bond_dist,
1529 gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1530 gmx_bool bRenumRes, gmx_bool bRTPresname)
1536 t_params plist[F_NRE];
1543 gmx_residuetype_t*rt;
1546 gmx_residuetype_init(&rt);
1550 print_resall(debug, atoms->nres, restp, atype);
1551 dump_hb(debug, atoms->nres, hb);
1555 at2bonds(&(plist[F_BONDS]), hb,
1557 long_bond_dist, short_bond_dist);
1559 /* specbonds: disulphide bonds & heme-his */
1560 do_ssbonds(&(plist[F_BONDS]),
1561 atoms, nssbonds, ssbonds,
1564 nmissat = name2type(atoms, &cgnr, atype, restp, rt);
1569 fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1574 gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1579 /* Cleanup bonds (sort and rm doubles) */
1580 clean_bonds(&(plist[F_BONDS]));
1582 snew(vsite_type, atoms->nr);
1583 for (i = 0; i < atoms->nr; i++)
1585 vsite_type[i] = NOTSET;
1589 /* determine which atoms will be vsites and add dummy masses
1590 also renumber atom numbers in plist[0..F_NRE]! */
1591 do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1592 &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1595 /* Make Angles and Dihedrals */
1596 fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1597 snew(excls, atoms->nr);
1598 init_nnb(&nnb, atoms->nr, 4);
1599 gen_nnb(&nnb, plist);
1600 print_nnb(&nnb, "NNB");
1601 gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1607 gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1608 if (plist[F_CMAP].nr > 0)
1610 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1615 /* set mass of all remaining hydrogen atoms */
1618 do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1622 /* Cleanup bonds (sort and rm doubles) */
1623 /* clean_bonds(&(plist[F_BONDS]));*/
1626 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1627 " %4d pairs, %4d bonds and %4d virtual sites\n",
1628 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1629 plist[F_LJ14].nr, plist[F_BONDS].nr,
1630 plist[F_VSITE2].nr +
1631 plist[F_VSITE3].nr +
1632 plist[F_VSITE3FD].nr +
1633 plist[F_VSITE3FAD].nr +
1634 plist[F_VSITE3OUT].nr +
1635 plist[F_VSITE4FD].nr +
1636 plist[F_VSITE4FDN].nr );
1638 print_sums(atoms, FALSE);
1640 if (FALSE == bChargeGroups)
1642 scrub_charge_groups(cgnr, atoms->nr);
1647 for (i = 0; i < atoms->nres; i++)
1649 atoms->resinfo[i].nr = i + 1;
1650 atoms->resinfo[i].ic = ' ';
1656 fprintf(stderr, "Writing topology\n");
1657 /* We can copy the bonded types from the first restp,
1658 * since the types have to be identical for all residues in one molecule.
1660 for (i = 0; i < ebtsNR; i++)
1662 bts[i] = restp[0].rb[i].type;
1664 write_top(top_file, posre_fn, molname,
1666 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1670 free_t_hackblock(atoms->nres, &hb);
1671 free_t_restp(atoms->nres, &restp);
1672 gmx_residuetype_destroy(rt);
1674 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1676 for (i = 0; i < F_NRE; i++)
1678 sfree(plist[i].param);
1680 for (i = 0; i < atoms->nr; i++)