1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
50 #include "gmx_fatal.h"
52 #include "gpp_nextnb.h"
65 #include "gen_vsite.h"
68 #include "fflibutil.h"
72 #include "gromacs/utility/exceptions.h"
73 #include "gromacs/utility/programinfo.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 = 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, gpp_atomtype_t atype,
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 = get_atomtype_massA(restp[resind].atom[j].type,
497 cg = restp[resind].cgnr[j];
498 /* A charge group number -1 signals a separate charge group
501 if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) )
510 fprintf(debug, "atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
511 i+1, *(at->atomname[i]), curcg, qt, is_int(qt));
522 at->atom[i].typeB = at->atom[i].type;
523 at->atom[i].qB = at->atom[i].q;
524 at->atom[i].mB = at->atom[i].m;
526 nmissat += missing_atoms(&restp[resind], resind, at, i0, i);
531 static void print_top_heavy_H(FILE *out, real mHmult)
535 fprintf(out, "; Using deuterium instead of hydrogen\n\n");
537 else if (mHmult == 4.0)
539 fprintf(out, "#define HEAVY_H\n\n");
541 else if (mHmult != 1.0)
543 fprintf(stderr, "WARNING: unsupported proton mass multiplier (%g) "
544 "in pdb2top\n", mHmult);
548 void print_top_comment(FILE *out,
549 const char *filename,
553 char ffdir_parent[STRLEN];
556 nice_header(out, filename);
557 fprintf(out, ";\tThis is a %s topology file\n;\n", bITP ? "include" : "standalone");
560 gmx::BinaryInformationSettings settings;
561 settings.generatedByHeader(true);
562 settings.linePrefix(";\t");
563 gmx::printBinaryInformation(out, gmx::ProgramInfo::getInstance(),
566 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
568 if (strchr(ffdir, '/') == NULL)
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 == NULL) ? 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 == NULL) ? 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, 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);
740 static atom_id search_res_atom(const char *type, int resind,
742 const char *bondtype, gmx_bool bAllowMissing)
746 for (i = 0; (i < atoms->nr); i++)
748 if (atoms->atom[i].resind == resind)
750 return search_atom(type, i, atoms, bondtype, bAllowMissing);
757 static void do_ssbonds(t_params *ps, t_atoms *atoms,
758 int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
763 for (i = 0; (i < nssbonds); i++)
765 ri = ssbonds[i].res1;
766 rj = ssbonds[i].res2;
767 ai = search_res_atom(ssbonds[i].a1, ri, atoms,
768 "special bond", bAllowMissing);
769 aj = search_res_atom(ssbonds[i].a2, rj, atoms,
770 "special bond", bAllowMissing);
771 if ((ai == NO_ATID) || (aj == NO_ATID))
773 gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
774 ssbonds[i].a1, ssbonds[i].a2);
776 add_param(ps, ai, aj, NULL, NULL);
780 static void at2bonds(t_params *psb, t_hackblock *hb,
783 real long_bond_dist, real short_bond_dist)
787 real dist2, long_bond_dist2, short_bond_dist2;
790 long_bond_dist2 = sqr(long_bond_dist);
791 short_bond_dist2 = sqr(short_bond_dist);
802 fprintf(stderr, "Making bonds...\n");
804 for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
806 /* add bonds from list of bonded interactions */
807 for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
809 /* Unfortunately we can not issue errors or warnings
810 * for missing atoms in bonds, as the hydrogens and terminal atoms
811 * have not been added yet.
813 ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[0], i, atoms,
815 aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[1], i, atoms,
817 if (ai != NO_ATID && aj != NO_ATID)
819 dist2 = distance2(x[ai], x[aj]);
820 if (dist2 > long_bond_dist2)
822 fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
823 ai+1, aj+1, sqrt(dist2));
825 else if (dist2 < short_bond_dist2)
827 fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
828 ai+1, aj+1, sqrt(dist2));
830 add_param(psb, ai, aj, NULL, hb[resind].rb[ebtsBONDS].b[j].s);
833 /* add bonds from list of hacks (each added atom gets a bond) */
834 while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
836 for (j = 0; j < hb[resind].nhack; j++)
838 if ( ( hb[resind].hack[j].tp > 0 ||
839 hb[resind].hack[j].oname == NULL ) &&
840 strcmp(hb[resind].hack[j].a[0], *(atoms->atomname[i])) == 0)
842 switch (hb[resind].hack[j].tp)
844 case 9: /* COOH terminus */
845 add_param(psb, i, i+1, NULL, NULL); /* C-O */
846 add_param(psb, i, i+2, NULL, NULL); /* C-OA */
847 add_param(psb, i+2, i+3, NULL, NULL); /* OA-H */
850 for (k = 0; (k < hb[resind].hack[j].nr); k++)
852 add_param(psb, i, i+k+1, NULL, NULL);
859 /* we're now at the start of the next residue */
863 static int pcompar(const void *a, const void *b)
870 d = pa->a[0] - pb->a[0];
873 d = pa->a[1] - pb->a[1];
877 return strlen(pb->s) - strlen(pa->s);
885 static void clean_bonds(t_params *ps)
892 /* swap atomnumbers in bond if first larger than second: */
893 for (i = 0; (i < ps->nr); i++)
895 if (ps->param[i].a[1] < ps->param[i].a[0])
897 a = ps->param[i].a[0];
898 ps->param[i].a[0] = ps->param[i].a[1];
899 ps->param[i].a[1] = a;
904 qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
906 /* remove doubles, keep the first one always. */
908 for (i = 1; (i < ps->nr); i++)
910 if ((ps->param[i].a[0] != ps->param[j-1].a[0]) ||
911 (ps->param[i].a[1] != ps->param[j-1].a[1]) )
915 cp_param(&(ps->param[j]), &(ps->param[i]));
920 fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
925 fprintf(stderr, "No bonds\n");
929 void print_sums(t_atoms *atoms, gmx_bool bSystem)
937 where = " in system";
946 for (i = 0; (i < atoms->nr); i++)
948 m += atoms->atom[i].m;
949 qtot += atoms->atom[i].q;
952 fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
953 fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
956 static void check_restp_type(const char *name, int t1, int t2)
960 gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
964 static void check_restp_types(t_restp *r0, t_restp *r1)
968 check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
969 check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
970 check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
971 check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
973 for (i = 0; i < ebtsNR; i++)
975 check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
979 void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
983 const char *Hnum = "123456";
987 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
989 * restp->atomname[at_start], resnr, restp->resname);
991 strcpy(buf, hack->nname);
992 buf[strlen(buf)+1] = '\0';
995 buf[strlen(buf)] = '-';
998 restp->natom += hack->nr;
999 srenew(restp->atom, restp->natom);
1000 srenew(restp->atomname, restp->natom);
1001 srenew(restp->cgnr, restp->natom);
1002 /* shift the rest */
1003 for (k = restp->natom-1; k > at_start+hack->nr; k--)
1006 restp->atom [k - hack->nr];
1007 restp->atomname[k] =
1008 restp->atomname[k - hack->nr];
1010 restp->cgnr [k - hack->nr];
1013 for (k = 0; k < hack->nr; k++)
1015 /* set counter in atomname */
1018 buf[strlen(buf)-1] = Hnum[k];
1020 snew( restp->atomname[at_start+1+k], 1);
1021 restp->atom [at_start+1+k] = *hack->atom;
1022 *restp->atomname[at_start+1+k] = strdup(buf);
1023 if (hack->cgnr != NOTSET)
1025 restp->cgnr [at_start+1+k] = hack->cgnr;
1029 restp->cgnr [at_start+1+k] = restp->cgnr[at_start];
1034 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1035 int nrtp, t_restp rtp[],
1036 int nres, t_resinfo *resinfo,
1038 t_hackblock **ntdb, t_hackblock **ctdb,
1049 /* first the termini */
1050 for (i = 0; i < nterpairs; i++)
1052 if (rn[i] >= 0 && ntdb[i] != NULL)
1054 copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1056 if (rc[i] >= 0 && ctdb[i] != NULL)
1058 merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1062 /* then the whole rtp */
1063 for (i = 0; i < nres; i++)
1065 /* Here we allow a mismatch of one character when looking for the rtp entry.
1066 * For such a mismatch there should be only one mismatching name.
1067 * This is mainly useful for small molecules such as ions.
1068 * Note that this will usually not work for protein, DNA and RNA,
1069 * since there the residue names should be listed in residuetypes.dat
1070 * and an error will have been generated earlier in the process.
1072 key = *resinfo[i].rtp;
1073 snew(resinfo[i].rtp, 1);
1074 *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1075 res = get_restp(*resinfo[i].rtp, nrtp, rtp);
1076 copy_t_restp(res, &(*restp)[i]);
1078 /* Check that we do not have different bonded types in one molecule */
1079 check_restp_types(&(*restp)[0], &(*restp)[i]);
1082 for (j = 0; j < nterpairs && tern == -1; j++)
1090 for (j = 0; j < nterpairs && terc == -1; j++)
1097 bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1099 if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1100 (terc >= 0 && ctdb[terc] == NULL)))
1102 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).");
1104 if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1105 (terc >= 0 && ctdb[terc]->nhack == 0)))
1107 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.");
1111 /* now perform t_hack's on t_restp's,
1112 i.e. add's and deletes from termini database will be
1113 added to/removed from residue topology
1114 we'll do this on one big dirty loop, so it won't make easy reading! */
1115 for (i = 0; i < nres; i++)
1117 for (j = 0; j < (*hb)[i].nhack; j++)
1119 if ( (*hb)[i].hack[j].nr)
1121 /* find atom in restp */
1122 for (l = 0; l < (*restp)[i].natom; l++)
1124 if ( ( (*hb)[i].hack[j].oname == NULL &&
1125 strcmp((*hb)[i].hack[j].a[0], *(*restp)[i].atomname[l]) == 0 ) ||
1126 ( (*hb)[i].hack[j].oname != NULL &&
1127 strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1132 if (l == (*restp)[i].natom)
1134 /* If we are doing an atom rename only, we don't need
1135 * to generate a fatal error if the old name is not found
1138 /* Deleting can happen also only on the input atoms,
1139 * not necessarily always on the rtp entry.
1141 if (!((*hb)[i].hack[j].oname != NULL &&
1142 (*hb)[i].hack[j].nname != NULL) &&
1143 !((*hb)[i].hack[j].oname != NULL &&
1144 (*hb)[i].hack[j].nname == NULL))
1147 "atom %s not found in buiding block %d%s "
1148 "while combining tdb and rtp",
1149 (*hb)[i].hack[j].oname != NULL ?
1150 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].a[0],
1151 i+1, *resinfo[i].rtp);
1156 if ( (*hb)[i].hack[j].oname == NULL)
1159 add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1164 if ( (*hb)[i].hack[j].nname == NULL)
1166 /* we're deleting */
1169 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1170 *(*restp)[i].atomname[l],
1171 i+1, (*restp)[i].resname);
1173 /* shift the rest */
1174 (*restp)[i].natom--;
1175 for (k = l; k < (*restp)[i].natom; k++)
1177 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
1178 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1179 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
1181 /* give back space */
1182 srenew((*restp)[i].atom, (*restp)[i].natom);
1183 srenew((*restp)[i].atomname, (*restp)[i].natom);
1184 srenew((*restp)[i].cgnr, (*restp)[i].natom);
1186 else /* nname != NULL */
1187 { /* we're replacing */
1190 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1191 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1192 i+1, (*restp)[i].resname);
1194 snew( (*restp)[i].atomname[l], 1);
1195 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
1196 *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1197 if ( (*hb)[i].hack[j].cgnr != NOTSET)
1199 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
1209 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1216 return (gmx_strcasecmp(anm, hack->nname) == 0);
1220 if (isdigit(anm[strlen(anm)-1]))
1222 *nr = anm[strlen(anm)-1] - '0';
1228 if (*nr <= 0 || *nr > hack->nr)
1234 return (strlen(anm) == strlen(hack->nname) + 1 &&
1235 gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1240 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1241 t_restp *rptr, t_hackblock *hbr,
1246 char *oldnm, *newnm;
1248 char *start_at, buf[STRLEN];
1250 gmx_bool bReplaceReplace, bFoundInAdd;
1253 oldnm = *pdba->atomname[atind];
1254 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1257 for (j = 0; j < hbr->nhack; j++)
1259 if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1260 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1262 /* This is a replace entry. */
1263 /* Check if we are not replacing a replaced atom. */
1264 bReplaceReplace = FALSE;
1265 for (k = 0; k < hbr->nhack; k++)
1268 hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1269 gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1271 /* The replace in hack[j] replaces an atom that
1272 * was already replaced in hack[k], we do not want
1273 * second or higher level replaces at this stage.
1275 bReplaceReplace = TRUE;
1278 if (bReplaceReplace)
1280 /* Skip this replace. */
1284 /* This atom still has the old name, rename it */
1285 newnm = hbr->hack[j].nname;
1286 for (k = 0; k < rptr->natom; k++)
1288 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1293 if (k == rptr->natom)
1295 /* The new name is not present in the rtp.
1296 * We need to apply the replace also to the rtp entry.
1299 /* We need to find the add hack that can add this atom
1300 * to find out after which atom it should be added.
1302 bFoundInAdd = FALSE;
1303 for (k = 0; k < hbr->nhack; k++)
1305 if (hbr->hack[k].oname == NULL &&
1306 hbr->hack[k].nname != NULL &&
1307 atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1311 start_at = hbr->hack[k].a[0];
1315 sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1318 for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1320 if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1325 if (start_nr == rptr->natom)
1327 gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1328 start_at, rptr->resname, newnm);
1330 /* We can add the atom after atom start_nr */
1331 add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1339 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'",
1341 hbr->hack[j].oname, hbr->hack[j].nname,
1348 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1349 oldnm, rptr->resname, resnr, newnm);
1351 /* Rename the atom in pdba */
1352 snew(pdba->atomname[atind], 1);
1353 *pdba->atomname[atind] = strdup(newnm);
1355 else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1356 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1358 /* This is a delete entry, check if this atom is present
1359 * in the rtp entry of this residue.
1361 for (k = 0; k < rptr->natom; k++)
1363 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1368 if (k == rptr->natom)
1370 /* This atom is not present in the rtp entry,
1371 * delete is now from the input pdba.
1375 printf("Deleting atom '%s' in residue '%s' %d\n",
1376 oldnm, rptr->resname, resnr);
1378 /* We should free the atom name,
1379 * but it might be used multiple times in the symtab.
1380 * sfree(pdba->atomname[atind]);
1382 for (k = atind+1; k < pdba->nr; k++)
1384 pdba->atom[k-1] = pdba->atom[k];
1385 pdba->atomname[k-1] = pdba->atomname[k];
1386 copy_rvec(x[k], x[k-1]);
1397 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1398 t_atoms *pdba, rvec *x,
1405 for (i = 0; i < pdba->nr; i++)
1407 oldnm = *pdba->atomname[i];
1408 rptr = &restp[pdba->atom[i].resind];
1409 for (j = 0; (j < rptr->natom); j++)
1411 if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1416 if (j == rptr->natom)
1418 /* Not found yet, check if we have to rename this atom */
1419 if (match_atomnames_with_rtp_atom(pdba, x, i,
1420 rptr, &(hb[pdba->atom[i].resind]),
1423 /* We deleted this atom, decrease the atom counter by 1. */
1430 #define NUM_CMAP_ATOMS 5
1431 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t rt)
1433 int residx, i, j, k;
1435 t_resinfo *resinfo = atoms->resinfo;
1436 int nres = atoms->nres;
1438 atom_id cmap_atomid[NUM_CMAP_ATOMS];
1439 int cmap_chainnum = -1, this_residue_index;
1450 fprintf(stderr, "Making cmap torsions...");
1452 /* End loop at nres-1, since the very last residue does not have a +N atom, and
1453 * therefore we get a valgrind invalid 4 byte read error with atom am */
1454 for (residx = 0; residx < nres-1; residx++)
1456 /* Add CMAP terms from the list of CMAP interactions */
1457 for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1460 /* Loop over atoms in a candidate CMAP interaction and
1461 * check that they exist, are from the same chain and are
1462 * from residues labelled as protein. */
1463 for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1465 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
1466 i, atoms, ptr, TRUE);
1467 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1470 /* This break is necessary, because cmap_atomid[k]
1471 * == NO_ATID cannot be safely used as an index
1472 * into the atom array. */
1475 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1478 cmap_chainnum = resinfo[this_residue_index].chainnum;
1482 /* Does the residue for this atom have the same
1483 * chain number as the residues for previous
1485 bAddCMAP = bAddCMAP &&
1486 cmap_chainnum == resinfo[this_residue_index].chainnum;
1488 bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt, *(resinfo[this_residue_index].name));
1493 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);
1497 if (residx < nres-1)
1499 while (atoms->atom[i].resind < residx+1)
1506 /* Start the next residue */
1510 scrub_charge_groups(int *cgnr, int natoms)
1514 for (i = 0; i < natoms; i++)
1521 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1522 t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1523 int nrtp, t_restp rtp[],
1524 t_restp *restp, t_hackblock *hb,
1525 gmx_bool bAllowMissing,
1526 gmx_bool bVsites, gmx_bool bVsiteAromatics,
1529 int nssbonds, t_ssbond *ssbonds,
1530 real long_bond_dist, real short_bond_dist,
1531 gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1532 gmx_bool bRenumRes, gmx_bool bRTPresname)
1538 t_params plist[F_NRE];
1545 gmx_residuetype_t rt;
1548 gmx_residuetype_init(&rt);
1552 print_resall(debug, atoms->nres, restp, atype);
1553 dump_hb(debug, atoms->nres, hb);
1557 at2bonds(&(plist[F_BONDS]), hb,
1559 long_bond_dist, short_bond_dist);
1561 /* specbonds: disulphide bonds & heme-his */
1562 do_ssbonds(&(plist[F_BONDS]),
1563 atoms, nssbonds, ssbonds,
1566 nmissat = name2type(atoms, &cgnr, atype, restp, rt);
1571 fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1576 gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1581 /* Cleanup bonds (sort and rm doubles) */
1582 clean_bonds(&(plist[F_BONDS]));
1584 snew(vsite_type, atoms->nr);
1585 for (i = 0; i < atoms->nr; i++)
1587 vsite_type[i] = NOTSET;
1591 /* determine which atoms will be vsites and add dummy masses
1592 also renumber atom numbers in plist[0..F_NRE]! */
1593 do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1594 &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1597 /* Make Angles and Dihedrals */
1598 fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1599 snew(excls, atoms->nr);
1600 init_nnb(&nnb, atoms->nr, 4);
1601 gen_nnb(&nnb, plist);
1602 print_nnb(&nnb, "NNB");
1603 gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1609 gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1610 if (plist[F_CMAP].nr > 0)
1612 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1617 /* set mass of all remaining hydrogen atoms */
1620 do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1624 /* Cleanup bonds (sort and rm doubles) */
1625 /* clean_bonds(&(plist[F_BONDS]));*/
1628 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1629 " %4d pairs, %4d bonds and %4d virtual sites\n",
1630 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1631 plist[F_LJ14].nr, plist[F_BONDS].nr,
1632 plist[F_VSITE2].nr +
1633 plist[F_VSITE3].nr +
1634 plist[F_VSITE3FD].nr +
1635 plist[F_VSITE3FAD].nr +
1636 plist[F_VSITE3OUT].nr +
1637 plist[F_VSITE4FD].nr +
1638 plist[F_VSITE4FDN].nr );
1640 print_sums(atoms, FALSE);
1642 if (FALSE == bChargeGroups)
1644 scrub_charge_groups(cgnr, atoms->nr);
1649 for (i = 0; i < atoms->nres; i++)
1651 atoms->resinfo[i].nr = i + 1;
1652 atoms->resinfo[i].ic = ' ';
1658 fprintf(stderr, "Writing topology\n");
1659 /* We can copy the bonded types from the first restp,
1660 * since the types have to be identical for all residues in one molecule.
1662 for (i = 0; i < ebtsNR; i++)
1664 bts[i] = restp[0].rb[i].type;
1666 write_top(top_file, posre_fn, molname,
1668 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1672 free_t_hackblock(atoms->nres, &hb);
1673 free_t_restp(atoms->nres, &restp);
1674 gmx_residuetype_destroy(rt);
1676 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1678 for (i = 0; i < F_NRE; i++)
1680 sfree(plist[i].param);
1682 for (i = 0; i < atoms->nr; i++)