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)
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]);
300 pret = fgets(buf, STRLEN, stdin);
304 sel = strtol(buf, NULL, 10);
308 while (pret == NULL || (sel < 0) || (sel >= nff));
310 /* Check for a current limitation of the fflib code.
311 * It will always read from the first ff directory in the list.
312 * This check assumes that the order of ffs matches the order
313 * in which fflib_open searches ff library files.
315 for (i = 0; i < sel; i++)
317 if (strcmp(ffs[i], ffs[sel]) == 0)
319 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.",
320 ffs[sel], fflib_forcefield_dir_ext());
329 if (strlen(ffs[sel]) >= (size_t)ff_maxlen)
331 gmx_fatal(FARGS, "Length of force field name (%d) >= maxlen (%d)",
332 strlen(ffs[sel]), ff_maxlen);
334 strcpy(forcefield, ffs[sel]);
336 if (strlen(ffdirs[sel]) >= (size_t)ffdir_maxlen)
338 gmx_fatal(FARGS, "Length of force field dir (%d) >= maxlen (%d)",
339 strlen(ffdirs[sel]), ffdir_maxlen);
341 strcpy(ffdir, ffdirs[sel]);
343 for (i = 0; (i < nff); i++)
354 void choose_watermodel(const char *wmsel, const char *ffdir,
357 const char *fn_watermodels = "watermodels.dat";
358 char fn_list[STRLEN];
365 if (strcmp(wmsel, "none") == 0)
371 else if (strcmp(wmsel, "select") != 0)
373 *watermodel = strdup(wmsel);
378 sprintf(fn_list, "%s%c%s", ffdir, DIR_SEPARATOR, fn_watermodels);
380 if (!fflib_fexist(fn_list))
382 fprintf(stderr, "No file '%s' found, will not include a water model\n",
389 fp = fflib_open(fn_list);
390 printf("\nSelect the Water Model:\n");
393 while (get_a_line(fp, buf, STRLEN))
395 srenew(model, nwm+1);
396 snew(model[nwm], STRLEN);
397 sscanf(buf, "%s%n", model[nwm], &i);
401 fprintf(stderr, "%2d: %s\n", nwm+1, buf+i);
410 fprintf(stderr, "%2d: %s\n", nwm+1, "None");
415 pret = fgets(buf, STRLEN, stdin);
419 sel = strtol(buf, NULL, 10);
423 while (pret == NULL || sel < 0 || sel > nwm);
431 *watermodel = strdup(model[sel]);
434 for (i = 0; i < nwm; i++)
441 static int name2type(t_atoms *at, int **cgnr, gpp_atomtype_t atype,
442 t_restp restp[], gmx_residuetype_t rt)
444 int i, j, prevresind, resind, i0, prevcg, cg, curcg;
462 for (i = 0; (i < at->nr); i++)
465 if (at->atom[i].resind != resind)
468 resind = at->atom[i].resind;
469 bProt = gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name));
470 bNterm = bProt && (resind == 0);
473 nmissat += missing_atoms(&restp[prevresind], prevresind, at, i0, i);
477 if (at->atom[i].m == 0)
481 fprintf(debug, "atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
482 i+1, *(at->atomname[i]), curcg, prevcg,
483 j == NOTSET ? NOTSET : restp[resind].cgnr[j]);
487 name = *(at->atomname[i]);
488 j = search_jtype(&restp[resind], name, bNterm);
489 at->atom[i].type = restp[resind].atom[j].type;
490 at->atom[i].q = restp[resind].atom[j].q;
491 at->atom[i].m = get_atomtype_massA(restp[resind].atom[j].type,
493 cg = restp[resind].cgnr[j];
494 /* A charge group number -1 signals a separate charge group
497 if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) )
506 fprintf(debug, "atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
507 i+1, *(at->atomname[i]), curcg, qt, is_int(qt));
518 at->atom[i].typeB = at->atom[i].type;
519 at->atom[i].qB = at->atom[i].q;
520 at->atom[i].mB = at->atom[i].m;
522 nmissat += missing_atoms(&restp[resind], resind, at, i0, i);
527 static void print_top_heavy_H(FILE *out, real mHmult)
531 fprintf(out, "; Using deuterium instead of hydrogen\n\n");
533 else if (mHmult == 4.0)
535 fprintf(out, "#define HEAVY_H\n\n");
537 else if (mHmult != 1.0)
539 fprintf(stderr, "WARNING: unsupported proton mass multiplier (%g) "
540 "in pdb2top\n", mHmult);
544 void print_top_comment(FILE *out,
545 const char *filename,
546 const char *generator,
550 char ffdir_parent[STRLEN];
553 nice_header(out, filename);
554 fprintf(out, ";\tThis is a %s topology file\n;\n", bITP ? "include" : "standalone");
555 fprintf(out, ";\tIt was generated using program:\n;\t%s\n;\n",
556 (NULL == generator) ? "unknown" : generator);
557 fprintf(out, ";\tCommand line was:\n;\t%s\n;\n", command_line());
559 if (strchr(ffdir, '/') == NULL)
561 fprintf(out, ";\tForce field was read from the standard Gromacs share directory.\n;\n\n");
563 else if (ffdir[0] == '.')
565 fprintf(out, ";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
569 strncpy(ffdir_parent, ffdir, STRLEN-1);
570 ffdir_parent[STRLEN-1] = '\0'; /*make sure it is 0-terminated even for long string*/
571 p = strrchr(ffdir_parent, '/');
576 ";\tForce field data was read from:\n"
580 ";\tThis might be a non-standard force field location. When you use this topology, the\n"
581 ";\tforce field must either be present in the current directory, or the location\n"
582 ";\tspecified in the GMXLIB path variable or with the 'include' mdp file option.\n;\n\n",
587 void print_top_header(FILE *out, const char *filename,
588 const char *title, gmx_bool bITP, const char *ffdir, real mHmult)
592 print_top_comment(out, filename, title, ffdir, bITP);
594 print_top_heavy_H(out, mHmult);
595 fprintf(out, "; Include forcefield parameters\n");
597 p = strrchr(ffdir, '/');
598 p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
600 fprintf(out, "#include \"%s/%s\"\n\n", p, fflib_forcefield_itp());
603 static void print_top_posre(FILE *out, const char *pr)
605 fprintf(out, "; Include Position restraint file\n");
606 fprintf(out, "#ifdef POSRES\n");
607 fprintf(out, "#include \"%s\"\n", pr);
608 fprintf(out, "#endif\n\n");
611 static void print_top_water(FILE *out, const char *ffdir, const char *water)
616 fprintf(out, "; Include water topology\n");
618 p = strrchr(ffdir, '/');
619 p = (ffdir[0] == '.' || p == NULL) ? ffdir : p+1;
620 fprintf(out, "#include \"%s/%s.itp\"\n", p, water);
623 fprintf(out, "#ifdef POSRES_WATER\n");
624 fprintf(out, "; Position restraint for each water oxygen\n");
625 fprintf(out, "[ position_restraints ]\n");
626 fprintf(out, ";%3s %5s %9s %10s %10s\n", "i", "funct", "fcx", "fcy", "fcz");
627 fprintf(out, "%4d %4d %10g %10g %10g\n", 1, 1, 1000.0, 1000.0, 1000.0);
628 fprintf(out, "#endif\n");
631 sprintf(buf, "%s/ions.itp", p);
633 if (fflib_fexist(buf))
635 fprintf(out, "; Include topology for ions\n");
636 fprintf(out, "#include \"%s\"\n", buf);
641 static void print_top_system(FILE *out, const char *title)
643 fprintf(out, "[ %s ]\n", dir2str(d_system));
644 fprintf(out, "; Name\n");
645 fprintf(out, "%s\n\n", title[0] ? title : "Protein");
648 void print_top_mols(FILE *out,
649 const char *title, const char *ffdir, const char *water,
650 int nincl, char **incls, int nmol, t_mols *mols)
657 fprintf(out, "; Include chain topologies\n");
658 for (i = 0; (i < nincl); i++)
660 incl = strrchr(incls[i], DIR_SEPARATOR);
667 /* Remove the path from the include name */
670 fprintf(out, "#include \"%s\"\n", incl);
677 print_top_water(out, ffdir, water);
679 print_top_system(out, title);
683 fprintf(out, "[ %s ]\n", dir2str(d_molecules));
684 fprintf(out, "; %-15s %5s\n", "Compound", "#mols");
685 for (i = 0; (i < nmol); i++)
687 fprintf(out, "%-15s %5d\n", mols[i].name, mols[i].nr);
692 void write_top(FILE *out, char *pr, char *molname,
693 t_atoms *at, gmx_bool bRTPresname,
694 int bts[], t_params plist[], t_excls excls[],
695 gpp_atomtype_t atype, int *cgnr, int nrexcl)
696 /* NOTE: nrexcl is not the size of *excl! */
698 if (at && atype && cgnr)
700 fprintf(out, "[ %s ]\n", dir2str(d_moleculetype));
701 fprintf(out, "; %-15s %5s\n", "Name", "nrexcl");
702 fprintf(out, "%-15s %5d\n\n", molname ? molname : "Protein", nrexcl);
704 print_atoms(out, atype, at, cgnr, bRTPresname);
705 print_bondeds(out, at->nr, d_bonds, F_BONDS, bts[ebtsBONDS], plist);
706 print_bondeds(out, at->nr, d_constraints, F_CONSTR, 0, plist);
707 print_bondeds(out, at->nr, d_constraints, F_CONSTRNC, 0, plist);
708 print_bondeds(out, at->nr, d_pairs, F_LJ14, 0, plist);
709 print_excl(out, at->nr, excls);
710 print_bondeds(out, at->nr, d_angles, F_ANGLES, bts[ebtsANGLES], plist);
711 print_bondeds(out, at->nr, d_dihedrals, F_PDIHS, bts[ebtsPDIHS], plist);
712 print_bondeds(out, at->nr, d_dihedrals, F_IDIHS, bts[ebtsIDIHS], plist);
713 print_bondeds(out, at->nr, d_cmap, F_CMAP, bts[ebtsCMAP], plist);
714 print_bondeds(out, at->nr, d_polarization, F_POLARIZATION, 0, plist);
715 print_bondeds(out, at->nr, d_thole_polarization, F_THOLE_POL, 0, plist);
716 print_bondeds(out, at->nr, d_vsites2, F_VSITE2, 0, plist);
717 print_bondeds(out, at->nr, d_vsites3, F_VSITE3, 0, plist);
718 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FD, 0, plist);
719 print_bondeds(out, at->nr, d_vsites3, F_VSITE3FAD, 0, plist);
720 print_bondeds(out, at->nr, d_vsites3, F_VSITE3OUT, 0, plist);
721 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FD, 0, plist);
722 print_bondeds(out, at->nr, d_vsites4, F_VSITE4FDN, 0, plist);
726 print_top_posre(out, pr);
731 static atom_id search_res_atom(const char *type, int resind,
733 const char *bondtype, gmx_bool bAllowMissing)
737 for (i = 0; (i < atoms->nr); i++)
739 if (atoms->atom[i].resind == resind)
741 return search_atom(type, i, atoms, bondtype, bAllowMissing);
748 static void do_ssbonds(t_params *ps, t_atoms *atoms,
749 int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
754 for (i = 0; (i < nssbonds); i++)
756 ri = ssbonds[i].res1;
757 rj = ssbonds[i].res2;
758 ai = search_res_atom(ssbonds[i].a1, ri, atoms,
759 "special bond", bAllowMissing);
760 aj = search_res_atom(ssbonds[i].a2, rj, atoms,
761 "special bond", bAllowMissing);
762 if ((ai == NO_ATID) || (aj == NO_ATID))
764 gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
765 ssbonds[i].a1, ssbonds[i].a2);
767 add_param(ps, ai, aj, NULL, NULL);
771 static void at2bonds(t_params *psb, t_hackblock *hb,
774 real long_bond_dist, real short_bond_dist)
778 real dist2, long_bond_dist2, short_bond_dist2;
781 long_bond_dist2 = sqr(long_bond_dist);
782 short_bond_dist2 = sqr(short_bond_dist);
793 fprintf(stderr, "Making bonds...\n");
795 for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
797 /* add bonds from list of bonded interactions */
798 for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
800 /* Unfortunately we can not issue errors or warnings
801 * for missing atoms in bonds, as the hydrogens and terminal atoms
802 * have not been added yet.
804 ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[0], i, atoms,
806 aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[1], i, atoms,
808 if (ai != NO_ATID && aj != NO_ATID)
810 dist2 = distance2(x[ai], x[aj]);
811 if (dist2 > long_bond_dist2)
813 fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
814 ai+1, aj+1, sqrt(dist2));
816 else if (dist2 < short_bond_dist2)
818 fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
819 ai+1, aj+1, sqrt(dist2));
821 add_param(psb, ai, aj, NULL, hb[resind].rb[ebtsBONDS].b[j].s);
824 /* add bonds from list of hacks (each added atom gets a bond) */
825 while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
827 for (j = 0; j < hb[resind].nhack; j++)
829 if ( ( hb[resind].hack[j].tp > 0 ||
830 hb[resind].hack[j].oname == NULL ) &&
831 strcmp(hb[resind].hack[j].a[0], *(atoms->atomname[i])) == 0)
833 switch (hb[resind].hack[j].tp)
835 case 9: /* COOH terminus */
836 add_param(psb, i, i+1, NULL, NULL); /* C-O */
837 add_param(psb, i, i+2, NULL, NULL); /* C-OA */
838 add_param(psb, i+2, i+3, NULL, NULL); /* OA-H */
841 for (k = 0; (k < hb[resind].hack[j].nr); k++)
843 add_param(psb, i, i+k+1, NULL, NULL);
850 /* we're now at the start of the next residue */
854 static int pcompar(const void *a, const void *b)
861 d = pa->a[0] - pb->a[0];
864 d = pa->a[1] - pb->a[1];
868 return strlen(pb->s) - strlen(pa->s);
876 static void clean_bonds(t_params *ps)
883 /* swap atomnumbers in bond if first larger than second: */
884 for (i = 0; (i < ps->nr); i++)
886 if (ps->param[i].a[1] < ps->param[i].a[0])
888 a = ps->param[i].a[0];
889 ps->param[i].a[0] = ps->param[i].a[1];
890 ps->param[i].a[1] = a;
895 qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
897 /* remove doubles, keep the first one always. */
899 for (i = 1; (i < ps->nr); i++)
901 if ((ps->param[i].a[0] != ps->param[j-1].a[0]) ||
902 (ps->param[i].a[1] != ps->param[j-1].a[1]) )
906 cp_param(&(ps->param[j]), &(ps->param[i]));
911 fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
916 fprintf(stderr, "No bonds\n");
920 void print_sums(t_atoms *atoms, gmx_bool bSystem)
928 where = " in system";
937 for (i = 0; (i < atoms->nr); i++)
939 m += atoms->atom[i].m;
940 qtot += atoms->atom[i].q;
943 fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
944 fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
947 static void check_restp_type(const char *name, int t1, int t2)
951 gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
955 static void check_restp_types(t_restp *r0, t_restp *r1)
959 check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
960 check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
961 check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
962 check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
964 for (i = 0; i < ebtsNR; i++)
966 check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
970 void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
974 const char *Hnum = "123456";
978 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
980 * restp->atomname[at_start], resnr, restp->resname);
982 strcpy(buf, hack->nname);
983 buf[strlen(buf)+1] = '\0';
986 buf[strlen(buf)] = '-';
989 restp->natom += hack->nr;
990 srenew(restp->atom, restp->natom);
991 srenew(restp->atomname, restp->natom);
992 srenew(restp->cgnr, restp->natom);
994 for (k = restp->natom-1; k > at_start+hack->nr; k--)
997 restp->atom [k - hack->nr];
999 restp->atomname[k - hack->nr];
1001 restp->cgnr [k - hack->nr];
1004 for (k = 0; k < hack->nr; k++)
1006 /* set counter in atomname */
1009 buf[strlen(buf)-1] = Hnum[k];
1011 snew( restp->atomname[at_start+1+k], 1);
1012 restp->atom [at_start+1+k] = *hack->atom;
1013 *restp->atomname[at_start+1+k] = strdup(buf);
1014 if (hack->cgnr != NOTSET)
1016 restp->cgnr [at_start+1+k] = hack->cgnr;
1020 restp->cgnr [at_start+1+k] = restp->cgnr[at_start];
1025 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1026 int nrtp, t_restp rtp[],
1027 int nres, t_resinfo *resinfo,
1029 t_hackblock **ntdb, t_hackblock **ctdb,
1040 /* first the termini */
1041 for (i = 0; i < nterpairs; i++)
1043 if (rn[i] >= 0 && ntdb[i] != NULL)
1045 copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1047 if (rc[i] >= 0 && ctdb[i] != NULL)
1049 merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1053 /* then the whole rtp */
1054 for (i = 0; i < nres; i++)
1056 /* Here we allow a mismatch of one character when looking for the rtp entry.
1057 * For such a mismatch there should be only one mismatching name.
1058 * This is mainly useful for small molecules such as ions.
1059 * Note that this will usually not work for protein, DNA and RNA,
1060 * since there the residue names should be listed in residuetypes.dat
1061 * and an error will have been generated earlier in the process.
1063 key = *resinfo[i].rtp;
1064 snew(resinfo[i].rtp, 1);
1065 *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1066 res = get_restp(*resinfo[i].rtp, nrtp, rtp);
1067 copy_t_restp(res, &(*restp)[i]);
1069 /* Check that we do not have different bonded types in one molecule */
1070 check_restp_types(&(*restp)[0], &(*restp)[i]);
1073 for (j = 0; j < nterpairs && tern == -1; j++)
1081 for (j = 0; j < nterpairs && terc == -1; j++)
1088 bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1090 if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
1091 (terc >= 0 && ctdb[terc] == NULL)))
1093 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).");
1095 if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1096 (terc >= 0 && ctdb[terc]->nhack == 0)))
1098 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.");
1102 /* now perform t_hack's on t_restp's,
1103 i.e. add's and deletes from termini database will be
1104 added to/removed from residue topology
1105 we'll do this on one big dirty loop, so it won't make easy reading! */
1106 for (i = 0; i < nres; i++)
1108 for (j = 0; j < (*hb)[i].nhack; j++)
1110 if ( (*hb)[i].hack[j].nr)
1112 /* find atom in restp */
1113 for (l = 0; l < (*restp)[i].natom; l++)
1115 if ( ( (*hb)[i].hack[j].oname == NULL &&
1116 strcmp((*hb)[i].hack[j].a[0], *(*restp)[i].atomname[l]) == 0 ) ||
1117 ( (*hb)[i].hack[j].oname != NULL &&
1118 strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1123 if (l == (*restp)[i].natom)
1125 /* If we are doing an atom rename only, we don't need
1126 * to generate a fatal error if the old name is not found
1129 /* Deleting can happen also only on the input atoms,
1130 * not necessarily always on the rtp entry.
1132 if (!((*hb)[i].hack[j].oname != NULL &&
1133 (*hb)[i].hack[j].nname != NULL) &&
1134 !((*hb)[i].hack[j].oname != NULL &&
1135 (*hb)[i].hack[j].nname == NULL))
1138 "atom %s not found in buiding block %d%s "
1139 "while combining tdb and rtp",
1140 (*hb)[i].hack[j].oname != NULL ?
1141 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].a[0],
1142 i+1, *resinfo[i].rtp);
1147 if ( (*hb)[i].hack[j].oname == NULL)
1150 add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1155 if ( (*hb)[i].hack[j].nname == NULL)
1157 /* we're deleting */
1160 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1161 *(*restp)[i].atomname[l],
1162 i+1, (*restp)[i].resname);
1164 /* shift the rest */
1165 (*restp)[i].natom--;
1166 for (k = l; k < (*restp)[i].natom; k++)
1168 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
1169 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1170 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
1172 /* give back space */
1173 srenew((*restp)[i].atom, (*restp)[i].natom);
1174 srenew((*restp)[i].atomname, (*restp)[i].natom);
1175 srenew((*restp)[i].cgnr, (*restp)[i].natom);
1177 else /* nname != NULL */
1178 { /* we're replacing */
1181 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1182 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1183 i+1, (*restp)[i].resname);
1185 snew( (*restp)[i].atomname[l], 1);
1186 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
1187 *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
1188 if ( (*hb)[i].hack[j].cgnr != NOTSET)
1190 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
1200 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1207 return (gmx_strcasecmp(anm, hack->nname) == 0);
1211 if (isdigit(anm[strlen(anm)-1]))
1213 *nr = anm[strlen(anm)-1] - '0';
1219 if (*nr <= 0 || *nr > hack->nr)
1225 return (strlen(anm) == strlen(hack->nname) + 1 &&
1226 gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1231 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1232 t_restp *rptr, t_hackblock *hbr,
1237 char *oldnm, *newnm;
1239 char *start_at, buf[STRLEN];
1241 gmx_bool bReplaceReplace, bFoundInAdd;
1244 oldnm = *pdba->atomname[atind];
1245 resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1248 for (j = 0; j < hbr->nhack; j++)
1250 if (hbr->hack[j].oname != NULL && hbr->hack[j].nname != NULL &&
1251 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1253 /* This is a replace entry. */
1254 /* Check if we are not replacing a replaced atom. */
1255 bReplaceReplace = FALSE;
1256 for (k = 0; k < hbr->nhack; k++)
1259 hbr->hack[k].oname != NULL && hbr->hack[k].nname != NULL &&
1260 gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1262 /* The replace in hack[j] replaces an atom that
1263 * was already replaced in hack[k], we do not want
1264 * second or higher level replaces at this stage.
1266 bReplaceReplace = TRUE;
1269 if (bReplaceReplace)
1271 /* Skip this replace. */
1275 /* This atom still has the old name, rename it */
1276 newnm = hbr->hack[j].nname;
1277 for (k = 0; k < rptr->natom; k++)
1279 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1284 if (k == rptr->natom)
1286 /* The new name is not present in the rtp.
1287 * We need to apply the replace also to the rtp entry.
1290 /* We need to find the add hack that can add this atom
1291 * to find out after which atom it should be added.
1293 bFoundInAdd = FALSE;
1294 for (k = 0; k < hbr->nhack; k++)
1296 if (hbr->hack[k].oname == NULL &&
1297 hbr->hack[k].nname != NULL &&
1298 atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1302 start_at = hbr->hack[k].a[0];
1306 sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1309 for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1311 if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1316 if (start_nr == rptr->natom)
1318 gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1319 start_at, rptr->resname, newnm);
1321 /* We can add the atom after atom start_nr */
1322 add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1330 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'",
1332 hbr->hack[j].oname, hbr->hack[j].nname,
1339 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1340 oldnm, rptr->resname, resnr, newnm);
1342 /* Rename the atom in pdba */
1343 snew(pdba->atomname[atind], 1);
1344 *pdba->atomname[atind] = strdup(newnm);
1346 else if (hbr->hack[j].oname != NULL && hbr->hack[j].nname == NULL &&
1347 gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1349 /* This is a delete entry, check if this atom is present
1350 * in the rtp entry of this residue.
1352 for (k = 0; k < rptr->natom; k++)
1354 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1359 if (k == rptr->natom)
1361 /* This atom is not present in the rtp entry,
1362 * delete is now from the input pdba.
1366 printf("Deleting atom '%s' in residue '%s' %d\n",
1367 oldnm, rptr->resname, resnr);
1369 /* We should free the atom name,
1370 * but it might be used multiple times in the symtab.
1371 * sfree(pdba->atomname[atind]);
1373 for (k = atind+1; k < pdba->nr; k++)
1375 pdba->atom[k-1] = pdba->atom[k];
1376 pdba->atomname[k-1] = pdba->atomname[k];
1377 copy_rvec(x[k], x[k-1]);
1388 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1389 t_atoms *pdba, rvec *x,
1396 for (i = 0; i < pdba->nr; i++)
1398 oldnm = *pdba->atomname[i];
1399 rptr = &restp[pdba->atom[i].resind];
1400 for (j = 0; (j < rptr->natom); j++)
1402 if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1407 if (j == rptr->natom)
1409 /* Not found yet, check if we have to rename this atom */
1410 if (match_atomnames_with_rtp_atom(pdba, x, i,
1411 rptr, &(hb[pdba->atom[i].resind]),
1414 /* We deleted this atom, decrease the atom counter by 1. */
1421 #define NUM_CMAP_ATOMS 5
1422 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms, gmx_residuetype_t rt)
1424 int residx, i, j, k;
1426 t_resinfo *resinfo = atoms->resinfo;
1427 int nres = atoms->nres;
1429 atom_id cmap_atomid[NUM_CMAP_ATOMS];
1430 int cmap_chainnum = -1, this_residue_index;
1441 fprintf(stderr, "Making cmap torsions...");
1443 /* End loop at nres-1, since the very last residue does not have a +N atom, and
1444 * therefore we get a valgrind invalid 4 byte read error with atom am */
1445 for (residx = 0; residx < nres-1; residx++)
1447 /* Add CMAP terms from the list of CMAP interactions */
1448 for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1451 /* Loop over atoms in a candidate CMAP interaction and
1452 * check that they exist, are from the same chain and are
1453 * from residues labelled as protein. */
1454 for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1456 cmap_atomid[k] = search_atom(restp[residx].rb[ebtsCMAP].b[j].a[k],
1457 i, atoms, ptr, TRUE);
1458 bAddCMAP = bAddCMAP && (cmap_atomid[k] != NO_ATID);
1461 /* This break is necessary, because cmap_atomid[k]
1462 * == NO_ATID cannot be safely used as an index
1463 * into the atom array. */
1466 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1469 cmap_chainnum = resinfo[this_residue_index].chainnum;
1473 /* Does the residue for this atom have the same
1474 * chain number as the residues for previous
1476 bAddCMAP = bAddCMAP &&
1477 cmap_chainnum == resinfo[this_residue_index].chainnum;
1479 bAddCMAP = bAddCMAP && gmx_residuetype_is_protein(rt, *(resinfo[this_residue_index].name));
1484 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);
1488 if (residx < nres-1)
1490 while (atoms->atom[i].resind < residx+1)
1497 /* Start the next residue */
1501 scrub_charge_groups(int *cgnr, int natoms)
1505 for (i = 0; i < natoms; i++)
1512 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1513 t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1514 int nrtp, t_restp rtp[],
1515 t_restp *restp, t_hackblock *hb,
1516 gmx_bool bAllowMissing,
1517 gmx_bool bVsites, gmx_bool bVsiteAromatics,
1520 int nssbonds, t_ssbond *ssbonds,
1521 real long_bond_dist, real short_bond_dist,
1522 gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1523 gmx_bool bRenumRes, gmx_bool bRTPresname)
1529 t_params plist[F_NRE];
1536 gmx_residuetype_t rt;
1539 gmx_residuetype_init(&rt);
1543 print_resall(debug, atoms->nres, restp, atype);
1544 dump_hb(debug, atoms->nres, hb);
1548 at2bonds(&(plist[F_BONDS]), hb,
1550 long_bond_dist, short_bond_dist);
1552 /* specbonds: disulphide bonds & heme-his */
1553 do_ssbonds(&(plist[F_BONDS]),
1554 atoms, nssbonds, ssbonds,
1557 nmissat = name2type(atoms, &cgnr, atype, restp, rt);
1562 fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1567 gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1572 /* Cleanup bonds (sort and rm doubles) */
1573 clean_bonds(&(plist[F_BONDS]));
1575 snew(vsite_type, atoms->nr);
1576 for (i = 0; i < atoms->nr; i++)
1578 vsite_type[i] = NOTSET;
1582 /* determine which atoms will be vsites and add dummy masses
1583 also renumber atom numbers in plist[0..F_NRE]! */
1584 do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1585 &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1588 /* Make Angles and Dihedrals */
1589 fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1590 snew(excls, atoms->nr);
1591 init_nnb(&nnb, atoms->nr, 4);
1592 gen_nnb(&nnb, plist);
1593 print_nnb(&nnb, "NNB");
1594 gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1600 gen_cmap(&(plist[F_CMAP]), restp, atoms, rt);
1601 if (plist[F_CMAP].nr > 0)
1603 fprintf(stderr, "There are %4d cmap torsion pairs\n",
1608 /* set mass of all remaining hydrogen atoms */
1611 do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1615 /* Cleanup bonds (sort and rm doubles) */
1616 /* clean_bonds(&(plist[F_BONDS]));*/
1619 "There are %4d dihedrals, %4d impropers, %4d angles\n"
1620 " %4d pairs, %4d bonds and %4d virtual sites\n",
1621 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1622 plist[F_LJ14].nr, plist[F_BONDS].nr,
1623 plist[F_VSITE2].nr +
1624 plist[F_VSITE3].nr +
1625 plist[F_VSITE3FD].nr +
1626 plist[F_VSITE3FAD].nr +
1627 plist[F_VSITE3OUT].nr +
1628 plist[F_VSITE4FD].nr +
1629 plist[F_VSITE4FDN].nr );
1631 print_sums(atoms, FALSE);
1633 if (FALSE == bChargeGroups)
1635 scrub_charge_groups(cgnr, atoms->nr);
1640 for (i = 0; i < atoms->nres; i++)
1642 atoms->resinfo[i].nr = i + 1;
1643 atoms->resinfo[i].ic = ' ';
1649 fprintf(stderr, "Writing topology\n");
1650 /* We can copy the bonded types from the first restp,
1651 * since the types have to be identical for all residues in one molecule.
1653 for (i = 0; i < ebtsNR; i++)
1655 bts[i] = restp[0].rb[i].type;
1657 write_top(top_file, posre_fn, molname,
1659 bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1663 free_t_hackblock(atoms->nres, &hb);
1664 free_t_restp(atoms->nres, &restp);
1665 gmx_residuetype_destroy(rt);
1667 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1669 for (i = 0; i < F_NRE; i++)
1671 sfree(plist[i].param);
1673 for (i = 0; i < atoms->nr; i++)