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"
71 /* this must correspond to enum in pdb2top.h */
72 const char *hh[ehisNR] = { "HISD", "HISE", "HISH", "HIS1" };
74 static int missing_atoms(t_restp *rp, int resind, t_atoms *at, int i0, int i)
78 gmx_bool bFound, bRet;
81 for (j = 0; j < rp->natom; j++)
83 name = *(rp->atomname[j]);
85 for (k = i0; k < i; k++)
87 bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]), name));
92 fprintf(stderr, "\nWARNING: "
93 "atom %s is missing in residue %s %d in the pdb file\n",
94 name, *(at->resinfo[resind].name), at->resinfo[resind].nr);
95 if (name[0] == 'H' || name[0] == 'h')
97 fprintf(stderr, " You might need to add atom %s to the hydrogen database of building block %s\n"
98 " in the file %s.hdb (see the manual)\n",
99 name, *(at->resinfo[resind].rtp), rp->filebase);
101 fprintf(stderr, "\n");
108 gmx_bool is_int(double x)
110 const double tol = 1e-4;
119 return (fabs(x-ix) < tol);
122 static void swap_strings(char **s, int i, int j)
132 choose_ff(const char *ffsel,
133 char *forcefield, int ff_maxlen,
134 char *ffdir, int ffdir_maxlen)
137 char **ffdirs, **ffs, **ffs_dir, *ptr;
138 int i, j, sel, cwdsel, nfound;
139 char buf[STRLEN], **desc;
143 nff = fflib_search_file_in_dirend(fflib_forcefield_itp(),
144 fflib_forcefield_dir_ext(),
149 gmx_fatal(FARGS, "No force fields found (files with name '%s' in subdirectories ending on '%s')",
150 fflib_forcefield_itp(), fflib_forcefield_dir_ext());
153 /* Replace with unix path separators */
154 if (DIR_SEPARATOR != '/')
156 for (i = 0; i < nff; i++)
158 while ( (ptr = strchr(ffdirs[i], DIR_SEPARATOR)) != NULL)
165 /* Store the force field names in ffs */
168 for (i = 0; i < nff; i++)
170 /* Remove the path from the ffdir name - use our unix standard here! */
171 ptr = strrchr(ffdirs[i], '/');
174 ffs[i] = strdup(ffdirs[i]);
175 ffs_dir[i] = low_gmxlibfn(ffdirs[i], FALSE, FALSE);
176 if (ffs_dir[i] == NULL)
178 gmx_fatal(FARGS, "Can no longer find file '%s'", ffdirs[i]);
183 ffs[i] = strdup(ptr+1);
184 ffs_dir[i] = strdup(ffdirs[i]);
186 ffs_dir[i][strlen(ffs_dir[i])-strlen(ffs[i])-1] = '\0';
187 /* Remove the extension from the ffdir name */
188 ffs[i][strlen(ffs[i])-strlen(fflib_forcefield_dir_ext())] = '\0';
196 for (i = 0; i < nff; i++)
198 if (strcmp(ffs[i], ffsel) == 0)
200 /* Matching ff name */
204 if (strncmp(ffs_dir[i], ".", 1) == 0)
221 "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
222 "current directory. Use interactive selection (not the -ff option) if\n"
223 "you would prefer a different one.\n", ffsel, nfound);
228 "Force field '%s' occurs in %d places, but not in the current directory.\n"
229 "Run without the -ff switch and select the force field interactively.", ffsel, nfound);
232 else if (nfound == 0)
234 gmx_fatal(FARGS, "Could not find force field '%s' in current directory, install tree or GMXDATA path.", ffsel);
240 for (i = 0; (i < nff); i++)
242 sprintf(buf, "%s%c%s%s%c%s",
243 ffs_dir[i], DIR_SEPARATOR,
244 ffs[i], fflib_forcefield_dir_ext(), DIR_SEPARATOR,
245 fflib_forcefield_doc());
248 /* We don't use fflib_open, because we don't want printf's */
249 fp = ffopen(buf, "r");
250 snew(desc[i], STRLEN);
251 get_a_line(fp, desc[i], STRLEN);
256 desc[i] = strdup(ffs[i]);
259 /* Order force fields from the same dir alphabetically
260 * and put deprecated force fields at the end.
262 for (i = 0; (i < nff); i++)
264 for (j = i+1; (j < nff); j++)
266 if (strcmp(ffs_dir[i], ffs_dir[j]) == 0 &&
267 ((desc[i][0] == '[' && desc[j][0] != '[') ||
268 ((desc[i][0] == '[' || desc[j][0] != '[') &&
269 gmx_strcasecmp(desc[i], desc[j]) > 0)))
271 swap_strings(ffdirs, i, j);
272 swap_strings(ffs, i, j);
273 swap_strings(desc, i, j);
278 printf("\nSelect the Force Field:\n");
279 for (i = 0; (i < nff); i++)
281 if (i == 0 || strcmp(ffs_dir[i-1], ffs_dir[i]) != 0)
283 if (strcmp(ffs_dir[i], ".") == 0)
285 printf("From current directory:\n");
289 printf("From '%s':\n", ffs_dir[i]);
292 printf("%2d: %s\n", i+1, desc[i]);
299 pret = fgets(buf, STRLEN, stdin);
303 sscanf(buf, "%d", &sel);
307 while (pret == NULL || (sel < 0) || (sel >= nff));
309 /* Check for a current limitation of the fflib code.
310 * It will always read from the first ff directory in the list.
311 * This check assumes that the order of ffs matches the order
312 * in which fflib_open searches ff library files.
314 for (i = 0; i < sel; i++)
316 if (strcmp(ffs[i], ffs[sel]) == 0)
318 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.",
319 ffs[sel], fflib_forcefield_dir_ext());
328 if (strlen(ffs[sel]) >= (size_t)ff_maxlen)
330 gmx_fatal(FARGS, "Length of force field name (%d) >= maxlen (%d)",
331 strlen(ffs[sel]), ff_maxlen);
333 strcpy(forcefield, ffs[sel]);
335 if (strlen(ffdirs[sel]) >= (size_t)ffdir_maxlen)
337 gmx_fatal(FARGS, "Length of force field dir (%d) >= maxlen (%d)",
338 strlen(ffdirs[sel]), ffdir_maxlen);
340 strcpy(ffdir, ffdirs[sel]);
342 for (i = 0; (i < nff); i++)
353 void choose_watermodel(const char *wmsel, const char *ffdir,
356 const char *fn_watermodels = "watermodels.dat";
357 char fn_list[STRLEN];
364 if (strcmp(wmsel, "none") == 0)
370 else if (strcmp(wmsel, "select") != 0)
372 *watermodel = strdup(wmsel);
377 sprintf(fn_list, "%s%c%s", ffdir, DIR_SEPARATOR, fn_watermodels);
379 if (!fflib_fexist(fn_list))
381 fprintf(stderr, "No file '%s' found, will not include a water model\n",
388 fp = fflib_open(fn_list);
389 printf("\nSelect the Water Model:\n");
392 while (get_a_line(fp, buf, STRLEN))
394 srenew(model, nwm+1);
395 snew(model[nwm], STRLEN);
396 sscanf(buf, "%s%n", model[nwm], &i);
400 fprintf(stderr, "%2d: %s\n", nwm+1, buf+i);
409 fprintf(stderr, "%2d: %s\n", nwm+1, "None");
413 pret = fgets(buf, STRLEN, stdin);
417 sscanf(buf, "%d", &sel);
421 while (pret == NULL || sel < 0 || sel > nwm);
429 *watermodel = strdup(model[sel]);
432 for (i = 0; i < nwm; i++)
439 static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
440 t_restp restp[], gmx_residuetype_t rt)
442 int i, j, prevresind, resind, i0, prevcg, cg, curcg;
444 gmx_bool bProt, bNterm;
461 for (i = 0; (i < at->nr); i++)
464 if (at->atom[i].resind != resind)
466 resind = at->atom[i].resind;
467 bProt = gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name));
468 bNterm = bProt && (resind == 0);
471 nmissat += missing_atoms(&restp[prevresind], prevresind, at, i0, i);
475 if (at->atom[i].m == 0)
479 fprintf(debug, "atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
480 i+1, *(at->atomname[i]), curcg, prevcg,
481 j == NOTSET ? NOTSET : restp[resind].cgnr[j]);
485 name = *(at->atomname[i]);
486 j = search_jtype(&restp[resind], name, bNterm);
487 at->atom[i].type = restp[resind].atom[j].type;
488 at->atom[i].q = restp[resind].atom[j].q;
489 at->atom[i].m = get_atomtype_massA(restp[resind].atom[j].type,
491 cg = restp[resind].cgnr[j];
492 /* A charge group number -1 signals a separate charge group
495 if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) )
504 fprintf(debug, "atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
505 i+1, *(at->atomname[i]), curcg, qt, is_int(qt));
516 at->atom[i].typeB = at->atom[i].type;
517 at->atom[i].qB = at->atom[i].q;
518 at->atom[i].mB = at->atom[i].m;
520 nmissat += missing_atoms(&restp[resind], resind, at, i0, i);
525 static void print_top_heavy_H(FILE *out, real mHmult)
529 fprintf(out, "; Using deuterium instead of hydrogen\n\n");
531 else if (mHmult == 4.0)
533 fprintf(out, "#define HEAVY_H\n\n");
535 else if (mHmult != 1.0)
537 fprintf(stderr, "WARNING: unsupported proton mass multiplier (%g) "
538 "in pdb2top\n", mHmult);
542 void print_top_comment(FILE *out,
543 const char *filename,
544 const char *generator,
549 char ffdir_parent[STRLEN];
552 nice_header(out, filename);
553 fprintf(out, ";\tThis is a %s topology file\n;\n", bITP ? "include" : "standalone");
554 fprintf(out, ";\tIt was generated using program:\n;\t%s\n;\n",
555 (NULL == generator) ? "unknown" : generator);
556 fprintf(out, ";\tCommand line was:\n;\t%s\n;\n", command_line());
558 if (strchr(ffdir, '/') == NULL)
560 fprintf(out, ";\tForce field was read from the standard Gromacs share directory.\n;\n\n");
562 else if (ffdir[0] == '.')
564 fprintf(out, ";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
568 strncpy(ffdir_parent, ffdir, STRLEN-1);
569 ffdir_parent[STRLEN-1] = '\0'; /*make sure it is 0-terminated even for long string*/
570 p = strrchr(ffdir_parent, '/');
575 ";\tForce field data was read from:\n"
579 ";\tThis might be a non-standard force field location. When you use this topology, the\n"
580 ";\tforce field must either be present in the current directory, or the location\n"
581 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
586 void print_top_header(FILE *out, const char *filename,
587 const char *title, gmx_bool bITP, const char *ffdir, real mHmult)
591 print_top_comment(out, filename, title, ffdir, bITP);
593 print_top_heavy_H(out, mHmult);
594 fprintf(out, "; Include forcefield parameters\n");
596 p = strrchr(ffdir, '/');
597 p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
599 fprintf(out, "#include \"%s/%s\"\n\n", p, fflib_forcefield_itp());
602 static void print_top_posre(FILE *out, const char *pr)
604 fprintf(out, "; Include Position restraint file\n");
605 fprintf(out, "#ifdef POSRES\n");
606 fprintf(out, "#include \"%s\"\n", pr);
607 fprintf(out, "#endif\n\n");
610 static void print_top_water(FILE *out, const char *ffdir, const char *water)
615 fprintf(out, "; Include water topology\n");
617 p = strrchr(ffdir, '/');
618 p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
619 fprintf(out, "#include \"%s/%s.itp\"\n", p, water);
622 fprintf(out, "#ifdef POSRES_WATER\n");
623 fprintf(out, "; Position restraint for each water oxygen\n");
624 fprintf(out, "[ position_restraints ]\n");
625 fprintf(out, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
626 fprintf(out, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
627 fprintf(out, "#endif\n");
630 sprintf(buf, "%s/ions.itp", p);
632 if (fflib_fexist(buf))
634 fprintf(out, "; Include topology for ions\n");
635 fprintf(out, "#include \"%s\"\n", buf);
640 static void print_top_system(FILE *out, const char *title)
642 fprintf(out, "[ %s ]\n", dir2str(d_system));
643 fprintf(out, "; Name\n");
644 fprintf(out, "%s\n\n", title[0] ? title : "Protein");
647 void print_top_mols(FILE *out,
648 const char *title, const char *ffdir, const char *water,
649 int nincl, char **incls, int nmol, t_mols *mols)
656 fprintf(out, "; Include chain topologies\n");
657 for (i = 0; (i < nincl); i++)
659 incl = strrchr(incls[i], DIR_SEPARATOR);
666 /* Remove the path from the include name */
669 fprintf(out, "#include \"%s\"\n", incl);
676 print_top_water(out, ffdir, water);
678 print_top_system(out, title);
682 fprintf(out, "[ %s ]\n", dir2str(d_molecules));
683 fprintf(out, "; %-15s %5s\n", "Compound", "#mols");
684 for (i = 0; (i < nmol); i++)
686 fprintf(out, "%-15s %5d\n", mols[i].name, mols[i].nr);
691 void write_top(FILE *out, char *pr, char *molname,
692 t_atoms *at, gmx_bool bRTPresname,
693 int bts[], t_params plist[], t_excls excls[],
694 gpp_atomtype_t atype, int *cgnr, int nrexcl)
695 /* NOTE: nrexcl is not the size of *excl! */
697 if (at && atype && cgnr)
699 fprintf(out, "[ %s ]\n", dir2str(d_moleculetype));
700 fprintf(out, "; %-15s %5s\n", "Name", "nrexcl");
701 fprintf(out, "%-15s %5d\n\n", molname ? molname : "Protein", nrexcl);
703 print_atoms(out, atype, at, cgnr, bRTPresname);
704 print_bondeds(out, at->nr, d_bonds, F_BONDS, bts[ebtsBONDS], plist);
705 print_bondeds(out, at->nr, d_constraints, F_CONSTR, 0, plist);
706 print_bondeds(out, at->nr, d_constraints, F_CONSTRNC, 0, plist);
707 print_bondeds(out, at->nr, d_pairs, F_LJ14, 0, plist);
708 print_excl(out, at->nr, excls);
709 print_bondeds(out, at->nr, d_angles, F_ANGLES, bts[ebtsANGLES], plist);
710 print_bondeds(out, at->nr, d_dihedrals, F_PDIHS, bts[ebtsPDIHS], plist);
711 print_bondeds(out, at->nr, d_dihedrals, F_IDIHS, bts[ebtsIDIHS], plist);
712 print_bondeds(out, at->nr, d_cmap, F_CMAP, bts[ebtsCMAP], plist);
713 print_bondeds(out, at->nr, d_polarization, F_POLARIZATION, 0, plist);
714 print_bondeds(out, at->nr, d_thole_polarization, F_THOLE_POL, 0, plist);
715 print_bondeds(out, at->nr, d_vsites2, F_VSITE2, 0, plist);
716 print_bondeds(out, at->nr, d_vsites3, F_VSITE3, 0, plist);
717 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FD, 0, plist);
718 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FAD, 0, plist);
719 print_bondeds(out, at->nr, d_vsites3, F_VSITE3OUT, 0, plist);
720 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FD, 0, plist);
721 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FDN, 0, plist);
725 print_top_posre(out, pr);
730 static atom_id search_res_atom(const char *type, int resind,
732 const char *bondtype, gmx_bool bAllowMissing)
736 for (i = 0; (i < atoms->nr); i++)
738 if (atoms->atom[i].resind == resind)
740 return search_atom(type, i, atoms, bondtype, bAllowMissing);
747 static void do_ssbonds(t_params *ps, t_atoms *atoms,
748 int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
753 for (i = 0; (i < nssbonds); i++)
755 ri = ssbonds[i].res1;
756 rj = ssbonds[i].res2;
757 ai = search_res_atom(ssbonds[i].a1, ri, atoms,
758 "special bond", bAllowMissing);
759 aj = search_res_atom(ssbonds[i].a2, rj, atoms,
760 "special bond", bAllowMissing);
761 if ((ai == NO_ATID) || (aj == NO_ATID))
763 gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
764 ssbonds[i].a1, ssbonds[i].a2);
766 add_param(ps, ai, aj, NULL, NULL);
770 static void at2bonds(t_params *psb, t_hackblock *hb,
773 real long_bond_dist, real short_bond_dist)
777 real dist2, long_bond_dist2, short_bond_dist2;
780 long_bond_dist2 = sqr(long_bond_dist);
781 short_bond_dist2 = sqr(short_bond_dist);
792 fprintf(stderr, "Making bonds...\n");
794 for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
796 /* add bonds from list of bonded interactions */
797 for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
799 /* Unfortunately we can not issue errors or warnings
800 * for missing atoms in bonds, as the hydrogens and terminal atoms
801 * have not been added yet.
803 ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].AI, i, atoms,
805 aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].AJ, i, atoms,
807 if (ai != NO_ATID && aj != NO_ATID)
809 dist2 = distance2(x[ai], x[aj]);
810 if (dist2 > long_bond_dist2)
812 fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
813 ai+1, aj+1, sqrt(dist2));
815 else if (dist2 < short_bond_dist2)
817 fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
818 ai+1, aj+1, sqrt(dist2));
820 add_param(psb, ai, aj, NULL, hb[resind].rb[ebtsBONDS].b[j].s);
823 /* add bonds from list of hacks (each added atom gets a bond) */
824 while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
826 for (j = 0; j < hb[resind].nhack; j++)
828 if ( ( hb[resind].hack[j].tp > 0 ||
829 hb[resind].hack[j].oname == NULL ) &&
830 strcmp(hb[resind].hack[j].AI, *(atoms->atomname[i])) == 0)
832 switch (hb[resind].hack[j].tp)
834 case 9: /* COOH terminus */
835 add_param(psb, i, i+1, NULL, NULL); /* C-O */
836 add_param(psb, i, i+2, NULL, NULL); /* C-OA */
837 add_param(psb, i+2, i+3, NULL, NULL); /* OA-H */
840 for (k = 0; (k < hb[resind].hack[j].nr); k++)
842 add_param(psb, i, i+k+1, NULL, NULL);
849 /* we're now at the start of the next residue */
853 static int pcompar(const void *a, const void *b)
867 return strlen(pb->s) - strlen(pa->s);
875 static void clean_bonds(t_params *ps)
882 /* swap atomnumbers in bond if first larger than second: */
883 for (i = 0; (i < ps->nr); i++)
885 if (ps->param[i].AJ < ps->param[i].AI)
888 ps->param[i].AI = ps->param[i].AJ;
894 qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
896 /* remove doubles, keep the first one always. */
898 for (i = 1; (i < ps->nr); i++)
900 if ((ps->param[i].AI != ps->param[j-1].AI) ||
901 (ps->param[i].AJ != ps->param[j-1].AJ) )
905 cp_param(&(ps->param[j]), &(ps->param[i]));
910 fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
915 fprintf(stderr, "No bonds\n");
919 void print_sums(t_atoms *atoms, gmx_bool bSystem)
927 where = " in system";
936 for (i = 0; (i < atoms->nr); i++)
938 m += atoms->atom[i].m;
939 qtot += atoms->atom[i].q;
942 fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
943 fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
946 static void check_restp_type(const char *name, int t1, int t2)
950 gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
954 static void check_restp_types(t_restp *r0, t_restp *r1)
958 check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
959 check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
960 check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
961 check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
963 for (i = 0; i < ebtsNR; i++)
965 check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
969 void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
973 const char *Hnum = "123456";
977 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
979 * restp->atomname[at_start], resnr, restp->resname);
981 strcpy(buf, hack->nname);
982 buf[strlen(buf)+1] = '\0';
985 buf[strlen(buf)] = '-';
988 restp->natom += hack->nr;
989 srenew(restp->atom, restp->natom);
990 srenew(restp->atomname, restp->natom);
991 srenew(restp->cgnr, restp->natom);
993 for (k = restp->natom-1; k > at_start+hack->nr; k--)
996 restp->atom [k - hack->nr];
998 restp->atomname[k - hack->nr];
1000 restp->cgnr [k - hack->nr];
1003 for (k = 0; k < hack->nr; k++)
1005 /* set counter in atomname */
1008 buf[strlen(buf)-1] = Hnum[k];
1010 snew( restp->atomname[at_start+1+k], 1);
1011 restp->atom [at_start+1+k] = *hack->atom;
1012 *restp->atomname[at_start+1+k] = strdup(buf);
1013 if (hack->cgnr != NOTSET)
1015 restp->cgnr [at_start+1+k] = hack->cgnr;
1019 restp->cgnr [at_start+1+k] = restp->cgnr[at_start];
1024 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1025 int nrtp, t_restp rtp[],
1026 int nres, t_resinfo *resinfo,
1028 t_hackblock **ntdb, t_hackblock **ctdb,
1035 const char *Hnum = "123456";
1037 gmx_bool bN, bC, bRM;
1041 /* first the termini */
1042 for (i = 0; i < nterpairs; i++)
1044 if (rn[i] >= 0 && ntdb[i] != NULL)
1046 copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1048 if (rc[i] >= 0 && ctdb[i] != NULL)
1050 merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1054 /* then the whole rtp */
1055 for (i = 0; i < nres; i++)
1057 /* Here we allow a mismatch of one character when looking for the rtp entry.
1058 * For such a mismatch there should be only one mismatching name.
1059 * This is mainly useful for small molecules such as ions.
1060 * Note that this will usually not work for protein, DNA and RNA,
1061 * since there the residue names should be listed in residuetypes.dat
1062 * and an error will have been generated earlier in the process.
1064 key = *resinfo[i].rtp;
1065 snew(resinfo[i].rtp, 1);
1066 *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1067 res = get_restp(*resinfo[i].rtp, nrtp, rtp);
1068 copy_t_restp(res, &(*restp)[i]);
1070 /* Check that we do not have different bonded types in one molecule */
1071 check_restp_types(&(*restp)[0], &(*restp)[i]);
1074 for (j = 0; j < nterpairs && tern == -1; j++)
1082 for (j = 0; j < nterpairs && terc == -1; j++)
1089 bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1091 if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1092 (terc >= 0 && ctdb[terc] == NULL)))
1094 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).");
1096 if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1097 (terc >= 0 && ctdb[terc]->nhack == 0)))
1099 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.");
1103 /* now perform t_hack's on t_restp's,
1104 i.e. add's and deletes from termini database will be
1105 added to/removed from residue topology
1106 we'll do this on one big dirty loop, so it won't make easy reading! */
1107 for (i = 0; i < nres; i++)
1109 for (j = 0; j < (*hb)[i].nhack; j++)
1111 if ( (*hb)[i].hack[j].nr)
1113 /* find atom in restp */
1114 for (l = 0; l < (*restp)[i].natom; l++)
1116 if ( ( (*hb)[i].hack[j].oname == NULL &&
1117 strcmp((*hb)[i].hack[j].AI, *(*restp)[i].atomname[l]) == 0 ) ||
1118 ( (*hb)[i].hack[j].oname != NULL &&
1119 strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1124 if (l == (*restp)[i].natom)
1126 /* If we are doing an atom rename only, we don't need
1127 * to generate a fatal error if the old name is not found
1130 /* Deleting can happen also only on the input atoms,
1131 * not necessarily always on the rtp entry.
1133 if (!((*hb)[i].hack[j].oname != NULL &&
1134 (*hb)[i].hack[j].nname != NULL) &&
1135 !((*hb)[i].hack[j].oname != NULL &&
1136 (*hb)[i].hack[j].nname == NULL))
1139 "atom %s not found in buiding block %d%s "
1140 "while combining tdb and rtp",
1141 (*hb)[i].hack[j].oname != NULL ?
1142 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].AI,
1143 i+1, *resinfo[i].rtp);
1148 if ( (*hb)[i].hack[j].oname == NULL)
1151 add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1156 if ( (*hb)[i].hack[j].nname == NULL)
1158 /* we're deleting */
1161 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1162 *(*restp)[i].atomname[l],
1163 i+1, (*restp)[i].resname);
1165 /* shift the rest */
1166 (*restp)[i].natom--;
1167 for (k = l; k < (*restp)[i].natom; k++)
1169 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
1170 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1171 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
1173 /* give back space */
1174 srenew((*restp)[i].atom, (*restp)[i].natom);
1175 srenew((*restp)[i].atomname, (*restp)[i].natom);
1176 srenew((*restp)[i].cgnr, (*restp)[i].natom);
1178 else /* nname != NULL */
1179 { /* we're replacing */
1182 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1183 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1184 i+1, (*restp)[i].resname);
1186 snew( (*restp)[i].atomname[l], 1);
1187 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
1188 *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1189 if ( (*hb)[i].hack[j].cgnr != NOTSET)
1191 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
1201 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1208 return (gmx_strcasecmp(anm, hack->nname) == 0);
1212 if (isdigit(anm[strlen(anm)-1]))
1214 *nr = anm[strlen(anm)-1] - '0';
1220 if (*nr <= 0 || *nr > hack->nr)
1226 return (strlen(anm) == strlen(hack->nname) + 1 &&
1227 gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1232 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1233 t_restp *rptr, t_hackblock *hbr,
1238 char *oldnm, *newnm;
1240 char *start_at, buf[STRLEN];
1242 gmx_bool bReplaceReplace, bFoundInAdd;
1245 oldnm = *pdba->atomname[atind];
1246 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1249 for (j = 0; j < hbr->nhack; j++)
1251 if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1252 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1254 /* This is a replace entry. */
1255 /* Check if we are not replacing a replaced atom. */
1256 bReplaceReplace = FALSE;
1257 for (k = 0; k < hbr->nhack; k++)
1260 hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1261 gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1263 /* The replace in hack[j] replaces an atom that
1264 * was already replaced in hack[k], we do not want
1265 * second or higher level replaces at this stage.
1267 bReplaceReplace = TRUE;
1270 if (bReplaceReplace)
1272 /* Skip this replace. */
1276 /* This atom still has the old name, rename it */
1277 newnm = hbr->hack[j].nname;
1278 for (k = 0; k < rptr->natom; k++)
1280 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1285 if (k == rptr->natom)
1287 /* The new name is not present in the rtp.
1288 * We need to apply the replace also to the rtp entry.
1291 /* We need to find the add hack that can add this atom
1292 * to find out after which atom it should be added.
1294 bFoundInAdd = FALSE;
1295 for (k = 0; k < hbr->nhack; k++)
1297 if (hbr->hack[k].oname == NULL &&
1298 hbr->hack[k].nname != NULL &&
1299 atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1303 start_at = hbr->hack[k].a[0];
1307 sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1310 for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1312 if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1317 if (start_nr == rptr->natom)
1319 gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1320 start_at, rptr->resname, newnm);
1322 /* We can add the atom after atom start_nr */
1323 add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1331 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'",
1333 hbr->hack[j].oname, hbr->hack[j].nname,
1340 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1341 oldnm, rptr->resname, resnr, newnm);
1343 /* Rename the atom in pdba */
1344 snew(pdba->atomname[atind], 1);
1345 *pdba->atomname[atind] = strdup(newnm);
1347 else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1348 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1350 /* This is a delete entry, check if this atom is present
1351 * in the rtp entry of this residue.
1353 for (k = 0; k < rptr->natom; k++)
1355 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1360 if (k == rptr->natom)
1362 /* This atom is not present in the rtp entry,
1363 * delete is now from the input pdba.
1367 printf("Deleting atom '%s' in residue '%s' %d\n",
1368 oldnm, rptr->resname, resnr);
1370 /* We should free the atom name,
1371 * but it might be used multiple times in the symtab.
1372 * sfree(pdba->atomname[atind]);
1374 for (k = atind+1; k < pdba->nr; k++)
1376 pdba->atom[k-1] = pdba->atom[k];
1377 pdba->atomname[k-1] = pdba->atomname[k];
1378 copy_rvec(x[k], x[k-1]);
1389 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1390 t_atoms *pdba, rvec *x,
1394 char *oldnm, *newnm;
1399 char *start_at, buf[STRLEN];
1401 gmx_bool bFoundInAdd;
1403 for (i = 0; i < pdba->nr; i++)
1405 oldnm = *pdba->atomname[i];
1406 resnr = pdba->resinfo[pdba->atom[i].resind].nr;
1407 rptr = &restp[pdba->atom[i].resind];
1408 for (j = 0; (j < rptr->natom); j++)
1410 if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1415 if (j == rptr->natom)
1417 /* Not found yet, check if we have to rename this atom */
1418 if (match_atomnames_with_rtp_atom(pdba, x, i,
1419 rptr, &(hb[pdba->atom[i].resind]),
1422 /* We deleted this atom, decrease the atom counter by 1. */
1429 #define NUM_CMAP_ATOMS 5
1430 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t rt)
1432 int residx, i, j, k;
1434 t_resinfo *resinfo = atoms->resinfo;
1435 int nres = atoms->nres;
1437 atom_id cmap_atomid[NUM_CMAP_ATOMS];
1438 int cmap_chainnum = -1, this_residue_index;
1449 fprintf(stderr, "Making cmap torsions...");
1451 /* End loop at nres-1, since the very last residue does not have a +N atom, and
1452 * therefore we get a valgrind invalid 4 byte read error with atom am */
1453 for (residx = 0; residx < nres-1; residx++)
1455 /* Add CMAP terms from the list of CMAP interactions */
1456 for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1459 /* Loop over atoms in a candidate CMAP interaction and
1460 * check that they exist, are from the same chain and are
1461 * from residues labelled as protein. */
1462 for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1464 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
1465 i, atoms, ptr, TRUE);
1466 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1469 /* This break is necessary, because cmap_atomid[k]
1470 * == NO_ATID cannot be safely used as an index
1471 * into the atom array. */
1474 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1477 cmap_chainnum = resinfo[this_residue_index].chainnum;
1481 /* Does the residue for this atom have the same
1482 * chain number as the residues for previous
1484 bAddCMAP = bAddCMAP &&
1485 cmap_chainnum == resinfo[this_residue_index].chainnum;
1487 bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt, *(resinfo[this_residue_index].name));
1492 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);
1496 if (residx < nres-1)
1498 while (atoms->atom[i].resind < residx+1)
1505 /* Start the next residue */
1509 scrub_charge_groups(int *cgnr, int natoms)
1513 for (i = 0; i < natoms; i++)
1520 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1521 t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1522 int nrtp, t_restp rtp[],
1523 t_restp *restp, t_hackblock *hb,
1524 gmx_bool bAllowMissing,
1525 gmx_bool bVsites, gmx_bool bVsiteAromatics,
1528 int nssbonds, t_ssbond *ssbonds,
1529 real long_bond_dist, real short_bond_dist,
1530 gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1531 gmx_bool bRenumRes, gmx_bool bRTPresname)
1537 t_params plist[F_NRE];
1544 gmx_residuetype_t rt;
1547 gmx_residuetype_init(&rt);
1551 print_resall(debug, atoms->nres, restp, atype);
1552 dump_hb(debug, atoms->nres, hb);
1556 at2bonds(&(plist[F_BONDS]), hb,
1558 long_bond_dist, short_bond_dist);
1560 /* specbonds: disulphide bonds & heme-his */
1561 do_ssbonds(&(plist[F_BONDS]),
1562 atoms, nssbonds, ssbonds,
1565 nmissat = name2type(atoms, &cgnr, atype, restp, rt);
1570 fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1575 gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1580 /* Cleanup bonds (sort and rm doubles) */
1581 clean_bonds(&(plist[F_BONDS]));
1583 snew(vsite_type, atoms->nr);
1584 for (i = 0; i < atoms->nr; i++)
1586 vsite_type[i] = NOTSET;
1590 /* determine which atoms will be vsites and add dummy masses
1591 also renumber atom numbers in plist[0..F_NRE]! */
1592 do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1593 &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1596 /* Make Angles and Dihedrals */
1597 fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1598 snew(excls, atoms->nr);
1599 init_nnb(&nnb, atoms->nr, 4);
1600 gen_nnb(&nnb, plist);
1601 print_nnb(&nnb, "NNB");
1602 gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1608 gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1609 if (plist[F_CMAP].nr > 0)
1611 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1616 /* set mass of all remaining hydrogen atoms */
1619 do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1623 /* Cleanup bonds (sort and rm doubles) */
1624 /* clean_bonds(&(plist[F_BONDS]));*/
1627 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1628 " %4d pairs, %4d bonds and %4d virtual sites\n",
1629 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1630 plist[F_LJ14].nr, plist[F_BONDS].nr,
1631 plist[F_VSITE2].nr +
1632 plist[F_VSITE3].nr +
1633 plist[F_VSITE3FD].nr +
1634 plist[F_VSITE3FAD].nr +
1635 plist[F_VSITE3OUT].nr +
1636 plist[F_VSITE4FD].nr +
1637 plist[F_VSITE4FDN].nr );
1639 print_sums(atoms, FALSE);
1641 if (FALSE == bChargeGroups)
1643 scrub_charge_groups(cgnr, atoms->nr);
1648 for (i = 0; i < atoms->nres; i++)
1650 atoms->resinfo[i].nr = i + 1;
1651 atoms->resinfo[i].ic = ' ';
1657 fprintf(stderr, "Writing topology\n");
1658 /* We can copy the bonded types from the first restp,
1659 * since the types have to be identical for all residues in one molecule.
1661 for (i = 0; i < ebtsNR; i++)
1663 bts[i] = restp[0].rb[i].type;
1665 write_top(top_file, posre_fn, molname,
1667 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1671 free_t_hackblock(atoms->nres, &hb);
1672 free_t_restp(atoms->nres, &restp);
1673 gmx_residuetype_destroy(rt);
1675 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1677 for (i = 0; i < F_NRE; i++)
1679 sfree(plist[i].param);
1681 for (i = 0; i < atoms->nr; i++)