clang-tidy: readability-non-const-parameter (2/2)
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / pdb2top.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5  * Copyright (c) 2001-2004, The GROMACS development team.
6  * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8  * and including many others, as listed in the AUTHORS file in the
9  * top-level source directory and at http://www.gromacs.org.
10  *
11  * GROMACS is free software; you can redistribute it and/or
12  * modify it under the terms of the GNU Lesser General Public License
13  * as published by the Free Software Foundation; either version 2.1
14  * of the License, or (at your option) any later version.
15  *
16  * GROMACS is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19  * Lesser General Public License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with GROMACS; if not, see
23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
25  *
26  * If you want to redistribute modifications to GROMACS, please
27  * consider that scientific software is very special. Version
28  * control is crucial - bugs must be traceable. We will be happy to
29  * consider code for inclusion in the official distribution, but
30  * derived work must not be called official GROMACS. Details are found
31  * in the README & COPYING files - if they are missing, get the
32  * official version at http://www.gromacs.org.
33  *
34  * To help us fund GROMACS development, we humbly ask that you cite
35  * the research papers on the package. Check out http://www.gromacs.org.
36  */
37 #include "gmxpre.h"
38
39 #include "pdb2top.h"
40
41 #include <ctype.h>
42 #include <stdio.h>
43 #include <string.h>
44
45 #include <cmath>
46
47 #include <algorithm>
48 #include <string>
49 #include <vector>
50
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/gmxpreprocess/add_par.h"
53 #include "gromacs/gmxpreprocess/fflibutil.h"
54 #include "gromacs/gmxpreprocess/gen_ad.h"
55 #include "gromacs/gmxpreprocess/gen_vsite.h"
56 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
57 #include "gromacs/gmxpreprocess/h_db.h"
58 #include "gromacs/gmxpreprocess/notset.h"
59 #include "gromacs/gmxpreprocess/pgutil.h"
60 #include "gromacs/gmxpreprocess/resall.h"
61 #include "gromacs/gmxpreprocess/topdirs.h"
62 #include "gromacs/gmxpreprocess/topio.h"
63 #include "gromacs/gmxpreprocess/toputil.h"
64 #include "gromacs/math/functions.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/topology/residuetypes.h"
67 #include "gromacs/topology/symtab.h"
68 #include "gromacs/utility/binaryinformation.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/dir_separator.h"
71 #include "gromacs/utility/exceptions.h"
72 #include "gromacs/utility/fatalerror.h"
73 #include "gromacs/utility/futil.h"
74 #include "gromacs/utility/niceheader.h"
75 #include "gromacs/utility/path.h"
76 #include "gromacs/utility/programcontext.h"
77 #include "gromacs/utility/smalloc.h"
78 #include "gromacs/utility/strdb.h"
79 #include "gromacs/utility/stringutil.h"
80 #include "gromacs/utility/textwriter.h"
81
82 /* this must correspond to enum in pdb2top.h */
83 const char *hh[ehisNR]   = { "HISD", "HISE", "HISH", "HIS1" };
84
85 static int missing_atoms(t_restp *rp, int resind, t_atoms *at, int i0, int i)
86 {
87     int      j, k, nmiss;
88     char    *name;
89     gmx_bool bFound;
90
91     nmiss = 0;
92     for (j = 0; j < rp->natom; j++)
93     {
94         name   = *(rp->atomname[j]);
95         bFound = FALSE;
96         for (k = i0; k < i; k++)
97         {
98             bFound = (bFound || !gmx_strcasecmp(*(at->atomname[k]), name));
99         }
100         if (!bFound)
101         {
102             nmiss++;
103             fprintf(stderr, "\nWARNING: "
104                     "atom %s is missing in residue %s %d in the pdb file\n",
105                     name, *(at->resinfo[resind].name), at->resinfo[resind].nr);
106             if (name[0] == 'H' || name[0] == 'h')
107             {
108                 fprintf(stderr, "         You might need to add atom %s to the hydrogen database of building block %s\n"
109                         "         in the file %s.hdb (see the manual)\n",
110                         name, *(at->resinfo[resind].rtp), rp->filebase);
111             }
112             fprintf(stderr, "\n");
113         }
114     }
115
116     return nmiss;
117 }
118
119 gmx_bool is_int(double x)
120 {
121     const double tol = 1e-4;
122     int          ix;
123
124     if (x < 0)
125     {
126         x = -x;
127     }
128     ix = std::round(x);
129
130     return (fabs(x-ix) < tol);
131 }
132
133 static void
134 choose_ff_impl(const char *ffsel,
135                char *forcefield, int ff_maxlen,
136                char *ffdir, int ffdir_maxlen)
137 {
138     std::vector<gmx::DataFileInfo> ffdirs = fflib_enumerate_forcefields();
139     const int nff = static_cast<int>(ffdirs.size());
140
141     /* Replace with unix path separators */
142     if (DIR_SEPARATOR != '/')
143     {
144         for (int i = 0; i < nff; ++i)
145         {
146             std::replace(ffdirs[i].dir.begin(), ffdirs[i].dir.end(), DIR_SEPARATOR, '/');
147         }
148     }
149
150     /* Store the force field names in ffs */
151     std::vector<std::string> ffs;
152     ffs.reserve(ffdirs.size());
153     for (int i = 0; i < nff; ++i)
154     {
155         ffs.push_back(gmx::stripSuffixIfPresent(ffdirs[i].name,
156                                                 fflib_forcefield_dir_ext()));
157     }
158
159     int sel;
160     if (ffsel != nullptr)
161     {
162         sel         = -1;
163         int cwdsel  = -1;
164         int nfound  = 0;
165         for (int i = 0; i < nff; ++i)
166         {
167             if (ffs[i] == ffsel)
168             {
169                 /* Matching ff name */
170                 sel = i;
171                 nfound++;
172
173                 if (ffdirs[i].dir == ".")
174                 {
175                     cwdsel = i;
176                 }
177             }
178         }
179
180         if (cwdsel != -1)
181         {
182             sel = cwdsel;
183         }
184
185         if (nfound > 1)
186         {
187             if (cwdsel != -1)
188             {
189                 fprintf(stderr,
190                         "Force field '%s' occurs in %d places. pdb2gmx is using the one in the\n"
191                         "current directory. Use interactive selection (not the -ff option) if\n"
192                         "you would prefer a different one.\n", ffsel, nfound);
193             }
194             else
195             {
196                 std::string message = gmx::formatString(
197                             "Force field '%s' occurs in %d places, but not in "
198                             "the current directory.\n"
199                             "Run without the -ff switch and select the force "
200                             "field interactively.", ffsel, nfound);
201                 GMX_THROW(gmx::InconsistentInputError(message));
202             }
203         }
204         else if (nfound == 0)
205         {
206             std::string message = gmx::formatString(
207                         "Could not find force field '%s' in current directory, "
208                         "install tree or GMXLIB path.", ffsel);
209             GMX_THROW(gmx::InconsistentInputError(message));
210         }
211     }
212     else if (nff > 1)
213     {
214         std::vector<std::string> desc;
215         desc.reserve(ffdirs.size());
216         for (int i = 0; i < nff; ++i)
217         {
218             std::string docFileName(
219                     gmx::Path::join(ffdirs[i].dir, ffdirs[i].name,
220                                     fflib_forcefield_doc()));
221             // TODO: Just try to open the file with a method that does not
222             // throw/bail out with a fatal error instead of multiple checks.
223             if (gmx::File::exists(docFileName, gmx::File::returnFalseOnError))
224             {
225                 // TODO: Use a C++ API without such an intermediate/fixed-length buffer.
226                 char  buf[STRLEN];
227                 /* We don't use fflib_open, because we don't want printf's */
228                 FILE *fp = gmx_ffopen(docFileName.c_str(), "r");
229                 get_a_line(fp, buf, STRLEN);
230                 gmx_ffclose(fp);
231                 desc.emplace_back(buf);
232             }
233             else
234             {
235                 desc.push_back(ffs[i]);
236             }
237         }
238         /* Order force fields from the same dir alphabetically
239          * and put deprecated force fields at the end.
240          */
241         for (int i = 0; i < nff; ++i)
242         {
243             for (int j = i + 1; j < nff; ++j)
244             {
245                 if (ffdirs[i].dir == ffdirs[j].dir &&
246                     ((desc[i][0] == '[' && desc[j][0] != '[') ||
247                      ((desc[i][0] == '[' || desc[j][0] != '[') &&
248                       gmx_strcasecmp(desc[i].c_str(), desc[j].c_str()) > 0)))
249                 {
250                     std::swap(ffdirs[i].name, ffdirs[j].name);
251                     std::swap(ffs[i], ffs[j]);
252                     std::swap(desc[i], desc[j]);
253                 }
254             }
255         }
256
257         printf("\nSelect the Force Field:\n");
258         for (int i = 0; i < nff; ++i)
259         {
260             if (i == 0 || ffdirs[i-1].dir != ffdirs[i].dir)
261             {
262                 if (ffdirs[i].dir == ".")
263                 {
264                     printf("From current directory:\n");
265                 }
266                 else
267                 {
268                     printf("From '%s':\n", ffdirs[i].dir.c_str());
269                 }
270             }
271             printf("%2d: %s\n", i+1, desc[i].c_str());
272         }
273
274         sel = -1;
275         // TODO: Add a C++ API for this.
276         char   buf[STRLEN];
277         char  *pret;
278         do
279         {
280             pret = fgets(buf, STRLEN, stdin);
281
282             if (pret != nullptr)
283             {
284                 sel = strtol(buf, nullptr, 10);
285                 sel--;
286             }
287         }
288         while (pret == nullptr || (sel < 0) || (sel >= nff));
289
290         /* Check for a current limitation of the fflib code.
291          * It will always read from the first ff directory in the list.
292          * This check assumes that the order of ffs matches the order
293          * in which fflib_open searches ff library files.
294          */
295         for (int i = 0; i < sel; i++)
296         {
297             if (ffs[i] == ffs[sel])
298             {
299                 std::string message = gmx::formatString(
300                             "Can only select the first of multiple force "
301                             "field entries with directory name '%s%s' in "
302                             "the list. If you want to use the next entry, "
303                             "run pdb2gmx in a different directory, set GMXLIB "
304                             "to point to the desired force field first, and/or "
305                             "rename or move the force field directory present "
306                             "in the current working directory.",
307                             ffs[sel].c_str(), fflib_forcefield_dir_ext());
308                 GMX_THROW(gmx::NotImplementedError(message));
309             }
310         }
311     }
312     else
313     {
314         sel = 0;
315     }
316
317     if (ffs[sel].length() >= static_cast<size_t>(ff_maxlen))
318     {
319         std::string message = gmx::formatString(
320                     "Length of force field name (%d) >= maxlen (%d)",
321                     static_cast<int>(ffs[sel].length()), ff_maxlen);
322         GMX_THROW(gmx::InvalidInputError(message));
323     }
324     strcpy(forcefield, ffs[sel].c_str());
325
326     std::string ffpath;
327     if (ffdirs[sel].bFromDefaultDir)
328     {
329         ffpath = ffdirs[sel].name;
330     }
331     else
332     {
333         ffpath = gmx::Path::join(ffdirs[sel].dir, ffdirs[sel].name);
334     }
335     if (ffpath.length() >= static_cast<size_t>(ffdir_maxlen))
336     {
337         std::string message = gmx::formatString(
338                     "Length of force field dir (%d) >= maxlen (%d)",
339                     static_cast<int>(ffpath.length()), ffdir_maxlen);
340         GMX_THROW(gmx::InvalidInputError(message));
341     }
342     strcpy(ffdir, ffpath.c_str());
343 }
344
345 void
346 choose_ff(const char *ffsel,
347           char *forcefield, int ff_maxlen,
348           char *ffdir, int ffdir_maxlen)
349 {
350     try
351     {
352         choose_ff_impl(ffsel, forcefield, ff_maxlen, ffdir, ffdir_maxlen);
353     }
354     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
355 }
356
357 void choose_watermodel(const char *wmsel, const char *ffdir,
358                        char **watermodel)
359 {
360     const char *fn_watermodels = "watermodels.dat";
361     char        fn_list[STRLEN];
362     FILE       *fp;
363     char        buf[STRLEN];
364     int         nwm, sel, i;
365     char      **model;
366     char       *pret;
367
368     if (strcmp(wmsel, "none") == 0)
369     {
370         *watermodel = nullptr;
371
372         return;
373     }
374     else if (strcmp(wmsel, "select") != 0)
375     {
376         *watermodel = gmx_strdup(wmsel);
377
378         return;
379     }
380
381     sprintf(fn_list, "%s%c%s", ffdir, DIR_SEPARATOR, fn_watermodels);
382
383     if (!fflib_fexist(fn_list))
384     {
385         fprintf(stderr, "No file '%s' found, will not include a water model\n",
386                 fn_watermodels);
387         *watermodel = nullptr;
388
389         return;
390     }
391
392     fp = fflib_open(fn_list);
393     printf("\nSelect the Water Model:\n");
394     nwm   = 0;
395     model = nullptr;
396     while (get_a_line(fp, buf, STRLEN))
397     {
398         srenew(model, nwm+1);
399         snew(model[nwm], STRLEN);
400         sscanf(buf, "%s%n", model[nwm], &i);
401         if (i > 0)
402         {
403             ltrim(buf+i);
404             fprintf(stderr, "%2d: %s\n", nwm+1, buf+i);
405             nwm++;
406         }
407         else
408         {
409             sfree(model[nwm]);
410         }
411     }
412     gmx_ffclose(fp);
413     fprintf(stderr, "%2d: %s\n", nwm+1, "None");
414
415     sel = -1;
416     do
417     {
418         pret = fgets(buf, STRLEN, stdin);
419
420         if (pret != nullptr)
421         {
422             sel = strtol(buf, nullptr, 10);
423             sel--;
424         }
425     }
426     while (pret == nullptr || sel < 0 || sel > nwm);
427
428     if (sel == nwm)
429     {
430         *watermodel = nullptr;
431     }
432     else
433     {
434         *watermodel = gmx_strdup(model[sel]);
435     }
436
437     for (i = 0; i < nwm; i++)
438     {
439         sfree(model[i]);
440     }
441     sfree(model);
442 }
443
444 static int name2type(t_atoms *at, int **cgnr,
445                      t_restp restp[], gmx_residuetype_t *rt)
446 {
447     int         i, j, prevresind, resind, i0, prevcg, cg, curcg;
448     char       *name;
449     gmx_bool    bNterm;
450     double      qt;
451     int         nmissat;
452
453     nmissat = 0;
454
455     resind = -1;
456     bNterm = FALSE;
457     i0     = 0;
458     snew(*cgnr, at->nr);
459     qt     = 0;
460     prevcg = NOTSET;
461     curcg  = 0;
462     cg     = -1;
463     j      = NOTSET;
464
465     for (i = 0; (i < at->nr); i++)
466     {
467         prevresind = resind;
468         if (at->atom[i].resind != resind)
469         {
470             gmx_bool bProt;
471             resind = at->atom[i].resind;
472             bProt  = gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name));
473             bNterm = bProt && (resind == 0);
474             if (resind > 0)
475             {
476                 nmissat += missing_atoms(&restp[prevresind], prevresind, at, i0, i);
477             }
478             i0 = i;
479         }
480         if (at->atom[i].m == 0)
481         {
482             if (debug)
483             {
484                 fprintf(debug, "atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
485                         i+1, *(at->atomname[i]), curcg, prevcg,
486                         j == NOTSET ? NOTSET : restp[resind].cgnr[j]);
487             }
488             qt               = 0;
489             prevcg           = cg;
490             name             = *(at->atomname[i]);
491             j                = search_jtype(&restp[resind], name, bNterm);
492             at->atom[i].type = restp[resind].atom[j].type;
493             at->atom[i].q    = restp[resind].atom[j].q;
494             at->atom[i].m    = restp[resind].atom[j].m;
495             cg               = restp[resind].cgnr[j];
496             /* A charge group number -1 signals a separate charge group
497              * for this atom.
498              */
499             if ( (cg == -1) || (cg != prevcg) || (resind != prevresind) )
500             {
501                 curcg++;
502             }
503         }
504         else
505         {
506             if (debug)
507             {
508                 fprintf(debug, "atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
509                         i+1, *(at->atomname[i]), curcg, qt, is_int(qt));
510             }
511             cg = -1;
512             if (is_int(qt))
513             {
514                 qt = 0;
515                 curcg++;
516             }
517             qt += at->atom[i].q;
518         }
519         (*cgnr)[i]        = curcg;
520         at->atom[i].typeB = at->atom[i].type;
521         at->atom[i].qB    = at->atom[i].q;
522         at->atom[i].mB    = at->atom[i].m;
523     }
524     nmissat += missing_atoms(&restp[resind], resind, at, i0, i);
525
526     return nmissat;
527 }
528
529 static void print_top_heavy_H(FILE *out, real mHmult)
530 {
531     if (mHmult == 2.0)
532     {
533         fprintf(out, "; Using deuterium instead of hydrogen\n\n");
534     }
535     else if (mHmult == 4.0)
536     {
537         fprintf(out, "#define HEAVY_H\n\n");
538     }
539     else if (mHmult != 1.0)
540     {
541         fprintf(stderr, "WARNING: unsupported proton mass multiplier (%g) "
542                 "in pdb2top\n", mHmult);
543     }
544 }
545
546 void print_top_comment(FILE       *out,
547                        const char *filename,
548                        const char *ffdir,
549                        gmx_bool    bITP)
550 {
551     char  ffdir_parent[STRLEN];
552     char *p;
553
554     try
555     {
556         gmx::TextWriter writer(out);
557         gmx::niceHeader(&writer, filename, ';');
558         writer.writeLine(gmx::formatString(";\tThis is a %s topology file", bITP ? "include" : "standalone"));
559         writer.writeLine(";");
560
561         gmx::BinaryInformationSettings settings;
562         settings.generatedByHeader(true);
563         settings.linePrefix(";\t");
564         gmx::printBinaryInformation(&writer, gmx::getProgramContext(), settings);
565     }
566     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
567
568     if (strchr(ffdir, '/') == nullptr)
569     {
570         fprintf(out, ";\tForce field was read from the standard GROMACS share directory.\n;\n\n");
571     }
572     else if (ffdir[0] == '.')
573     {
574         fprintf(out, ";\tForce field was read from current directory or a relative path - path added.\n;\n\n");
575     }
576     else
577     {
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, '/');
581
582         *p = '\0';
583
584         fprintf(out,
585                 ";\tForce field data was read from:\n"
586                 ";\t%s\n"
587                 ";\n"
588                 ";\tNote:\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",
592                 ffdir_parent);
593     }
594 }
595
596 void print_top_header(FILE *out, const char *filename,
597                       gmx_bool bITP, const char *ffdir, real mHmult)
598 {
599     const char *p;
600
601     print_top_comment(out, filename, ffdir, bITP);
602
603     print_top_heavy_H(out, mHmult);
604     fprintf(out, "; Include forcefield parameters\n");
605
606     p = strrchr(ffdir, '/');
607     p = (ffdir[0] == '.' || p == nullptr) ? ffdir : p+1;
608
609     fprintf(out, "#include \"%s/%s\"\n\n", p, fflib_forcefield_itp());
610 }
611
612 static void print_top_posre(FILE *out, const char *pr)
613 {
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");
618 }
619
620 static void print_top_water(FILE *out, const char *ffdir, const char *water)
621 {
622     const char *p;
623     char        buf[STRLEN];
624
625     fprintf(out, "; Include water topology\n");
626
627     p = strrchr(ffdir, '/');
628     p = (ffdir[0] == '.' || p == nullptr) ? ffdir : p+1;
629     fprintf(out, "#include \"%s/%s.itp\"\n", p, water);
630
631     fprintf(out, "\n");
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");
638     fprintf(out, "\n");
639
640     sprintf(buf, "%s/ions.itp", p);
641
642     if (fflib_fexist(buf))
643     {
644         fprintf(out, "; Include topology for ions\n");
645         fprintf(out, "#include \"%s\"\n", buf);
646         fprintf(out, "\n");
647     }
648 }
649
650 static void print_top_system(FILE *out, const char *title)
651 {
652     fprintf(out, "[ %s ]\n", dir2str(d_system));
653     fprintf(out, "; Name\n");
654     fprintf(out, "%s\n\n", title[0] ? title : "Protein");
655 }
656
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)
660 {
661     int   i;
662     char *incl;
663
664     if (nincl > 0)
665     {
666         fprintf(out, "; Include chain topologies\n");
667         for (i = 0; (i < nincl); i++)
668         {
669             incl = strrchr(incls[i], DIR_SEPARATOR);
670             if (incl == nullptr)
671             {
672                 incl = incls[i];
673             }
674             else
675             {
676                 /* Remove the path from the include name */
677                 incl = incl + 1;
678             }
679             fprintf(out, "#include \"%s\"\n", incl);
680         }
681         fprintf(out, "\n");
682     }
683
684     if (water)
685     {
686         print_top_water(out, ffdir, water);
687     }
688     print_top_system(out, title);
689
690     if (nmol)
691     {
692         fprintf(out, "[ %s ]\n", dir2str(d_molecules));
693         fprintf(out, "; %-15s %5s\n", "Compound", "#mols");
694         for (i = 0; (i < nmol); i++)
695         {
696             fprintf(out, "%-15s %5d\n", mols[i].name, mols[i].nr);
697         }
698     }
699 }
700
701 void write_top(FILE *out, char *pr, const 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! */
706 {
707     if (at && atype && cgnr)
708     {
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);
712
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);
732
733         if (pr)
734         {
735             print_top_posre(out, pr);
736         }
737     }
738 }
739
740
741
742 static void do_ssbonds(t_params *ps, t_atoms *atoms,
743                        int nssbonds, t_ssbond *ssbonds, gmx_bool bAllowMissing)
744 {
745     int     i, ri, rj;
746     int     ai, aj;
747
748     for (i = 0; (i < nssbonds); i++)
749     {
750         ri = ssbonds[i].res1;
751         rj = ssbonds[i].res2;
752         ai = search_res_atom(ssbonds[i].a1, ri, atoms,
753                              "special bond", bAllowMissing);
754         aj = search_res_atom(ssbonds[i].a2, rj, atoms,
755                              "special bond", bAllowMissing);
756         if ((ai == -1) || (aj == -1))
757         {
758             gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!",
759                       ssbonds[i].a1, ssbonds[i].a2);
760         }
761         add_param(ps, ai, aj, nullptr, nullptr);
762     }
763 }
764
765 static void at2bonds(t_params *psb, t_hackblock *hb,
766                      t_atoms *atoms,
767                      rvec x[],
768                      real long_bond_dist, real short_bond_dist)
769 {
770     int         resind, i, j, k;
771     int         ai, aj;
772     real        dist2, long_bond_dist2, short_bond_dist2;
773     const char *ptr;
774
775     long_bond_dist2  = gmx::square(long_bond_dist);
776     short_bond_dist2 = gmx::square(short_bond_dist);
777
778     if (debug)
779     {
780         ptr = "bond";
781     }
782     else
783     {
784         ptr = "check";
785     }
786
787     fprintf(stderr, "Making bonds...\n");
788     i = 0;
789     for (resind = 0; (resind < atoms->nres) && (i < atoms->nr); resind++)
790     {
791         /* add bonds from list of bonded interactions */
792         for (j = 0; j < hb[resind].rb[ebtsBONDS].nb; j++)
793         {
794             /* Unfortunately we can not issue errors or warnings
795              * for missing atoms in bonds, as the hydrogens and terminal atoms
796              * have not been added yet.
797              */
798             ai = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[0], i, atoms,
799                              ptr, TRUE);
800             aj = search_atom(hb[resind].rb[ebtsBONDS].b[j].a[1], i, atoms,
801                              ptr, TRUE);
802             if (ai != -1 && aj != -1)
803             {
804                 dist2 = distance2(x[ai], x[aj]);
805                 if (dist2 > long_bond_dist2)
806
807                 {
808                     fprintf(stderr, "Warning: Long Bond (%d-%d = %g nm)\n",
809                             ai+1, aj+1, std::sqrt(dist2));
810                 }
811                 else if (dist2 < short_bond_dist2)
812                 {
813                     fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
814                             ai+1, aj+1, std::sqrt(dist2));
815                 }
816                 add_param(psb, ai, aj, nullptr, hb[resind].rb[ebtsBONDS].b[j].s);
817             }
818         }
819         /* add bonds from list of hacks (each added atom gets a bond) */
820         while ( (i < atoms->nr) && (atoms->atom[i].resind == resind) )
821         {
822             for (j = 0; j < hb[resind].nhack; j++)
823             {
824                 if ( ( hb[resind].hack[j].tp > 0 ||
825                        hb[resind].hack[j].oname == nullptr ) &&
826                      strcmp(hb[resind].hack[j].a[0], *(atoms->atomname[i])) == 0)
827                 {
828                     switch (hb[resind].hack[j].tp)
829                     {
830                         case 9:                                         /* COOH terminus */
831                             add_param(psb, i, i+1, nullptr, nullptr);   /* C-O  */
832                             add_param(psb, i, i+2, nullptr, nullptr);   /* C-OA */
833                             add_param(psb, i+2, i+3, nullptr, nullptr); /* OA-H */
834                             break;
835                         default:
836                             for (k = 0; (k < hb[resind].hack[j].nr); k++)
837                             {
838                                 add_param(psb, i, i+k+1, nullptr, nullptr);
839                             }
840                     }
841                 }
842             }
843             i++;
844         }
845         /* we're now at the start of the next residue */
846     }
847 }
848
849 static int pcompar(const void *a, const void *b)
850 {
851     const t_param *pa, *pb;
852     int            d;
853     pa = static_cast<const t_param *>(a);
854     pb = static_cast<const t_param *>(b);
855
856     d = pa->a[0] - pb->a[0];
857     if (d == 0)
858     {
859         d = pa->a[1] - pb->a[1];
860     }
861     if (d == 0)
862     {
863         return strlen(pb->s) - strlen(pa->s);
864     }
865     else
866     {
867         return d;
868     }
869 }
870
871 static void clean_bonds(t_params *ps)
872 {
873     int     i, j;
874     int     a;
875
876     if (ps->nr > 0)
877     {
878         /* swap atomnumbers in bond if first larger than second: */
879         for (i = 0; (i < ps->nr); i++)
880         {
881             if (ps->param[i].a[1] < ps->param[i].a[0])
882             {
883                 a                 = ps->param[i].a[0];
884                 ps->param[i].a[0] = ps->param[i].a[1];
885                 ps->param[i].a[1] = a;
886             }
887         }
888
889         /* Sort bonds */
890         qsort(ps->param, ps->nr, (size_t)sizeof(ps->param[0]), pcompar);
891
892         /* remove doubles, keep the first one always. */
893         j = 1;
894         for (i = 1; (i < ps->nr); i++)
895         {
896             if ((ps->param[i].a[0] != ps->param[j-1].a[0]) ||
897                 (ps->param[i].a[1] != ps->param[j-1].a[1]) )
898             {
899                 if (j != i)
900                 {
901                     cp_param(&(ps->param[j]), &(ps->param[i]));
902                 }
903                 j++;
904             }
905         }
906         fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j);
907         ps->nr = j;
908     }
909     else
910     {
911         fprintf(stderr, "No bonds\n");
912     }
913 }
914
915 void print_sums(t_atoms *atoms, gmx_bool bSystem)
916 {
917     double      m, qtot;
918     int         i;
919     const char *where;
920
921     if (bSystem)
922     {
923         where = " in system";
924     }
925     else
926     {
927         where = "";
928     }
929
930     m    = 0;
931     qtot = 0;
932     for (i = 0; (i < atoms->nr); i++)
933     {
934         m    += atoms->atom[i].m;
935         qtot += atoms->atom[i].q;
936     }
937
938     fprintf(stderr, "Total mass%s %.3f a.m.u.\n", where, m);
939     fprintf(stderr, "Total charge%s %.3f e\n", where, qtot);
940 }
941
942 static void check_restp_type(const char *name, int t1, int t2)
943 {
944     if (t1 != t2)
945     {
946         gmx_fatal(FARGS, "Residues in one molecule have a different '%s' type: %d and %d", name, t1, t2);
947     }
948 }
949
950 static void check_restp_types(t_restp *r0, t_restp *r1)
951 {
952     int i;
953
954     check_restp_type("all dihedrals", r0->bKeepAllGeneratedDihedrals, r1->bKeepAllGeneratedDihedrals);
955     check_restp_type("nrexcl", r0->nrexcl, r1->nrexcl);
956     check_restp_type("HH14", r0->bGenerateHH14Interactions, r1->bGenerateHH14Interactions);
957     check_restp_type("remove dihedrals", r0->bRemoveDihedralIfWithImproper, r1->bRemoveDihedralIfWithImproper);
958
959     for (i = 0; i < ebtsNR; i++)
960     {
961         check_restp_type(btsNames[i], r0->rb[i].type, r1->rb[i].type);
962     }
963 }
964
965 static void add_atom_to_restp(t_restp *restp, int at_start, const t_hack *hack)
966 {
967     char        buf[STRLEN];
968     int         k;
969     const char *Hnum = "123456";
970
971     /*if (debug)
972        {
973         fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
974                 hack->nname,
975      * restp->atomname[at_start], resnr, restp->resname);
976                 }*/
977     strcpy(buf, hack->nname);
978     buf[strlen(buf)+1] = '\0';
979     if (hack->nr > 1)
980     {
981         buf[strlen(buf)] = '-';
982     }
983     /* make space */
984     restp->natom += hack->nr;
985     srenew(restp->atom,     restp->natom);
986     srenew(restp->atomname, restp->natom);
987     srenew(restp->cgnr,     restp->natom);
988     /* shift the rest */
989     for (k = restp->natom-1; k > at_start+hack->nr; k--)
990     {
991         restp->atom[k] =
992             restp->atom    [k - hack->nr];
993         restp->atomname[k] =
994             restp->atomname[k - hack->nr];
995         restp->cgnr[k] =
996             restp->cgnr    [k - hack->nr];
997     }
998     /* now add them */
999     for (k = 0; k < hack->nr; k++)
1000     {
1001         /* set counter in atomname */
1002         if (hack->nr > 1)
1003         {
1004             buf[strlen(buf)-1] = Hnum[k];
1005         }
1006         snew( restp->atomname[at_start+1+k], 1);
1007         restp->atom     [at_start+1+k] = *hack->atom;
1008         *restp->atomname[at_start+1+k] = gmx_strdup(buf);
1009         if (hack->cgnr != NOTSET)
1010         {
1011             restp->cgnr   [at_start+1+k] = hack->cgnr;
1012         }
1013         else
1014         {
1015             restp->cgnr   [at_start+1+k] = restp->cgnr[at_start];
1016         }
1017     }
1018 }
1019
1020 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
1021                         int nrtp, t_restp rtp[],
1022                         int nres, t_resinfo *resinfo,
1023                         int nterpairs,
1024                         t_hackblock **ntdb, t_hackblock **ctdb,
1025                         const int *rn, const int *rc,
1026                         gmx_bool bAllowMissing)
1027 {
1028     int         i, j, k, l;
1029     char       *key;
1030     t_restp    *res;
1031     int         tern, terc;
1032     gmx_bool    bRM;
1033
1034     snew(*hb, nres);
1035     snew(*restp, nres);
1036     /* first the termini */
1037     for (i = 0; i < nterpairs; i++)
1038     {
1039         if (rn[i] >= 0 && ntdb[i] != nullptr)
1040         {
1041             copy_t_hackblock(ntdb[i], &(*hb)[rn[i]]);
1042         }
1043         if (rc[i] >= 0 && ctdb[i] != nullptr)
1044         {
1045             merge_t_hackblock(ctdb[i], &(*hb)[rc[i]]);
1046         }
1047     }
1048
1049     /* then the whole rtp */
1050     for (i = 0; i < nres; i++)
1051     {
1052         /* Here we allow a mismatch of one character when looking for the rtp entry.
1053          * For such a mismatch there should be only one mismatching name.
1054          * This is mainly useful for small molecules such as ions.
1055          * Note that this will usually not work for protein, DNA and RNA,
1056          * since there the residue names should be listed in residuetypes.dat
1057          * and an error will have been generated earlier in the process.
1058          */
1059         key = *resinfo[i].rtp;
1060         snew(resinfo[i].rtp, 1);
1061         *resinfo[i].rtp = search_rtp(key, nrtp, rtp);
1062         res             = get_restp(*resinfo[i].rtp, nrtp, rtp);
1063         copy_t_restp(res, &(*restp)[i]);
1064
1065         /* Check that we do not have different bonded types in one molecule */
1066         check_restp_types(&(*restp)[0], &(*restp)[i]);
1067
1068         tern = -1;
1069         for (j = 0; j < nterpairs && tern == -1; j++)
1070         {
1071             if (i == rn[j])
1072             {
1073                 tern = j;
1074             }
1075         }
1076         terc = -1;
1077         for (j = 0; j < nterpairs && terc == -1; j++)
1078         {
1079             if (i == rc[j])
1080             {
1081                 terc = j;
1082             }
1083         }
1084         bRM = merge_t_bondeds(res->rb, (*hb)[i].rb, tern >= 0, terc >= 0);
1085
1086         if (bRM && ((tern >= 0 && ntdb[tern] == nullptr) ||
1087                     (terc >= 0 && ctdb[terc] == nullptr)))
1088         {
1089             const char *errString = "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).";
1090             if (bAllowMissing)
1091             {
1092                 fprintf(stderr, "%s\n", errString);
1093             }
1094             else
1095             {
1096                 gmx_fatal(FARGS, errString);
1097             }
1098         }
1099         else if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
1100                          (terc >= 0 && ctdb[terc]->nhack == 0)))
1101         {
1102             const char *errString = "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             if (bAllowMissing)
1104             {
1105                 fprintf(stderr, "%s\n", errString);
1106             }
1107             else
1108             {
1109                 gmx_fatal(FARGS, errString);
1110             }
1111         }
1112     }
1113
1114     /* now perform t_hack's on t_restp's,
1115        i.e. add's and deletes from termini database will be
1116        added to/removed from residue topology
1117        we'll do this on one big dirty loop, so it won't make easy reading! */
1118     for (i = 0; i < nres; i++)
1119     {
1120         for (j = 0; j < (*hb)[i].nhack; j++)
1121         {
1122             if ( (*hb)[i].hack[j].nr)
1123             {
1124                 /* find atom in restp */
1125                 for (l = 0; l < (*restp)[i].natom; l++)
1126                 {
1127                     if ( ( (*hb)[i].hack[j].oname == nullptr &&
1128                            strcmp((*hb)[i].hack[j].a[0], *(*restp)[i].atomname[l]) == 0 ) ||
1129                          ( (*hb)[i].hack[j].oname != nullptr &&
1130                            strcmp((*hb)[i].hack[j].oname, *(*restp)[i].atomname[l]) == 0 ) )
1131                     {
1132                         break;
1133                     }
1134                 }
1135                 if (l == (*restp)[i].natom)
1136                 {
1137                     /* If we are doing an atom rename only, we don't need
1138                      * to generate a fatal error if the old name is not found
1139                      * in the rtp.
1140                      */
1141                     /* Deleting can happen also only on the input atoms,
1142                      * not necessarily always on the rtp entry.
1143                      */
1144                     if (!((*hb)[i].hack[j].oname != nullptr &&
1145                           (*hb)[i].hack[j].nname != nullptr) &&
1146                         !((*hb)[i].hack[j].oname != nullptr &&
1147                           (*hb)[i].hack[j].nname == nullptr))
1148                     {
1149                         gmx_fatal(FARGS,
1150                                   "atom %s not found in buiding block %d%s "
1151                                   "while combining tdb and rtp",
1152                                   (*hb)[i].hack[j].oname != nullptr ?
1153                                   (*hb)[i].hack[j].oname : (*hb)[i].hack[j].a[0],
1154                                   i+1, *resinfo[i].rtp);
1155                     }
1156                 }
1157                 else
1158                 {
1159                     if ( (*hb)[i].hack[j].oname == nullptr)
1160                     {
1161                         /* we're adding: */
1162                         add_atom_to_restp(&(*restp)[i], l, &(*hb)[i].hack[j]);
1163                     }
1164                     else
1165                     {
1166                         /* oname != NULL */
1167                         if ( (*hb)[i].hack[j].nname == nullptr)
1168                         {
1169                             /* we're deleting */
1170                             if (debug)
1171                             {
1172                                 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
1173                                         *(*restp)[i].atomname[l],
1174                                         i+1, (*restp)[i].resname);
1175                             }
1176                             /* shift the rest */
1177                             (*restp)[i].natom--;
1178                             for (k = l; k < (*restp)[i].natom; k++)
1179                             {
1180                                 (*restp)[i].atom    [k] = (*restp)[i].atom    [k+1];
1181                                 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
1182                                 (*restp)[i].cgnr    [k] = (*restp)[i].cgnr    [k+1];
1183                             }
1184                             /* give back space */
1185                             srenew((*restp)[i].atom,     (*restp)[i].natom);
1186                             srenew((*restp)[i].atomname, (*restp)[i].natom);
1187                             srenew((*restp)[i].cgnr,     (*restp)[i].natom);
1188                         }
1189                         else /* nname != NULL */
1190                         {    /* we're replacing */
1191                             if (debug)
1192                             {
1193                                 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
1194                                         *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
1195                                         i+1, (*restp)[i].resname);
1196                             }
1197                             snew( (*restp)[i].atomname[l], 1);
1198                             (*restp)[i].atom[l]      =       *(*hb)[i].hack[j].atom;
1199                             *(*restp)[i].atomname[l] = gmx_strdup((*hb)[i].hack[j].nname);
1200                             if ( (*hb)[i].hack[j].cgnr != NOTSET)
1201                             {
1202                                 (*restp)[i].cgnr   [l] = (*hb)[i].hack[j].cgnr;
1203                             }
1204                         }
1205                     }
1206                 }
1207             }
1208         }
1209     }
1210 }
1211
1212 static gmx_bool atomname_cmp_nr(const char *anm, t_hack *hack, int *nr)
1213 {
1214
1215     if (hack->nr == 1)
1216     {
1217         *nr = 0;
1218
1219         return (gmx_strcasecmp(anm, hack->nname) == 0);
1220     }
1221     else
1222     {
1223         if (isdigit(anm[strlen(anm)-1]))
1224         {
1225             *nr = anm[strlen(anm)-1] - '0';
1226         }
1227         else
1228         {
1229             *nr = 0;
1230         }
1231         if (*nr <= 0 || *nr > hack->nr)
1232         {
1233             return FALSE;
1234         }
1235         else
1236         {
1237             return (strlen(anm) == strlen(hack->nname) + 1 &&
1238                     gmx_strncasecmp(anm, hack->nname, strlen(hack->nname)) == 0);
1239         }
1240     }
1241 }
1242
1243 static gmx_bool match_atomnames_with_rtp_atom(t_atoms *pdba, rvec *x, int atind,
1244                                               t_restp *rptr, t_hackblock *hbr,
1245                                               gmx_bool bVerbose)
1246 {
1247     int      resnr;
1248     int      j, k;
1249     char    *oldnm, *newnm;
1250     int      anmnr;
1251     char    *start_at, buf[STRLEN];
1252     int      start_nr;
1253     gmx_bool bReplaceReplace, bFoundInAdd;
1254     gmx_bool bDeleted;
1255
1256     oldnm = *pdba->atomname[atind];
1257     resnr = pdba->resinfo[pdba->atom[atind].resind].nr;
1258
1259     bDeleted = FALSE;
1260     for (j = 0; j < hbr->nhack; j++)
1261     {
1262         if (hbr->hack[j].oname != nullptr && hbr->hack[j].nname != nullptr &&
1263             gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1264         {
1265             /* This is a replace entry. */
1266             /* Check if we are not replacing a replaced atom. */
1267             bReplaceReplace = FALSE;
1268             for (k = 0; k < hbr->nhack; k++)
1269             {
1270                 if (k != j &&
1271                     hbr->hack[k].oname != nullptr && hbr->hack[k].nname != nullptr &&
1272                     gmx_strcasecmp(hbr->hack[k].nname, hbr->hack[j].oname) == 0)
1273                 {
1274                     /* The replace in hack[j] replaces an atom that
1275                      * was already replaced in hack[k], we do not want
1276                      * second or higher level replaces at this stage.
1277                      */
1278                     bReplaceReplace = TRUE;
1279                 }
1280             }
1281             if (bReplaceReplace)
1282             {
1283                 /* Skip this replace. */
1284                 continue;
1285             }
1286
1287             /* This atom still has the old name, rename it */
1288             newnm = hbr->hack[j].nname;
1289             for (k = 0; k < rptr->natom; k++)
1290             {
1291                 if (gmx_strcasecmp(newnm, *rptr->atomname[k]) == 0)
1292                 {
1293                     break;
1294                 }
1295             }
1296             if (k == rptr->natom)
1297             {
1298                 /* The new name is not present in the rtp.
1299                  * We need to apply the replace also to the rtp entry.
1300                  */
1301
1302                 /* We need to find the add hack that can add this atom
1303                  * to find out after which atom it should be added.
1304                  */
1305                 bFoundInAdd = FALSE;
1306                 for (k = 0; k < hbr->nhack; k++)
1307                 {
1308                     if (hbr->hack[k].oname == nullptr &&
1309                         hbr->hack[k].nname != nullptr &&
1310                         atomname_cmp_nr(newnm, &hbr->hack[k], &anmnr))
1311                     {
1312                         if (anmnr <= 1)
1313                         {
1314                             start_at = hbr->hack[k].a[0];
1315                         }
1316                         else
1317                         {
1318                             sprintf(buf, "%s%d", hbr->hack[k].nname, anmnr-1);
1319                             start_at = buf;
1320                         }
1321                         for (start_nr = 0; start_nr < rptr->natom; start_nr++)
1322                         {
1323                             if (gmx_strcasecmp(start_at, (*rptr->atomname[start_nr])) == 0)
1324                             {
1325                                 break;
1326                             }
1327                         }
1328                         if (start_nr == rptr->natom)
1329                         {
1330                             gmx_fatal(FARGS, "Could not find atom '%s' in residue building block '%s' to add atom '%s' to",
1331                                       start_at, rptr->resname, newnm);
1332                         }
1333                         /* We can add the atom after atom start_nr */
1334                         add_atom_to_restp(rptr, start_nr, &hbr->hack[j]);
1335
1336                         bFoundInAdd = TRUE;
1337                     }
1338                 }
1339
1340                 if (!bFoundInAdd)
1341                 {
1342                     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'",
1343                               newnm,
1344                               hbr->hack[j].oname, hbr->hack[j].nname,
1345                               rptr->resname);
1346                 }
1347             }
1348
1349             if (bVerbose)
1350             {
1351                 printf("Renaming atom '%s' in residue '%s' %d to '%s'\n",
1352                        oldnm, rptr->resname, resnr, newnm);
1353             }
1354             /* Rename the atom in pdba */
1355             snew(pdba->atomname[atind], 1);
1356             *pdba->atomname[atind] = gmx_strdup(newnm);
1357         }
1358         else if (hbr->hack[j].oname != nullptr && hbr->hack[j].nname == nullptr &&
1359                  gmx_strcasecmp(oldnm, hbr->hack[j].oname) == 0)
1360         {
1361             /* This is a delete entry, check if this atom is present
1362              * in the rtp entry of this residue.
1363              */
1364             for (k = 0; k < rptr->natom; k++)
1365             {
1366                 if (gmx_strcasecmp(oldnm, *rptr->atomname[k]) == 0)
1367                 {
1368                     break;
1369                 }
1370             }
1371             if (k == rptr->natom)
1372             {
1373                 /* This atom is not present in the rtp entry,
1374                  * delete is now from the input pdba.
1375                  */
1376                 if (bVerbose)
1377                 {
1378                     printf("Deleting atom '%s' in residue '%s' %d\n",
1379                            oldnm, rptr->resname, resnr);
1380                 }
1381                 /* We should free the atom name,
1382                  * but it might be used multiple times in the symtab.
1383                  * sfree(pdba->atomname[atind]);
1384                  */
1385                 for (k = atind+1; k < pdba->nr; k++)
1386                 {
1387                     pdba->atom[k-1]     = pdba->atom[k];
1388                     pdba->atomname[k-1] = pdba->atomname[k];
1389                     copy_rvec(x[k], x[k-1]);
1390                 }
1391                 pdba->nr--;
1392                 bDeleted = TRUE;
1393             }
1394         }
1395     }
1396
1397     return bDeleted;
1398 }
1399
1400 void match_atomnames_with_rtp(t_restp restp[], t_hackblock hb[],
1401                               t_atoms *pdba, rvec *x,
1402                               gmx_bool bVerbose)
1403 {
1404     int          i, j;
1405     char        *oldnm;
1406     t_restp     *rptr;
1407
1408     for (i = 0; i < pdba->nr; i++)
1409     {
1410         oldnm = *pdba->atomname[i];
1411         rptr  = &restp[pdba->atom[i].resind];
1412         for (j = 0; (j < rptr->natom); j++)
1413         {
1414             if (gmx_strcasecmp(oldnm, *(rptr->atomname[j])) == 0)
1415             {
1416                 break;
1417             }
1418         }
1419         if (j == rptr->natom)
1420         {
1421             /* Not found yet, check if we have to rename this atom */
1422             if (match_atomnames_with_rtp_atom(pdba, x, i,
1423                                               rptr, &(hb[pdba->atom[i].resind]),
1424                                               bVerbose))
1425             {
1426                 /* We deleted this atom, decrease the atom counter by 1. */
1427                 i--;
1428             }
1429         }
1430     }
1431 }
1432
1433 #define NUM_CMAP_ATOMS 5
1434 static void gen_cmap(t_params *psb, t_restp *restp, t_atoms *atoms)
1435 {
1436     int         residx, i, j, k;
1437     const char *ptr;
1438     const char *pname;
1439     t_resinfo  *resinfo = atoms->resinfo;
1440     int         nres    = atoms->nres;
1441     gmx_bool    bAddCMAP;
1442     int         cmap_atomid[NUM_CMAP_ATOMS];
1443     int         cmap_chainnum = -1, this_residue_index;
1444
1445     if (debug)
1446     {
1447         ptr = "cmap";
1448     }
1449     else
1450     {
1451         ptr = "check";
1452     }
1453
1454     fprintf(stderr, "Making cmap torsions...\n");
1455     i = 0;
1456     /* Most cmap entries use the N atom from the next residue, so the last
1457      * residue should not have its CMAP entry in that case, but for things like
1458      * dipeptides we sometimes define a complete CMAP entry inside a residue,
1459      * and in this case we need to process everything through the last residue.
1460      */
1461     for (residx = 0; residx < nres; residx++)
1462     {
1463         /* Add CMAP terms from the list of CMAP interactions */
1464         for (j = 0; j < restp[residx].rb[ebtsCMAP].nb; j++)
1465         {
1466             bAddCMAP = TRUE;
1467             /* Loop over atoms in a candidate CMAP interaction and
1468              * check that they exist, are from the same chain and are
1469              * from residues labelled as protein. */
1470             for (k = 0; k < NUM_CMAP_ATOMS && bAddCMAP; k++)
1471             {
1472                 /* Assign the pointer to the name of the next reference atom.
1473                  * This can use -/+ labels to refer to previous/next residue.
1474                  */
1475                 pname = restp[residx].rb[ebtsCMAP].b[j].a[k];
1476                 /* Skip this CMAP entry if it refers to residues before the
1477                  * first or after the last residue.
1478                  */
1479                 if (((strchr(pname, '-') != nullptr) && (residx == 0)) ||
1480                     ((strchr(pname, '+') != nullptr) && (residx == nres-1)))
1481                 {
1482                     bAddCMAP = FALSE;
1483                     break;
1484                 }
1485
1486                 cmap_atomid[k] = search_atom(pname,
1487                                              i, atoms, ptr, TRUE);
1488                 bAddCMAP = bAddCMAP && (cmap_atomid[k] != -1);
1489                 if (!bAddCMAP)
1490                 {
1491                     /* This break is necessary, because cmap_atomid[k]
1492                      * == -1 cannot be safely used as an index
1493                      * into the atom array. */
1494                     break;
1495                 }
1496                 this_residue_index = atoms->atom[cmap_atomid[k]].resind;
1497                 if (0 == k)
1498                 {
1499                     cmap_chainnum = resinfo[this_residue_index].chainnum;
1500                 }
1501                 else
1502                 {
1503                     /* Does the residue for this atom have the same
1504                      * chain number as the residues for previous
1505                      * atoms? */
1506                     bAddCMAP = bAddCMAP &&
1507                         cmap_chainnum == resinfo[this_residue_index].chainnum;
1508                 }
1509                 /* Here we used to check that the residuetype was protein and
1510                  * disable bAddCMAP if that was not the case. However, some
1511                  * special residues (say, alanine dipeptides) might not adhere
1512                  * to standard naming, and if we start calling them normal
1513                  * protein residues the user will be bugged to select termini.
1514                  *
1515                  * Instead, I believe that the right course of action is to
1516                  * keep the CMAP interaction if it is present in the RTP file
1517                  * and we correctly identified all atoms (which is the case
1518                  * if we got here).
1519                  */
1520             }
1521
1522             if (bAddCMAP)
1523             {
1524                 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);
1525             }
1526         }
1527
1528         if (residx < nres-1)
1529         {
1530             while (atoms->atom[i].resind < residx+1)
1531             {
1532                 i++;
1533             }
1534         }
1535     }
1536     /* Start the next residue */
1537 }
1538
1539 static void
1540 scrub_charge_groups(int *cgnr, int natoms)
1541 {
1542     int i;
1543
1544     for (i = 0; i < natoms; i++)
1545     {
1546         cgnr[i] = i+1;
1547     }
1548 }
1549
1550
1551 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
1552              t_atoms *atoms, rvec **x, gpp_atomtype_t atype, t_symtab *tab,
1553              int nrtp, t_restp rtp[],
1554              t_restp *restp, t_hackblock *hb,
1555              gmx_bool bAllowMissing,
1556              gmx_bool bVsites, gmx_bool bVsiteAromatics,
1557              const char *ffdir,
1558              real mHmult,
1559              int nssbonds, t_ssbond *ssbonds,
1560              real long_bond_dist, real short_bond_dist,
1561              gmx_bool bDeuterate, gmx_bool bChargeGroups, gmx_bool bCmap,
1562              gmx_bool bRenumRes, gmx_bool bRTPresname)
1563 {
1564     /*
1565        t_hackblock *hb;
1566        t_restp  *restp;
1567      */
1568     t_params          plist[F_NRE];
1569     t_excls          *excls;
1570     t_nextnb          nnb;
1571     int              *cgnr;
1572     int              *vsite_type;
1573     int               i, nmissat;
1574     int               bts[ebtsNR];
1575     gmx_residuetype_t*rt;
1576
1577     init_plist(plist);
1578     gmx_residuetype_init(&rt);
1579
1580     if (debug)
1581     {
1582         print_resall(debug, atoms->nres, restp, atype);
1583         dump_hb(debug, atoms->nres, hb);
1584     }
1585
1586     /* Make bonds */
1587     at2bonds(&(plist[F_BONDS]), hb,
1588              atoms, *x,
1589              long_bond_dist, short_bond_dist);
1590
1591     /* specbonds: disulphide bonds & heme-his */
1592     do_ssbonds(&(plist[F_BONDS]),
1593                atoms, nssbonds, ssbonds,
1594                bAllowMissing);
1595
1596     nmissat = name2type(atoms, &cgnr, restp, rt);
1597     if (nmissat)
1598     {
1599         if (bAllowMissing)
1600         {
1601             fprintf(stderr, "There were %d missing atoms in molecule %s\n",
1602                     nmissat, molname);
1603         }
1604         else
1605         {
1606             gmx_fatal(FARGS, "There were %d missing atoms in molecule %s, if you want to use this incomplete topology anyhow, use the option -missing",
1607                       nmissat, molname);
1608         }
1609     }
1610
1611     /* Cleanup bonds (sort and rm doubles) */
1612     clean_bonds(&(plist[F_BONDS]));
1613
1614     snew(vsite_type, atoms->nr);
1615     for (i = 0; i < atoms->nr; i++)
1616     {
1617         vsite_type[i] = NOTSET;
1618     }
1619     if (bVsites)
1620     {
1621         /* determine which atoms will be vsites and add dummy masses
1622            also renumber atom numbers in plist[0..F_NRE]! */
1623         do_vsites(nrtp, rtp, atype, atoms, tab, x, plist,
1624                   &vsite_type, &cgnr, mHmult, bVsiteAromatics, ffdir);
1625     }
1626
1627     /* Make Angles and Dihedrals */
1628     fprintf(stderr, "Generating angles, dihedrals and pairs...\n");
1629     snew(excls, atoms->nr);
1630     init_nnb(&nnb, atoms->nr, 4);
1631     gen_nnb(&nnb, plist);
1632     print_nnb(&nnb, "NNB");
1633     gen_pad(&nnb, atoms, restp, plist, excls, hb, bAllowMissing);
1634     done_nnb(&nnb);
1635
1636     /* Make CMAP */
1637     if (TRUE == bCmap)
1638     {
1639         gen_cmap(&(plist[F_CMAP]), restp, atoms);
1640         if (plist[F_CMAP].nr > 0)
1641         {
1642             fprintf(stderr, "There are %4d cmap torsion pairs\n",
1643                     plist[F_CMAP].nr);
1644         }
1645     }
1646
1647     /* set mass of all remaining hydrogen atoms */
1648     if (mHmult != 1.0)
1649     {
1650         do_h_mass(&(plist[F_BONDS]), vsite_type, atoms, mHmult, bDeuterate);
1651     }
1652     sfree(vsite_type);
1653
1654     /* Cleanup bonds (sort and rm doubles) */
1655     /* clean_bonds(&(plist[F_BONDS]));*/
1656
1657     fprintf(stderr,
1658             "There are %4d dihedrals, %4d impropers, %4d angles\n"
1659             "          %4d pairs,     %4d bonds and  %4d virtual sites\n",
1660             plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
1661             plist[F_LJ14].nr, plist[F_BONDS].nr,
1662             plist[F_VSITE2].nr +
1663             plist[F_VSITE3].nr +
1664             plist[F_VSITE3FD].nr +
1665             plist[F_VSITE3FAD].nr +
1666             plist[F_VSITE3OUT].nr +
1667             plist[F_VSITE4FD].nr +
1668             plist[F_VSITE4FDN].nr );
1669
1670     print_sums(atoms, FALSE);
1671
1672     if (FALSE == bChargeGroups)
1673     {
1674         scrub_charge_groups(cgnr, atoms->nr);
1675     }
1676
1677     if (bRenumRes)
1678     {
1679         for (i = 0; i < atoms->nres; i++)
1680         {
1681             atoms->resinfo[i].nr = i + 1;
1682             atoms->resinfo[i].ic = ' ';
1683         }
1684     }
1685
1686     if (top_file)
1687     {
1688         fprintf(stderr, "Writing topology\n");
1689         /* We can copy the bonded types from the first restp,
1690          * since the types have to be identical for all residues in one molecule.
1691          */
1692         for (i = 0; i < ebtsNR; i++)
1693         {
1694             bts[i] = restp[0].rb[i].type;
1695         }
1696         write_top(top_file, posre_fn, molname,
1697                   atoms, bRTPresname,
1698                   bts, plist, excls, atype, cgnr, restp[0].nrexcl);
1699     }
1700
1701     /* cleaning up */
1702     free_t_hackblock(atoms->nres, &hb);
1703     free_t_restp(atoms->nres, &restp);
1704     gmx_residuetype_destroy(rt);
1705
1706     /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
1707     sfree(cgnr);
1708     for (i = 0; i < F_NRE; i++)
1709     {
1710         sfree(plist[i].param);
1711     }
1712     for (i = 0; i < atoms->nr; i++)
1713     {
1714         sfree(excls[i].e);
1715     }
1716     sfree(excls);
1717 }