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.
46 #include "gromacs/utility/smalloc.h"
49 #include "gromacs/fileio/futil.h"
50 #include "gmx_fatal.h"
52 #include "gpp_nextnb.h"
59 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/fileio/pdbio.h"
63 #include "gromacs/fileio/filenm.h"
65 #include "gen_vsite.h"
68 #include "fflibutil.h"
71 #include "gromacs/fileio/strdb.h"
72 #include "gromacs/utility/exceptions.h"
73 #include "gromacs/utility/programcontext.h"
75 /* this must correspond to enum in pdb2top.h */
76 const char *hh[ehisNR] = { "HISD", "HISE", "HISH", "HIS1" };
78 static int missing_atoms(t_restp *rp, int resind, t_atoms *at, int i0, int i)
85 for (j = 0; j < rp->natom; j++)
87 name = *(rp->atomname[j]);
89 for (k = i0; k < i; k++)
91 bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]), name));
96 fprintf(stderr, "\nWARNING: "
97 "atom %s is missing in residue %s %d in the pdb file\n",
98 name, *(at->resinfo[resind].name), at->resinfo[resind].nr);
99 if (name[0] == 'H' || name[0] == 'h')
101 fprintf(stderr, " You might need to add atom %s to the hydrogen database of building block %s\n"
102 " in the file %s.hdb (see the manual)\n",
103 name, *(at->resinfo[resind].rtp), rp->filebase);
105 fprintf(stderr, "\n");
112 gmx_bool is_int(double x)
114 const double tol = 1e-4;
123 return (fabs(x-ix) < tol);
126 static void swap_strings(char **s, int i, int j)
136 choose_ff(const char *ffsel,
137 char *forcefield, int ff_maxlen,
138 char *ffdir, int ffdir_maxlen)
141 char **ffdirs, **ffs, **ffs_dir, *ptr;
142 int i, j, sel, cwdsel, nfound;
143 char buf[STRLEN], **desc;
147 nff = fflib_search_file_in_dirend(fflib_forcefield_itp(),
148 fflib_forcefield_dir_ext(),
153 gmx_fatal(FARGS, "No force fields found (files with name '%s' in subdirectories ending on '%s')",
154 fflib_forcefield_itp(), fflib_forcefield_dir_ext());
157 /* Replace with unix path separators */
158 if (DIR_SEPARATOR != '/')
160 for (i = 0; i < nff; i++)
162 while ( (ptr = strchr(ffdirs[i], DIR_SEPARATOR)) != NULL)
169 /* Store the force field names in ffs */
172 for (i = 0; i < nff; i++)
174 /* Remove the path from the ffdir name - use our unix standard here! */
175 ptr = strrchr(ffdirs[i], '/');
178 ffs[i] = strdup(ffdirs[i]);
179 ffs_dir[i] = low_gmxlibfn(ffdirs[i], FALSE, FALSE);
180 if (ffs_dir[i] == NULL)
182 gmx_fatal(FARGS, "Can no longer find file '%s'", ffdirs[i]);
187 ffs[i] = strdup(ptr+1);
188 ffs_dir[i] = strdup(ffdirs[i]);
190 ffs_dir[i][strlen(ffs_dir[i])-strlen(ffs[i])-1] = '\0';
191 /* Remove the extension from the ffdir name */
192 ffs[i][strlen(ffs[i])-strlen(fflib_forcefield_dir_ext())] = '\0';
200 for (i = 0; i < nff; i++)
202 if (strcmp(ffs[i], ffsel) == 0)
204 /* Matching ff name */
208 if (strncmp(ffs_dir[i], ".", 1) == 0)
225 "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
226 "current directory. Use interactive selection (not the -ff option) if\n"
227 "you would prefer a different one.\n", ffsel, nfound);
232 "Force field '%s' occurs in %d places, but not in the current directory.\n"
233 "Run without the -ff switch and select the force field interactively.", ffsel, nfound);
236 else if (nfound == 0)
238 gmx_fatal(FARGS, "Could not find force field '%s' in current directory, install tree or GMXDATA path.", ffsel);
244 for (i = 0; (i < nff); i++)
246 sprintf(buf, "%s%c%s%s%c%s",
247 ffs_dir[i], DIR_SEPARATOR,
248 ffs[i], fflib_forcefield_dir_ext(), DIR_SEPARATOR,
249 fflib_forcefield_doc());
252 /* We don't use fflib_open, because we don't want printf's */
253 fp = gmx_ffopen(buf, "r");
254 snew(desc[i], STRLEN);
255 get_a_line(fp, desc[i], STRLEN);
260 desc[i] = strdup(ffs[i]);
263 /* Order force fields from the same dir alphabetically
264 * and put deprecated force fields at the end.
266 for (i = 0; (i < nff); i++)
268 for (j = i+1; (j < nff); j++)
270 if (strcmp(ffs_dir[i], ffs_dir[j]) == 0 &&
271 ((desc[i][0] == '[' && desc[j][0] != '[') ||
272 ((desc[i][0] == '[' || desc[j][0] != '[') &&
273 gmx_strcasecmp(desc[i], desc[j]) > 0)))
275 swap_strings(ffdirs, i, j);
276 swap_strings(ffs, i, j);
277 swap_strings(desc, i, j);
282 printf("\nSelect the Force Field:\n");
283 for (i = 0; (i < nff); i++)
285 if (i == 0 || strcmp(ffs_dir[i-1], ffs_dir[i]) != 0)
287 if (strcmp(ffs_dir[i], ".") == 0)
289 printf("From current directory:\n");
293 printf("From '%s':\n", ffs_dir[i]);
296 printf("%2d: %s\n", i+1, desc[i]);
304 pret = fgets(buf, STRLEN, stdin);
308 sel = strtol(buf, NULL, 10);
312 while (pret == NULL || (sel < 0) || (sel >= nff));
314 /* Check for a current limitation of the fflib code.
315 * It will always read from the first ff directory in the list.
316 * This check assumes that the order of ffs matches the order
317 * in which fflib_open searches ff library files.
319 for (i = 0; i < sel; i++)
321 if (strcmp(ffs[i], ffs[sel]) == 0)
323 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.",
324 ffs[sel], fflib_forcefield_dir_ext());
333 if (strlen(ffs[sel]) >= (size_t)ff_maxlen)
335 gmx_fatal(FARGS, "Length of force field name (%d) >= maxlen (%d)",
336 strlen(ffs[sel]), ff_maxlen);
338 strcpy(forcefield, ffs[sel]);
340 if (strlen(ffdirs[sel]) >= (size_t)ffdir_maxlen)
342 gmx_fatal(FARGS, "Length of force field dir (%d) >= maxlen (%d)",
343 strlen(ffdirs[sel]), ffdir_maxlen);
345 strcpy(ffdir, ffdirs[sel]);
347 for (i = 0; (i < nff); i++)
358 void choose_watermodel(const char *wmsel, const char *ffdir,
361 const char *fn_watermodels = "watermodels.dat";
362 char fn_list[STRLEN];
369 if (strcmp(wmsel, "none") == 0)
375 else if (strcmp(wmsel, "select") != 0)
377 *watermodel = strdup(wmsel);
382 sprintf(fn_list, "%s%c%s", ffdir, DIR_SEPARATOR, fn_watermodels);
384 if (!fflib_fexist(fn_list))
386 fprintf(stderr, "No file '%s' found, will not include a water model\n",
393 fp = fflib_open(fn_list);
394 printf("\nSelect the Water Model:\n");
397 while (get_a_line(fp, buf, STRLEN))
399 srenew(model, nwm+1);
400 snew(model[nwm], STRLEN);
401 sscanf(buf, "%s%n", model[nwm], &i);
405 fprintf(stderr, "%2d: %s\n", nwm+1, buf+i);
414 fprintf(stderr, "%2d: %s\n", nwm+1, "None");
419 pret = fgets(buf, STRLEN, stdin);
423 sel = strtol(buf, NULL, 10);
427 while (pret == NULL || sel < 0 || sel > nwm);
435 *watermodel = strdup(model[sel]);
438 for (i = 0; i < nwm; i++)
445 static int name2type(t_atoms *at, int **cgnr,
446 t_restp restp[], gmx_residuetype_t rt)
448 int i, j, prevresind, resind, i0, prevcg, cg, curcg;
466 for (i = 0; (i < at->nr); i++)
469 if (at->atom[i].resind != resind)
472 resind = at->atom[i].resind;
473 bProt = gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name));
474 bNterm = bProt && (resind == 0);
477 nmissat += missing_atoms(&restp[prevresind], prevresind, at, i0, i);
481 if (at->atom[i].m == 0)
485 fprintf(debug, "atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
486 i+1, *(at->atomname[i]), curcg, prevcg,
487 j == NOTSET ? NOTSET : restp[resind].cgnr[j]);
491 name = *(at->atomname[i]);
492 j = search_jtype(&restp[resind], name, bNterm);
493 at->atom[i].type = restp[resind].atom[j].type;
494 at->atom[i].q = restp[resind].atom[j].q;
495 at->atom[i].m = restp[resind].atom[j].m;
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);
740 static void do_ssbonds(t_params *ps, t_atoms *atoms,
741 int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
746 for (i = 0; (i < nssbonds); i++)
748 ri = ssbonds[i].res1;
749 rj = ssbonds[i].res2;
750 ai = search_res_atom(ssbonds[i].a1, ri, atoms,
751 "special bond", bAllowMissing);
752 aj = search_res_atom(ssbonds[i].a2, rj, atoms,
753 "special bond", bAllowMissing);
754 if ((ai == NO_ATID) || (aj == NO_ATID))
756 gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
757 ssbonds[i].a1, ssbonds[i].a2);
759 add_param(ps, ai, aj, NULL, NULL);
763 static void at2bonds(t_params *psb, t_hackblock *hb,
766 real long_bond_dist, real short_bond_dist)
770 real dist2, long_bond_dist2, short_bond_dist2;
773 long_bond_dist2 = sqr(long_bond_dist);
774 short_bond_dist2 = sqr(short_bond_dist);
785 fprintf(stderr, "Making bonds...\n");
787 for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
789 /* add bonds from list of bonded interactions */
790 for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
792 /* Unfortunately we can not issue errors or warnings
793 * for missing atoms in bonds, as the hydrogens and terminal atoms
794 * have not been added yet.
796 ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[0], i, atoms,
798 aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[1], i, atoms,
800 if (ai != NO_ATID && aj != NO_ATID)
802 dist2 = distance2(x[ai], x[aj]);
803 if (dist2 > long_bond_dist2)
806 fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
807 ai+1, aj+1, sqrt(dist2));
809 else if (dist2 < short_bond_dist2)
811 fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
812 ai+1, aj+1, sqrt(dist2));
814 add_param(psb, ai, aj, NULL, hb[resind].rb[ebtsBONDS].b[j].s);
817 /* add bonds from list of hacks (each added atom gets a bond) */
818 while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
820 for (j = 0; j < hb[resind].nhack; j++)
822 if ( ( hb[resind].hack[j].tp > 0 ||
823 hb[resind].hack[j].oname == NULL ) &&
824 strcmp(hb[resind].hack[j].a[0], *(atoms->atomname[i])) == 0)
826 switch (hb[resind].hack[j].tp)
828 case 9: /* COOH terminus */
829 add_param(psb, i, i+1, NULL, NULL); /* C-O */
830 add_param(psb, i, i+2, NULL, NULL); /* C-OA */
831 add_param(psb, i+2, i+3, NULL, NULL); /* OA-H */
834 for (k = 0; (k < hb[resind].hack[j].nr); k++)
836 add_param(psb, i, i+k+1, NULL, NULL);
843 /* we're now at the start of the next residue */
847 static int pcompar(const void *a, const void *b)
854 d = pa->a[0] - pb->a[0];
857 d = pa->a[1] - pb->a[1];
861 return strlen(pb->s) - strlen(pa->s);
869 static void clean_bonds(t_params *ps)
876 /* swap atomnumbers in bond if first larger than second: */
877 for (i = 0; (i < ps->nr); i++)
879 if (ps->param[i].a[1] < ps->param[i].a[0])
881 a = ps->param[i].a[0];
882 ps->param[i].a[0] = ps->param[i].a[1];
883 ps->param[i].a[1] = a;
888 qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
890 /* remove doubles, keep the first one always. */
892 for (i = 1; (i < ps->nr); i++)
894 if ((ps->param[i].a[0] != ps->param[j-1].a[0]) ||
895 (ps->param[i].a[1] != ps->param[j-1].a[1]) )
899 cp_param(&(ps->param[j]), &(ps->param[i]));
904 fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
909 fprintf(stderr, "No bonds\n");
913 void print_sums(t_atoms *atoms, gmx_bool bSystem)
921 where = " in system";
930 for (i = 0; (i < atoms->nr); i++)
932 m += atoms->atom[i].m;
933 qtot += atoms->atom[i].q;
936 fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
937 fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
940 static void check_restp_type(const char *name, int t1, int t2)
944 gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
948 static void check_restp_types(t_restp *r0, t_restp *r1)
952 check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
953 check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
954 check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
955 check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
957 for (i = 0; i < ebtsNR; i++)
959 check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
963 void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
967 const char *Hnum = "123456";
971 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
973 * restp->atomname[at_start], resnr, restp->resname);
975 strcpy(buf, hack->nname);
976 buf[strlen(buf)+1] = '\0';
979 buf[strlen(buf)] = '-';
982 restp->natom += hack->nr;
983 srenew(restp->atom, restp->natom);
984 srenew(restp->atomname, restp->natom);
985 srenew(restp->cgnr, restp->natom);
987 for (k = restp->natom-1; k > at_start+hack->nr; k--)
990 restp->atom [k - hack->nr];
992 restp->atomname[k - hack->nr];
994 restp->cgnr [k - hack->nr];
997 for (k = 0; k < hack->nr; k++)
999 /* set counter in atomname */
1002 buf[strlen(buf)-1] = Hnum[k];
1004 snew( restp->atomname[at_start+1+k], 1);
1005 restp->atom [at_start+1+k] = *hack->atom;
1006 *restp->atomname[at_start+1+k] = strdup(buf);
1007 if (hack->cgnr != NOTSET)
1009 restp->cgnr [at_start+1+k] = hack->cgnr;
1013 restp->cgnr [at_start+1+k] = restp->cgnr[at_start];
1018 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1019 int nrtp, t_restp rtp[],
1020 int nres, t_resinfo *resinfo,
1022 t_hackblock **ntdb, t_hackblock **ctdb,
1033 /* first the termini */
1034 for (i = 0; i < nterpairs; i++)
1036 if (rn[i] >= 0 && ntdb[i] != NULL)
1038 copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1040 if (rc[i] >= 0 && ctdb[i] != NULL)
1042 merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1046 /* then the whole rtp */
1047 for (i = 0; i < nres; i++)
1049 /* Here we allow a mismatch of one character when looking for the rtp entry.
1050 * For such a mismatch there should be only one mismatching name.
1051 * This is mainly useful for small molecules such as ions.
1052 * Note that this will usually not work for protein, DNA and RNA,
1053 * since there the residue names should be listed in residuetypes.dat
1054 * and an error will have been generated earlier in the process.
1056 key = *resinfo[i].rtp;
1057 snew(resinfo[i].rtp, 1);
1058 *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1059 res = get_restp(*resinfo[i].rtp, nrtp, rtp);
1060 copy_t_restp(res, &(*restp)[i]);
1062 /* Check that we do not have different bonded types in one molecule */
1063 check_restp_types(&(*restp)[0], &(*restp)[i]);
1066 for (j = 0; j < nterpairs && tern == -1; j++)
1074 for (j = 0; j < nterpairs && terc == -1; j++)
1081 bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1083 if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1084 (terc >= 0 && ctdb[terc] == NULL)))
1086 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).");
1088 if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1089 (terc >= 0 && ctdb[terc]->nhack == 0)))
1091 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.");
1095 /* now perform t_hack's on t_restp's,
1096 i.e. add's and deletes from termini database will be
1097 added to/removed from residue topology
1098 we'll do this on one big dirty loop, so it won't make easy reading! */
1099 for (i = 0; i < nres; i++)
1101 for (j = 0; j < (*hb)[i].nhack; j++)
1103 if ( (*hb)[i].hack[j].nr)
1105 /* find atom in restp */
1106 for (l = 0; l < (*restp)[i].natom; l++)
1108 if ( ( (*hb)[i].hack[j].oname == NULL &&
1109 strcmp((*hb)[i].hack[j].a[0], *(*restp)[i].atomname[l]) == 0 ) ||
1110 ( (*hb)[i].hack[j].oname != NULL &&
1111 strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1116 if (l == (*restp)[i].natom)
1118 /* If we are doing an atom rename only, we don't need
1119 * to generate a fatal error if the old name is not found
1122 /* Deleting can happen also only on the input atoms,
1123 * not necessarily always on the rtp entry.
1125 if (!((*hb)[i].hack[j].oname != NULL &&
1126 (*hb)[i].hack[j].nname != NULL) &&
1127 !((*hb)[i].hack[j].oname != NULL &&
1128 (*hb)[i].hack[j].nname == NULL))
1131 "atom %s not found in buiding block %d%s "
1132 "while combining tdb and rtp",
1133 (*hb)[i].hack[j].oname != NULL ?
1134 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].a[0],
1135 i+1, *resinfo[i].rtp);
1140 if ( (*hb)[i].hack[j].oname == NULL)
1143 add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1148 if ( (*hb)[i].hack[j].nname == NULL)
1150 /* we're deleting */
1153 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1154 *(*restp)[i].atomname[l],
1155 i+1, (*restp)[i].resname);
1157 /* shift the rest */
1158 (*restp)[i].natom--;
1159 for (k = l; k < (*restp)[i].natom; k++)
1161 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
1162 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1163 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
1165 /* give back space */
1166 srenew((*restp)[i].atom, (*restp)[i].natom);
1167 srenew((*restp)[i].atomname, (*restp)[i].natom);
1168 srenew((*restp)[i].cgnr, (*restp)[i].natom);
1170 else /* nname != NULL */
1171 { /* we're replacing */
1174 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1175 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1176 i+1, (*restp)[i].resname);
1178 snew( (*restp)[i].atomname[l], 1);
1179 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
1180 *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1181 if ( (*hb)[i].hack[j].cgnr != NOTSET)
1183 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
1193 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1200 return (gmx_strcasecmp(anm, hack->nname) == 0);
1204 if (isdigit(anm[strlen(anm)-1]))
1206 *nr = anm[strlen(anm)-1] - '0';
1212 if (*nr <= 0 || *nr > hack->nr)
1218 return (strlen(anm) == strlen(hack->nname) + 1 &&
1219 gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1224 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1225 t_restp *rptr, t_hackblock *hbr,
1230 char *oldnm, *newnm;
1232 char *start_at, buf[STRLEN];
1234 gmx_bool bReplaceReplace, bFoundInAdd;
1237 oldnm = *pdba->atomname[atind];
1238 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1241 for (j = 0; j < hbr->nhack; j++)
1243 if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1244 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1246 /* This is a replace entry. */
1247 /* Check if we are not replacing a replaced atom. */
1248 bReplaceReplace = FALSE;
1249 for (k = 0; k < hbr->nhack; k++)
1252 hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1253 gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1255 /* The replace in hack[j] replaces an atom that
1256 * was already replaced in hack[k], we do not want
1257 * second or higher level replaces at this stage.
1259 bReplaceReplace = TRUE;
1262 if (bReplaceReplace)
1264 /* Skip this replace. */
1268 /* This atom still has the old name, rename it */
1269 newnm = hbr->hack[j].nname;
1270 for (k = 0; k < rptr->natom; k++)
1272 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1277 if (k == rptr->natom)
1279 /* The new name is not present in the rtp.
1280 * We need to apply the replace also to the rtp entry.
1283 /* We need to find the add hack that can add this atom
1284 * to find out after which atom it should be added.
1286 bFoundInAdd = FALSE;
1287 for (k = 0; k < hbr->nhack; k++)
1289 if (hbr->hack[k].oname == NULL &&
1290 hbr->hack[k].nname != NULL &&
1291 atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1295 start_at = hbr->hack[k].a[0];
1299 sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1302 for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1304 if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1309 if (start_nr == rptr->natom)
1311 gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1312 start_at, rptr->resname, newnm);
1314 /* We can add the atom after atom start_nr */
1315 add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1323 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'",
1325 hbr->hack[j].oname, hbr->hack[j].nname,
1332 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1333 oldnm, rptr->resname, resnr, newnm);
1335 /* Rename the atom in pdba */
1336 snew(pdba->atomname[atind], 1);
1337 *pdba->atomname[atind] = strdup(newnm);
1339 else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1340 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1342 /* This is a delete entry, check if this atom is present
1343 * in the rtp entry of this residue.
1345 for (k = 0; k < rptr->natom; k++)
1347 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1352 if (k == rptr->natom)
1354 /* This atom is not present in the rtp entry,
1355 * delete is now from the input pdba.
1359 printf("Deleting atom '%s' in residue '%s' %d\n",
1360 oldnm, rptr->resname, resnr);
1362 /* We should free the atom name,
1363 * but it might be used multiple times in the symtab.
1364 * sfree(pdba->atomname[atind]);
1366 for (k = atind+1; k < pdba->nr; k++)
1368 pdba->atom[k-1] = pdba->atom[k];
1369 pdba->atomname[k-1] = pdba->atomname[k];
1370 copy_rvec(x[k], x[k-1]);
1381 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1382 t_atoms *pdba, rvec *x,
1389 for (i = 0; i < pdba->nr; i++)
1391 oldnm = *pdba->atomname[i];
1392 rptr = &restp[pdba->atom[i].resind];
1393 for (j = 0; (j < rptr->natom); j++)
1395 if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1400 if (j == rptr->natom)
1402 /* Not found yet, check if we have to rename this atom */
1403 if (match_atomnames_with_rtp_atom(pdba, x, i,
1404 rptr, &(hb[pdba->atom[i].resind]),
1407 /* We deleted this atom, decrease the atom counter by 1. */
1414 #define NUM_CMAP_ATOMS 5
1415 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms)
1417 int residx, i, j, k;
1420 t_resinfo *resinfo = atoms->resinfo;
1421 int nres = atoms->nres;
1423 atom_id cmap_atomid[NUM_CMAP_ATOMS];
1424 int cmap_chainnum = -1, this_residue_index;
1435 fprintf(stderr, "Making cmap torsions...\n");
1437 /* Most cmap entries use the N atom from the next residue, so the last
1438 * residue should not have its CMAP entry in that case, but for things like
1439 * dipeptides we sometimes define a complete CMAP entry inside a residue,
1440 * and in this case we need to process everything through the last residue.
1442 for (residx = 0; residx < nres; residx++)
1444 /* Add CMAP terms from the list of CMAP interactions */
1445 for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1448 /* Loop over atoms in a candidate CMAP interaction and
1449 * check that they exist, are from the same chain and are
1450 * from residues labelled as protein. */
1451 for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1453 /* Assign the pointer to the name of the next reference atom.
1454 * This can use -/+ labels to refer to previous/next residue.
1456 pname = restp[residx].rb[ebtsCMAP].b[j].a[k];
1457 /* Skip this CMAP entry if it refers to residues before the
1458 * first or after the last residue.
1460 if (((strchr(pname, '-') != NULL) && (residx == 0)) ||
1461 ((strchr(pname, '+') != NULL) && (residx == nres-1)))
1467 cmap_atomid[k] = search_atom(pname,
1468 i, atoms, ptr, TRUE);
1469 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1472 /* This break is necessary, because cmap_atomid[k]
1473 * == NO_ATID cannot be safely used as an index
1474 * into the atom array. */
1477 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1480 cmap_chainnum = resinfo[this_residue_index].chainnum;
1484 /* Does the residue for this atom have the same
1485 * chain number as the residues for previous
1487 bAddCMAP = bAddCMAP &&
1488 cmap_chainnum == resinfo[this_residue_index].chainnum;
1490 /* Here we used to check that the residuetype was protein and
1491 * disable bAddCMAP if that was not the case. However, some
1492 * special residues (say, alanine dipeptides) might not adhere
1493 * to standard naming, and if we start calling them normal
1494 * protein residues the user will be bugged to select termini.
1496 * Instead, I believe that the right course of action is to
1497 * keep the CMAP interaction if it is present in the RTP file
1498 * and we correctly identified all atoms (which is the case
1505 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);
1509 if (residx < nres-1)
1511 while (atoms->atom[i].resind < residx+1)
1517 /* Start the next residue */
1521 scrub_charge_groups(int *cgnr, int natoms)
1525 for (i = 0; i < natoms; i++)
1532 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1533 t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1534 int nrtp, t_restp rtp[],
1535 t_restp *restp, t_hackblock *hb,
1536 gmx_bool bAllowMissing,
1537 gmx_bool bVsites, gmx_bool bVsiteAromatics,
1540 int nssbonds, t_ssbond *ssbonds,
1541 real long_bond_dist, real short_bond_dist,
1542 gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1543 gmx_bool bRenumRes, gmx_bool bRTPresname)
1549 t_params plist[F_NRE];
1556 gmx_residuetype_t rt;
1559 gmx_residuetype_init(&rt);
1563 print_resall(debug, atoms->nres, restp, atype);
1564 dump_hb(debug, atoms->nres, hb);
1568 at2bonds(&(plist[F_BONDS]), hb,
1570 long_bond_dist, short_bond_dist);
1572 /* specbonds: disulphide bonds & heme-his */
1573 do_ssbonds(&(plist[F_BONDS]),
1574 atoms, nssbonds, ssbonds,
1577 nmissat = name2type(atoms, &cgnr, restp, rt);
1582 fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1587 gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1592 /* Cleanup bonds (sort and rm doubles) */
1593 clean_bonds(&(plist[F_BONDS]));
1595 snew(vsite_type, atoms->nr);
1596 for (i = 0; i < atoms->nr; i++)
1598 vsite_type[i] = NOTSET;
1602 /* determine which atoms will be vsites and add dummy masses
1603 also renumber atom numbers in plist[0..F_NRE]! */
1604 do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1605 &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1608 /* Make Angles and Dihedrals */
1609 fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1610 snew(excls, atoms->nr);
1611 init_nnb(&nnb, atoms->nr, 4);
1612 gen_nnb(&nnb, plist);
1613 print_nnb(&nnb, "NNB");
1614 gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1620 gen_cmap(&(plist[F_CMAP]), restp, atoms);
1621 if (plist[F_CMAP].nr > 0)
1623 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1628 /* set mass of all remaining hydrogen atoms */
1631 do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1635 /* Cleanup bonds (sort and rm doubles) */
1636 /* clean_bonds(&(plist[F_BONDS]));*/
1639 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1640 " %4d pairs, %4d bonds and %4d virtual sites\n",
1641 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1642 plist[F_LJ14].nr, plist[F_BONDS].nr,
1643 plist[F_VSITE2].nr +
1644 plist[F_VSITE3].nr +
1645 plist[F_VSITE3FD].nr +
1646 plist[F_VSITE3FAD].nr +
1647 plist[F_VSITE3OUT].nr +
1648 plist[F_VSITE4FD].nr +
1649 plist[F_VSITE4FDN].nr );
1651 print_sums(atoms, FALSE);
1653 if (FALSE == bChargeGroups)
1655 scrub_charge_groups(cgnr, atoms->nr);
1660 for (i = 0; i < atoms->nres; i++)
1662 atoms->resinfo[i].nr = i + 1;
1663 atoms->resinfo[i].ic = ' ';
1669 fprintf(stderr, "Writing topology\n");
1670 /* We can copy the bonded types from the first restp,
1671 * since the types have to be identical for all residues in one molecule.
1673 for (i = 0; i < ebtsNR; i++)
1675 bts[i] = restp[0].rb[i].type;
1677 write_top(top_file, posre_fn, molname,
1679 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1683 free_t_hackblock(atoms->nres, &hb);
1684 free_t_restp(atoms->nres, &restp);
1685 gmx_residuetype_destroy(rt);
1687 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1689 for (i = 0; i < F_NRE; i++)
1691 sfree(plist[i].param);
1693 for (i = 0; i < atoms->nr; i++)