2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
52 #include "gromacs/commandline/cmdlineoptionsmodule.h"
53 #include "gromacs/fileio/confio.h"
54 #include "gromacs/fileio/filetypes.h"
55 #include "gromacs/fileio/gmxfio.h"
56 #include "gromacs/fileio/pdbio.h"
57 #include "gromacs/gmxlib/conformation_utilities.h"
58 #include "gromacs/gmxpreprocess/fflibutil.h"
59 #include "gromacs/gmxpreprocess/genhydro.h"
60 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
61 #include "gromacs/gmxpreprocess/grompp_impl.h"
62 #include "gromacs/gmxpreprocess/h_db.h"
63 #include "gromacs/gmxpreprocess/hizzie.h"
64 #include "gromacs/gmxpreprocess/pdb2top.h"
65 #include "gromacs/gmxpreprocess/pgutil.h"
66 #include "gromacs/gmxpreprocess/specbond.h"
67 #include "gromacs/gmxpreprocess/ter_db.h"
68 #include "gromacs/gmxpreprocess/toputil.h"
69 #include "gromacs/gmxpreprocess/xlate.h"
70 #include "gromacs/math/vec.h"
71 #include "gromacs/options/basicoptions.h"
72 #include "gromacs/options/filenameoption.h"
73 #include "gromacs/options/ioptionscontainer.h"
74 #include "gromacs/topology/atomprop.h"
75 #include "gromacs/topology/block.h"
76 #include "gromacs/topology/index.h"
77 #include "gromacs/topology/residuetypes.h"
78 #include "gromacs/topology/symtab.h"
79 #include "gromacs/topology/topology.h"
80 #include "gromacs/utility/dir_separator.h"
81 #include "gromacs/utility/exceptions.h"
82 #include "gromacs/utility/fatalerror.h"
83 #include "gromacs/utility/filestream.h"
84 #include "gromacs/utility/loggerbuilder.h"
85 #include "gromacs/utility/path.h"
86 #include "gromacs/utility/smalloc.h"
87 #include "gromacs/utility/strdb.h"
88 #include "gromacs/utility/stringutil.h"
90 #include "hackblock.h"
95 RtpRename(const char* newGmx, const char* newMain, const char* newNter, const char* newCter, const char* newBter) :
110 static const char* res2bb_notermini(const std::string& name, gmx::ArrayRef<const RtpRename> rr)
112 /* NOTE: This function returns the main building block name,
113 * it does not take terminal renaming into account.
115 auto found = std::find_if(rr.begin(), rr.end(), [&name](const auto& rename) {
116 return gmx::equalCaseInsensitive(name, rename.gmx);
118 return found != rr.end() ? found->main.c_str() : name.c_str();
121 static const char* select_res(int nr,
126 gmx::ArrayRef<const RtpRename> rr)
128 printf("Which %s type do you want for residue %d\n", title, resnr + 1);
129 for (int sel = 0; (sel < nr); sel++)
131 printf("%d. %s (%s)\n", sel, expl[sel], res2bb_notermini(name[sel], rr));
133 printf("\nType a number:");
137 if (scanf("%d", &userSelection) != 1)
139 gmx_fatal(FARGS, "Answer me for res %s %d!", title, resnr + 1);
142 return name[userSelection];
145 static const char* get_asptp(int resnr, gmx::ArrayRef<const RtpRename> rr)
153 const char* lh[easpNR] = { "ASP", "ASPH" };
154 const char* expl[easpNR] = { "Not protonated (charge -1)", "Protonated (charge 0)" };
156 return select_res(easpNR, resnr, lh, expl, "ASPARTIC ACID", rr);
159 static const char* get_glutp(int resnr, gmx::ArrayRef<const RtpRename> rr)
167 const char* lh[egluNR] = { "GLU", "GLUH" };
168 const char* expl[egluNR] = { "Not protonated (charge -1)", "Protonated (charge 0)" };
170 return select_res(egluNR, resnr, lh, expl, "GLUTAMIC ACID", rr);
173 static const char* get_glntp(int resnr, gmx::ArrayRef<const RtpRename> rr)
181 const char* lh[eglnNR] = { "GLN", "QLN" };
182 const char* expl[eglnNR] = { "Not protonated (charge 0)", "Protonated (charge +1)" };
184 return select_res(eglnNR, resnr, lh, expl, "GLUTAMINE", rr);
187 static const char* get_lystp(int resnr, gmx::ArrayRef<const RtpRename> rr)
195 const char* lh[elysNR] = { "LYSN", "LYS" };
196 const char* expl[elysNR] = { "Not protonated (charge 0)", "Protonated (charge +1)" };
198 return select_res(elysNR, resnr, lh, expl, "LYSINE", rr);
201 static const char* get_argtp(int resnr, gmx::ArrayRef<const RtpRename> rr)
209 const char* lh[eargNR] = { "ARGN", "ARG" };
210 const char* expl[eargNR] = { "Not protonated (charge 0)", "Protonated (charge +1)" };
212 return select_res(eargNR, resnr, lh, expl, "ARGININE", rr);
215 static const char* get_histp(int resnr, gmx::ArrayRef<const RtpRename> rr)
217 const char* expl[ehisNR] = { "H on ND1 only", "H on NE2 only", "H on ND1 and NE2",
220 return select_res(ehisNR, resnr, hh, expl, "HISTIDINE", rr);
223 static void read_rtprename(const char* fname, FILE* fp, std::vector<RtpRename>* rtprename)
225 char line[STRLEN], buf[STRLEN];
228 while (get_a_line(fp, line, STRLEN))
230 /* line is NULL-terminated and length<STRLEN, so final arg cannot overflow.
231 * For other args, we read up to 6 chars (so we can detect if the length is > 5).
232 * Note that the buffer length has been increased to 7 to allow this,
233 * so we just need to make sure the strings have zero-length initially.
245 int nc = sscanf(line, "%6s %6s %6s %6s %6s %s", gmx, main, nter, cter, bter, buf);
246 RtpRename newEntry(gmx, main, nter, cter, bter);
249 if (nc != 2 && nc != 5)
251 gmx_fatal(FARGS, "Residue renaming database '%s' has %d columns instead of %d or %d",
259 "A line in residue renaming database '%s' has %d columns, while previous "
260 "lines have %d columns",
266 /* This file does not have special termini names, copy them from main */
267 newEntry.nter = newEntry.main;
268 newEntry.cter = newEntry.main;
269 newEntry.bter = newEntry.main;
271 rtprename->push_back(newEntry);
275 static std::string search_resrename(gmx::ArrayRef<const RtpRename> rr,
279 bool bCompareFFRTPname)
281 auto found = std::find_if(rr.begin(), rr.end(), [&name, &bCompareFFRTPname](const auto& rename) {
282 return ((!bCompareFFRTPname && (name == rename.gmx))
283 || (bCompareFFRTPname && (name == rename.main)));
287 /* If found in the database, rename this residue's rtp building block,
288 * otherwise keep the old name.
290 if (found != rr.end())
294 newName = found->bter;
298 newName = found->nter;
302 newName = found->cter;
306 newName = found->main;
309 if (newName[0] == '-')
311 gmx_fatal(FARGS, "In the chosen force field there is no residue type for '%s'%s", name,
312 bStart ? (bEnd ? " as a standalone (starting & ending) residue" : " as a starting terminus")
313 : (bEnd ? " as an ending terminus" : ""));
320 static void rename_resrtp(t_atoms* pdba,
322 gmx::ArrayRef<const int> r_start,
323 gmx::ArrayRef<const int> r_end,
324 gmx::ArrayRef<const RtpRename> rr,
327 const gmx::MDLogger& logger)
329 bool bFFRTPTERRNM = (getenv("GMX_NO_FFRTP_TER_RENAME") == nullptr);
331 for (int r = 0; r < pdba->nres; r++)
335 for (int j = 0; j < nterpairs; j++)
342 for (int j = 0; j < nterpairs; j++)
350 std::string newName = search_resrename(rr, *pdba->resinfo[r].rtp, bStart, bEnd, false);
352 if (bFFRTPTERRNM && newName.empty() && (bStart || bEnd))
354 /* This is a terminal residue, but the residue name,
355 * currently stored in .rtp, is not a standard residue name,
356 * but probably a force field specific rtp name.
357 * Check if we need to rename it because it is terminal.
359 newName = search_resrename(rr, *pdba->resinfo[r].rtp, bStart, bEnd, true);
362 if (!newName.empty() && newName != *pdba->resinfo[r].rtp)
368 .appendTextFormatted("Changing rtp entry of residue %d %s to '%s'",
369 pdba->resinfo[r].nr, *pdba->resinfo[r].name,
372 pdba->resinfo[r].rtp = put_symtab(symtab, newName.c_str());
377 static void pdbres_to_gmxrtp(t_atoms* pdba)
381 for (i = 0; (i < pdba->nres); i++)
383 if (pdba->resinfo[i].rtp == nullptr)
385 pdba->resinfo[i].rtp = pdba->resinfo[i].name;
390 static void rename_pdbres(t_atoms* pdba, const char* oldnm, const char* newnm, bool bFullCompare, t_symtab* symtab)
395 for (i = 0; (i < pdba->nres); i++)
397 resnm = *pdba->resinfo[i].name;
398 if ((bFullCompare && (gmx::equalCaseInsensitive(resnm, oldnm)))
399 || (!bFullCompare && strstr(resnm, oldnm) != nullptr))
401 /* Rename the residue name (not the rtp name) */
402 pdba->resinfo[i].name = put_symtab(symtab, newnm);
407 static void rename_bb(t_atoms* pdba, const char* oldnm, const char* newnm, bool bFullCompare, t_symtab* symtab)
412 for (i = 0; (i < pdba->nres); i++)
414 /* We have not set the rtp name yes, use the residue name */
415 bbnm = *pdba->resinfo[i].name;
416 if ((bFullCompare && (gmx::equalCaseInsensitive(bbnm, oldnm)))
417 || (!bFullCompare && strstr(bbnm, oldnm) != nullptr))
419 /* Change the rtp builing block name */
420 pdba->resinfo[i].rtp = put_symtab(symtab, newnm);
425 static void rename_bbint(t_atoms* pdba,
427 const char* gettp(int, gmx::ArrayRef<const RtpRename>),
430 gmx::ArrayRef<const RtpRename> rr)
436 for (i = 0; i < pdba->nres; i++)
438 /* We have not set the rtp name yet, use the residue name */
439 bbnm = *pdba->resinfo[i].name;
440 if ((bFullCompare && (strcmp(bbnm, oldnm) == 0)) || (!bFullCompare && strstr(bbnm, oldnm) != nullptr))
443 pdba->resinfo[i].rtp = put_symtab(symtab, ptr);
448 static void check_occupancy(t_atoms* atoms, const char* filename, bool bVerbose, const gmx::MDLogger& logger)
454 ftp = fn2ftp(filename);
455 if (!atoms->pdbinfo || ((ftp != efPDB) && (ftp != efBRK) && (ftp != efENT)))
457 GMX_LOG(logger.warning).asParagraph().appendTextFormatted("No occupancies in %s", filename);
461 for (i = 0; (i < atoms->nr); i++)
463 if (atoms->pdbinfo[i].occup != 1)
467 GMX_LOG(logger.warning)
469 .appendTextFormatted("Occupancy for atom %s%d-%s is %f rather than 1",
470 *atoms->resinfo[atoms->atom[i].resind].name,
471 atoms->resinfo[atoms->atom[i].resind].nr,
472 *atoms->atomname[i], atoms->pdbinfo[i].occup);
474 if (atoms->pdbinfo[i].occup == 0)
484 if (nzero == atoms->nr)
486 GMX_LOG(logger.warning)
488 .appendTextFormatted(
489 "All occupancy fields zero. This is probably not an X-Ray structure");
491 else if ((nzero > 0) || (nnotone > 0))
493 GMX_LOG(logger.warning)
495 .appendTextFormatted(
496 "there were %d atoms with zero occupancy and %d atoms with "
497 " occupancy unequal to one (out of %d atoms). Check your pdb "
499 nzero, nnotone, atoms->nr);
503 GMX_LOG(logger.warning).asParagraph().appendTextFormatted("All occupancies are one");
508 static void write_posres(const char* fn, t_atoms* pdba, real fc)
513 fp = gmx_fio_fopen(fn, "w");
515 "; In this topology include file, you will find position restraint\n"
516 "; entries for all the heavy atoms in your original pdb file.\n"
517 "; This means that all the protons which were added by pdb2gmx are\n"
518 "; not restrained.\n"
520 "[ position_restraints ]\n"
521 "; %4s%6s%8s%8s%8s\n",
522 "atom", "type", "fx", "fy", "fz");
523 for (i = 0; (i < pdba->nr); i++)
525 if (!is_hydrogen(*pdba->atomname[i]) && !is_dummymass(*pdba->atomname[i]))
527 fprintf(fp, "%6d%6d %g %g %g\n", i + 1, 1, fc, fc, fc);
533 static int read_pdball(const char* inf,
547 /* Read a pdb file. (containing proteins) */
549 int natom, new_natom, i;
552 printf("Reading %s...\n", inf);
553 readConfAndAtoms(inf, symtab, title, atoms, pbcType, x, nullptr, box);
555 if (atoms->pdbinfo == nullptr)
557 snew(atoms->pdbinfo, atoms->nr);
559 if (fn2ftp(inf) == efPDB)
561 get_pdb_atomnumber(atoms, aps);
566 for (i = 0; i < atoms->nr; i++)
568 if (!is_hydrogen(*atoms->atomname[i]))
570 atoms->atom[new_natom] = atoms->atom[i];
571 atoms->atomname[new_natom] = atoms->atomname[i];
572 atoms->pdbinfo[new_natom] = atoms->pdbinfo[i];
573 copy_rvec((*x)[i], (*x)[new_natom]);
577 atoms->nr = new_natom;
584 printf(" '%s',", *title);
586 printf(" %d atoms\n", natom);
588 /* Rename residues */
589 rename_pdbres(atoms, "HOH", watres, false, symtab);
590 rename_pdbres(atoms, "SOL", watres, false, symtab);
591 rename_pdbres(atoms, "WAT", watres, false, symtab);
593 rename_atoms("xlateat.dat", nullptr, atoms, symtab, {}, true, rt, true, bVerbose);
601 write_sto_conf(outf, *title, atoms, *x, nullptr, *pbcType, box);
607 static void process_chain(t_atoms* pdba,
608 gmx::ArrayRef<gmx::RVec> x,
621 gmx::ArrayRef<const RtpRename> rr)
623 /* Rename aromatics, lys, asp and histidine */
626 rename_bb(pdba, "TYR", "TYRU", false, symtab);
630 rename_bb(pdba, "TRP", "TRPU", false, symtab);
634 rename_bb(pdba, "PHE", "PHEU", false, symtab);
638 rename_bbint(pdba, "LYS", get_lystp, false, symtab, rr);
642 rename_bbint(pdba, "ARG", get_argtp, false, symtab, rr);
646 rename_bbint(pdba, "GLN", get_glntp, false, symtab, rr);
650 rename_bbint(pdba, "ASP", get_asptp, false, symtab, rr);
654 rename_bb(pdba, "ASPH", "ASP", false, symtab);
658 rename_bbint(pdba, "GLU", get_glutp, false, symtab, rr);
662 rename_bb(pdba, "GLUH", "GLU", false, symtab);
667 set_histp(pdba, gmx::as_rvec_array(x.data()), symtab, angle, distance);
671 rename_bbint(pdba, "HIS", get_histp, true, symtab, rr);
674 /* Initialize the rtp builing block names with the residue names
675 * for the residues that have not been processed above.
677 pdbres_to_gmxrtp(pdba);
679 /* Now we have all rtp names set.
680 * The rtp names will conform to Gromacs naming,
681 * unless the input pdb file contained one or more force field specific
682 * rtp names as residue names.
686 /* struct for sorting the atoms from the pdb file */
689 int resnr; /* residue number */
690 int j; /* database order index */
691 int index; /* original atom number */
692 char anm1; /* second letter of atom name */
693 char altloc; /* alternate location indicator */
696 static bool pdbicomp(const t_pdbindex& a, const t_pdbindex& b)
698 int d = (a.resnr - b.resnr);
704 d = (a.anm1 - b.anm1);
707 d = (a.altloc - b.altloc);
714 static void sort_pdbatoms(gmx::ArrayRef<const PreprocessResidue> restp_chain,
717 t_atoms** newPdbAtoms,
718 std::vector<gmx::RVec>* x,
722 t_atoms* pdba = *pdbaptr;
723 std::vector<gmx::RVec> xnew;
730 for (int i = 0; i < natoms; i++)
732 atomnm = *pdba->atomname[i];
733 const PreprocessResidue* localPpResidue = &restp_chain[pdba->atom[i].resind];
735 std::find_if(localPpResidue->atomname.begin(), localPpResidue->atomname.end(),
736 [&atomnm](char** it) { return gmx::equalCaseInsensitive(atomnm, *it); });
737 if (found == localPpResidue->atomname.end())
739 std::string buf = gmx::formatString(
740 "Atom %s in residue %s %d was not found in rtp entry %s with %d atoms\n"
741 "while sorting atoms.\n%s",
742 atomnm, *pdba->resinfo[pdba->atom[i].resind].name,
743 pdba->resinfo[pdba->atom[i].resind].nr, localPpResidue->resname.c_str(),
744 localPpResidue->natom(),
746 ? "\nFor a hydrogen, this can be a different protonation state, or it\n"
747 "might have had a different number in the PDB file and was rebuilt\n"
748 "(it might for instance have been H3, and we only expected H1 & "
750 "Note that hydrogens might have been added to the entry for the "
752 "Remove this hydrogen or choose a different protonation state to "
754 "Option -ignh will ignore all hydrogens in the input."
756 gmx_fatal(FARGS, "%s", buf.c_str());
758 /* make shadow array to be sorted into indexgroup */
759 pdbi[i].resnr = pdba->atom[i].resind;
760 pdbi[i].j = std::distance(localPpResidue->atomname.begin(), found);
762 pdbi[i].anm1 = atomnm[1];
763 pdbi[i].altloc = pdba->pdbinfo[i].altloc;
765 std::sort(pdbi, pdbi + natoms, pdbicomp);
767 /* pdba is sorted in pdbnew using the pdbi index */
768 std::vector<int> a(natoms);
769 srenew(*newPdbAtoms, 1);
770 init_t_atoms((*newPdbAtoms), natoms, true);
771 (*newPdbAtoms)->nr = pdba->nr;
772 (*newPdbAtoms)->nres = pdba->nres;
773 srenew((*newPdbAtoms)->resinfo, pdba->nres);
774 std::copy(pdba->resinfo, pdba->resinfo + pdba->nres, (*newPdbAtoms)->resinfo);
775 for (int i = 0; i < natoms; i++)
777 (*newPdbAtoms)->atom[i] = pdba->atom[pdbi[i].index];
778 (*newPdbAtoms)->atomname[i] = pdba->atomname[pdbi[i].index];
779 (*newPdbAtoms)->pdbinfo[i] = pdba->pdbinfo[pdbi[i].index];
780 xnew.push_back(x->at(pdbi[i].index));
781 /* make indexgroup in block */
782 a[i] = pdbi[i].index;
787 /* copy the sorted pdbnew back to pdba */
788 *pdbaptr = *newPdbAtoms;
790 add_grp(block, gnames, a, "prot_sort");
794 static int remove_duplicate_atoms(t_atoms* pdba, gmx::ArrayRef<gmx::RVec> x, bool bVerbose, const gmx::MDLogger& logger)
796 int i, j, oldnatoms, ndel;
799 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Checking for duplicate atoms....");
800 oldnatoms = pdba->nr;
802 /* NOTE: pdba->nr is modified inside the loop */
803 for (i = 1; (i < pdba->nr); i++)
805 /* compare 'i' and 'i-1', throw away 'i' if they are identical
806 this is a 'while' because multiple alternate locations can be present */
807 while ((i < pdba->nr) && (pdba->atom[i - 1].resind == pdba->atom[i].resind)
808 && (strcmp(*pdba->atomname[i - 1], *pdba->atomname[i]) == 0))
813 ri = &pdba->resinfo[pdba->atom[i].resind];
816 .appendTextFormatted("deleting duplicate atom %4s %s%4d%c",
817 *pdba->atomname[i], *ri->name, ri->nr, ri->ic);
818 if (ri->chainid && (ri->chainid != ' '))
820 printf(" ch %c", ri->chainid);
824 if (pdba->pdbinfo[i].atomnr)
826 printf(" pdb nr %4d", pdba->pdbinfo[i].atomnr);
828 if (pdba->pdbinfo[i].altloc && (pdba->pdbinfo[i].altloc != ' '))
830 printf(" altloc %c", pdba->pdbinfo[i].altloc);
836 /* We can not free, since it might be in the symtab */
837 /* sfree(pdba->atomname[i]); */
838 for (j = i; j < pdba->nr; j++)
840 pdba->atom[j] = pdba->atom[j + 1];
841 pdba->atomname[j] = pdba->atomname[j + 1];
844 pdba->pdbinfo[j] = pdba->pdbinfo[j + 1];
846 copy_rvec(x[j + 1], x[j]);
848 srenew(pdba->atom, pdba->nr);
849 /* srenew(pdba->atomname, pdba->nr); */
850 srenew(pdba->pdbinfo, pdba->nr);
853 if (pdba->nr != oldnatoms)
857 .appendTextFormatted("Now there are %d atoms. Deleted %d duplicates.", pdba->nr, ndel);
863 static void checkResidueTypeSanity(t_atoms* pdba, int r0, int r1, ResidueType* rt)
865 std::string startResidueString =
866 gmx::formatString("%s%d", *pdba->resinfo[r0].name, pdba->resinfo[r0].nr);
867 std::string endResidueString =
868 gmx::formatString("%s%d", *pdba->resinfo[r1 - 1].name, pdba->resinfo[r1 - 1].nr);
870 // Check whether all residues in chain have the same chain ID.
871 bool allResiduesHaveSameChainID = true;
872 char chainID0 = pdba->resinfo[r0].chainid;
874 std::string residueString;
876 for (int i = r0 + 1; i < r1; i++)
878 chainID = pdba->resinfo[i].chainid;
879 if (chainID != chainID0)
881 allResiduesHaveSameChainID = false;
882 residueString = gmx::formatString("%s%d", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
887 if (!allResiduesHaveSameChainID)
890 "The chain covering the range %s--%s does not have a consistent chain ID. "
891 "The first residue has ID '%c', while residue %s has ID '%c'.",
892 startResidueString.c_str(), endResidueString.c_str(), chainID0,
893 residueString.c_str(), chainID);
896 // At this point all residues have the same ID. If they are also non-blank
897 // we can be a bit more aggressive and require the types match too.
900 bool allResiduesHaveSameType = true;
902 std::string restype0 = rt->typeOfNamedDatabaseResidue(*pdba->resinfo[r0].name);
904 for (int i = r0 + 1; i < r1; i++)
906 restype = rt->typeOfNamedDatabaseResidue(*pdba->resinfo[i].name);
907 if (!gmx::equalCaseInsensitive(restype, restype0))
909 allResiduesHaveSameType = false;
910 residueString = gmx::formatString("%s%d", *pdba->resinfo[i].name, pdba->resinfo[i].nr);
915 if (!allResiduesHaveSameType)
918 "The residues in the chain %s--%s do not have a consistent type. "
919 "The first residue has type '%s', while residue %s is of type '%s'. "
920 "Either there is a mistake in your chain, or it includes nonstandard "
921 "residue names that have not yet been added to the residuetypes.dat "
922 "file in the GROMACS library directory. If there are other molecules "
923 "such as ligands, they should not have the same chain ID as the "
924 "adjacent protein chain since it's a separate molecule.",
925 startResidueString.c_str(), endResidueString.c_str(), restype0.c_str(),
926 residueString.c_str(), restype.c_str());
931 static void find_nc_ter(t_atoms* pdba, int r0, int r1, int* r_start, int* r_end, ResidueType* rt, const gmx::MDLogger& logger)
934 std::optional<std::string> startrestype;
939 int startWarnings = 0;
943 // Check that all residues have the same chain identifier, and if it is
944 // non-blank we also require the residue types to match.
945 checkResidueTypeSanity(pdba, r0, r1, rt);
947 // If we return correctly from checkResidueTypeSanity(), the only
948 // remaining cases where we can have non-matching residue types is if
949 // the chain ID was blank, which could be the case e.g. for a structure
950 // read from a GRO file or other file types without chain information.
951 // In that case we need to be a bit more liberal and detect chains based
952 // on the residue type.
954 // If we get here, the chain ID must be identical for all residues
955 char chainID = pdba->resinfo[r0].chainid;
957 /* Find the starting terminus (typially N or 5') */
958 for (i = r0; i < r1 && *r_start == -1; i++)
960 startrestype = rt->optionalTypeOfNamedDatabaseResidue(*pdba->resinfo[i].name);
965 if (gmx::equalCaseInsensitive(*startrestype, "Protein")
966 || gmx::equalCaseInsensitive(*startrestype, "DNA")
967 || gmx::equalCaseInsensitive(*startrestype, "RNA"))
971 .appendTextFormatted("Identified residue %s%d as a starting terminus.",
972 *pdba->resinfo[i].name, pdba->resinfo[i].nr);
975 else if (gmx::equalCaseInsensitive(*startrestype, "Ion"))
981 .appendTextFormatted(
982 "Residue %s%d has type 'Ion', assuming it is not linked into a "
984 *pdba->resinfo[i].name, pdba->resinfo[i].nr);
990 .appendTextFormatted("Disabling further notes about ions.");
996 // Either no known residue type, or one not needing special handling
997 if (startWarnings < 5)
1001 GMX_LOG(logger.warning)
1003 .appendTextFormatted(
1004 "Starting residue %s%d in chain not identified as "
1006 "This chain lacks identifiers, which makes it impossible to do "
1008 "classification of the start/end residues. Here we need to "
1009 "guess this residue "
1010 "should not be part of the chain and instead introduce a "
1011 "break, but that will "
1012 "be catastrophic if they should in fact be linked. Please "
1013 "check your structure, "
1014 "and add %s to residuetypes.dat if this was not correct.",
1015 *pdba->resinfo[i].name, pdba->resinfo[i].nr, *pdba->resinfo[i].name);
1019 GMX_LOG(logger.warning)
1021 .appendTextFormatted(
1022 "No residues in chain starting at %s%d identified as "
1024 "This makes it impossible to link them into a molecule, which "
1026 "correct or a catastrophic error. Please check your structure, "
1028 "necessary residue names to residuetypes.dat if this was not "
1030 *pdba->resinfo[i].name, pdba->resinfo[i].nr);
1033 if (startWarnings == 4)
1035 GMX_LOG(logger.warning)
1037 .appendTextFormatted(
1038 "Disabling further warnings about unidentified residues at start "
1047 /* Go through the rest of the residues, check that they are the same class, and identify the ending terminus. */
1048 for (int i = *r_start; i < r1; i++)
1050 std::optional<std::string> restype =
1051 rt->optionalTypeOfNamedDatabaseResidue(*pdba->resinfo[i].name);
1056 if (gmx::equalCaseInsensitive(*restype, *startrestype) && endWarnings == 0)
1060 else if (gmx::equalCaseInsensitive(*startrestype, "Ion"))
1064 GMX_LOG(logger.info)
1066 .appendTextFormatted(
1067 "Residue %s%d has type 'Ion', assuming it is not linked into a "
1069 *pdba->resinfo[i].name, pdba->resinfo[i].nr);
1073 GMX_LOG(logger.info)
1075 .appendTextFormatted("Disabling further notes about ions.");
1081 // Either no known residue type, or one not needing special handling.
1082 GMX_RELEASE_ASSERT(chainID == ' ', "Chain ID must be blank");
1083 // Otherwise the call to checkResidueTypeSanity() will
1084 // have caught the problem.
1085 if (endWarnings < 5)
1087 GMX_LOG(logger.warning)
1089 .appendTextFormatted(
1090 "Residue %s%d in chain has different type ('%s') from "
1091 "residue %s%d ('%s'). This chain lacks identifiers, which "
1093 "it impossible to do strict classification of the start/end "
1094 "residues. Here we "
1095 "need to guess this residue should not be part of the chain "
1097 "introduce a break, but that will be catastrophic if they "
1098 "should in fact be "
1099 "linked. Please check your structure, and add %s to "
1101 "if this was not correct.",
1102 *pdba->resinfo[i].name, pdba->resinfo[i].nr, restype->c_str(),
1103 *pdba->resinfo[*r_start].name, pdba->resinfo[*r_start].nr,
1104 startrestype->c_str(), *pdba->resinfo[i].name);
1106 if (endWarnings == 4)
1108 GMX_LOG(logger.warning)
1110 .appendTextFormatted(
1111 "Disabling further warnings about unidentified residues at end "
1121 GMX_LOG(logger.info)
1123 .appendTextFormatted("Identified residue %s%d as a ending terminus.",
1124 *pdba->resinfo[*r_end].name, pdba->resinfo[*r_end].nr);
1128 enum class ChainSeparationType : int
1137 static const gmx::EnumerationArray<ChainSeparationType, const char*> c_chainSeparationTypeNames = {
1138 { "id_or_ter", "id_and_ter", "ter", "id", "interactive" }
1140 static const gmx::EnumerationArray<ChainSeparationType, const char*> c_chainSeparationTypeNotificationMessages = {
1141 { "Splitting chemical chains based on TER records or chain id changing.\n",
1142 "Splitting chemical chains based on TER records and chain id changing.\n",
1143 "Splitting chemical chains based on TER records only (ignoring chain id).\n",
1144 "Splitting chemical chains based on changing chain id only (ignoring TER records).\n",
1145 "Splitting chemical chains interactively.\n" }
1148 static void modify_chain_numbers(t_atoms* pdba, ChainSeparationType chainSeparation, const gmx::MDLogger& logger)
1151 char old_prev_chainid;
1152 char old_this_chainid;
1153 int old_prev_chainnum;
1154 int old_this_chainnum;
1156 char select[STRLEN];
1160 const char* prev_atomname;
1161 const char* this_atomname;
1162 const char* prev_resname;
1163 const char* this_resname;
1169 /* The default chain enumeration is based on TER records only */
1170 printf("%s", c_chainSeparationTypeNotificationMessages[chainSeparation]);
1172 old_prev_chainid = '?';
1173 old_prev_chainnum = -1;
1176 this_atomname = nullptr;
1178 this_resname = nullptr;
1182 for (i = 0; i < pdba->nres; i++)
1184 ri = &pdba->resinfo[i];
1185 old_this_chainid = ri->chainid;
1186 old_this_chainnum = ri->chainnum;
1188 prev_atomname = this_atomname;
1189 prev_atomnum = this_atomnum;
1190 prev_resname = this_resname;
1191 prev_resnum = this_resnum;
1192 prev_chainid = this_chainid;
1194 this_atomname = *(pdba->atomname[i]);
1195 this_atomnum = (pdba->pdbinfo != nullptr) ? pdba->pdbinfo[i].atomnr : i + 1;
1196 this_resname = *ri->name;
1197 this_resnum = ri->nr;
1198 this_chainid = ri->chainid;
1200 switch (chainSeparation)
1202 case ChainSeparationType::IdOrTer:
1203 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1209 case ChainSeparationType::IdAndTer:
1210 if (old_this_chainid != old_prev_chainid && old_this_chainnum != old_prev_chainnum)
1216 case ChainSeparationType::Id:
1217 if (old_this_chainid != old_prev_chainid)
1223 case ChainSeparationType::Ter:
1224 if (old_this_chainnum != old_prev_chainnum)
1229 case ChainSeparationType::Interactive:
1230 if (old_this_chainid != old_prev_chainid || old_this_chainnum != old_prev_chainnum)
1234 GMX_LOG(logger.info)
1236 .appendTextFormatted(
1237 "Split the chain (and introduce termini) between residue %s%d (chain id '%c', atom %d %s)\
1239 "and residue %s%d (chain id '%c', atom %d %s) ? [n/y]",
1240 prev_resname, prev_resnum, prev_chainid, prev_atomnum,
1241 prev_atomname, this_resname, this_resnum, this_chainid,
1242 this_atomnum, this_atomname);
1244 if (nullptr == fgets(select, STRLEN - 1, stdin))
1246 gmx_fatal(FARGS, "Error reading from stdin");
1249 if (i == 0 || select[0] == 'y')
1255 case ChainSeparationType::Count:
1256 gmx_fatal(FARGS, "Internal inconsistency - this shouldn't happen...");
1258 old_prev_chainid = old_this_chainid;
1259 old_prev_chainnum = old_this_chainnum;
1261 ri->chainnum = new_chainnum;
1268 char chainnum = ' ';
1271 bool bAllWat = false;
1273 std::vector<int> chainstart;
1280 bool bAllWat = false;
1282 std::vector<int> chainstart;
1283 std::vector<MoleculePatchDatabase*> ntdb;
1284 std::vector<MoleculePatchDatabase*> ctdb;
1285 std::vector<int> r_start;
1286 std::vector<int> r_end;
1288 std::vector<gmx::RVec> x;
1291 enum class VSitesType : int
1298 static const gmx::EnumerationArray<VSitesType, const char*> c_vsitesTypeNames = {
1299 { "none", "hydrogens", "aromatics" }
1302 enum class WaterType : int
1314 static const gmx::EnumerationArray<WaterType, const char*> c_waterTypeNames = {
1315 { "select", "none", "spc", "spce", "tip3p", "tip4p", "tip5p", "tips3p" }
1318 enum class MergeType : int
1325 static const gmx::EnumerationArray<MergeType, const char*> c_mergeTypeNames = { { "no", "all",
1334 class pdb2gmx : public ICommandLineOptionsModule
1340 bVsiteAromatics_(FALSE),
1341 chainSeparation_(ChainSeparationType::IdOrTer),
1342 vsitesType_(VSitesType::None),
1343 waterType_(WaterType::Select),
1344 mergeType_(MergeType::No),
1350 // From ICommandLineOptionsModule
1351 void init(CommandLineModuleSettings* /*settings*/) override {}
1353 void initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings) override;
1355 void optionsFinished() override;
1373 bool bAllowMissing_;
1377 bool bChargeGroups_;
1387 bool bVsiteAromatics_;
1391 real long_bond_dist_;
1392 real short_bond_dist_;
1394 std::string indexOutputFile_;
1395 std::string outputFile_;
1396 std::string topologyFile_;
1397 std::string includeTopologyFile_;
1398 std::string outputConfFile_;
1399 std::string inputConfFile_;
1400 std::string outFile_;
1403 ChainSeparationType chainSeparation_;
1404 VSitesType vsitesType_;
1405 WaterType waterType_;
1406 MergeType mergeType_;
1409 char forcefield_[STRLEN];
1410 char ffdir_[STRLEN];
1413 std::vector<std::string> incls_;
1414 std::vector<t_mols> mols_;
1416 std::unique_ptr<gmx::MDLogger> loggerPointer_;
1419 void pdb2gmx::initOptions(IOptionsContainer* options, ICommandLineOptionsModuleSettings* settings)
1421 const char* desc[] = {
1422 "[THISMODULE] reads a [REF].pdb[ref] (or [REF].gro[ref]) file, reads",
1423 "some database files, adds hydrogens to the molecules and generates",
1424 "coordinates in GROMACS (GROMOS), or optionally [REF].pdb[ref], format",
1425 "and a topology in GROMACS format.",
1426 "These files can subsequently be processed to generate a run input file.",
1428 "[THISMODULE] will search for force fields by looking for",
1429 "a [TT]forcefield.itp[tt] file in subdirectories [TT]<forcefield>.ff[tt]",
1430 "of the current working directory and of the GROMACS library directory",
1431 "as inferred from the path of the binary or the [TT]GMXLIB[tt] environment",
1433 "By default the forcefield selection is interactive,",
1434 "but you can use the [TT]-ff[tt] option to specify one of the short names",
1435 "in the list on the command line instead. In that case [THISMODULE] just looks",
1436 "for the corresponding [TT]<forcefield>.ff[tt] directory.",
1438 "After choosing a force field, all files will be read only from",
1439 "the corresponding force field directory.",
1440 "If you want to modify or add a residue types, you can copy the force",
1441 "field directory from the GROMACS library directory to your current",
1442 "working directory. If you want to add new protein residue types,",
1443 "you will need to modify [TT]residuetypes.dat[tt] in the library directory",
1444 "or copy the whole library directory to a local directory and set",
1445 "the environment variable [TT]GMXLIB[tt] to the name of that directory.",
1446 "Check Chapter 5 of the manual for more information about file formats.",
1449 "Note that a [REF].pdb[ref] file is nothing more than a file format, and it",
1450 "need not necessarily contain a protein structure. Every kind of",
1451 "molecule for which there is support in the database can be converted.",
1452 "If there is no support in the database, you can add it yourself.[PAR]",
1454 "The program has limited intelligence, it reads a number of database",
1455 "files, that allow it to make special bonds (Cys-Cys, Heme-His, etc.),",
1456 "if necessary this can be done manually. The program can prompt the",
1457 "user to select which kind of LYS, ASP, GLU, CYS or HIS residue is",
1458 "desired. For Lys the choice is between neutral (two protons on NZ) or",
1459 "protonated (three protons, default), for Asp and Glu unprotonated",
1460 "(default) or protonated, for His the proton can be either on ND1,",
1461 "on NE2 or on both. By default these selections are done automatically.",
1462 "For His, this is based on an optimal hydrogen bonding",
1463 "conformation. Hydrogen bonds are defined based on a simple geometric",
1464 "criterion, specified by the maximum hydrogen-donor-acceptor angle",
1465 "and donor-acceptor distance, which are set by [TT]-angle[tt] and",
1466 "[TT]-dist[tt] respectively.[PAR]",
1468 "The protonation state of N- and C-termini can be chosen interactively",
1469 "with the [TT]-ter[tt] flag. Default termini are ionized (NH3+ and COO-),",
1470 "respectively. Some force fields support zwitterionic forms for chains of",
1471 "one residue, but for polypeptides these options should NOT be selected.",
1472 "The AMBER force fields have unique forms for the terminal residues,",
1473 "and these are incompatible with the [TT]-ter[tt] mechanism. You need",
1474 "to prefix your N- or C-terminal residue names with \"N\" or \"C\"",
1475 "respectively to use these forms, making sure you preserve the format",
1476 "of the coordinate file. Alternatively, use named terminating residues",
1477 "(e.g. ACE, NME).[PAR]",
1479 "The separation of chains is not entirely trivial since the markup",
1480 "in user-generated PDB files frequently varies and sometimes it",
1481 "is desirable to merge entries across a TER record, for instance",
1482 "if you want a disulfide bridge or distance restraints between",
1483 "two protein chains or if you have a HEME group bound to a protein.",
1484 "In such cases multiple chains should be contained in a single",
1485 "[TT]moleculetype[tt] definition.",
1486 "To handle this, [THISMODULE] uses two separate options.",
1487 "First, [TT]-chainsep[tt] allows you to choose when a new chemical chain should",
1488 "start, and termini added when applicable. This can be done based on the",
1489 "existence of TER records, when the chain id changes, or combinations of either",
1490 "or both of these. You can also do the selection fully interactively.",
1491 "In addition, there is a [TT]-merge[tt] option that controls how multiple chains",
1492 "are merged into one moleculetype, after adding all the chemical termini (or not).",
1493 "This can be turned off (no merging), all non-water chains can be merged into a",
1494 "single molecule, or the selection can be done interactively.[PAR]",
1496 "[THISMODULE] will also check the occupancy field of the [REF].pdb[ref] file.",
1497 "If any of the occupancies are not one, indicating that the atom is",
1498 "not resolved well in the structure, a warning message is issued.",
1499 "When a [REF].pdb[ref] file does not originate from an X-ray structure determination",
1500 "all occupancy fields may be zero. Either way, it is up to the user",
1501 "to verify the correctness of the input data (read the article!).[PAR]",
1503 "During processing the atoms will be reordered according to GROMACS",
1504 "conventions. With [TT]-n[tt] an index file can be generated that",
1505 "contains one group reordered in the same way. This allows you to",
1506 "convert a GROMOS trajectory and coordinate file to GROMOS. There is",
1507 "one limitation: reordering is done after the hydrogens are stripped",
1508 "from the input and before new hydrogens are added. This means that",
1509 "you should not use [TT]-ignh[tt].[PAR]",
1511 "The [REF].gro[ref] and [TT].g96[tt] file formats do not support chain",
1512 "identifiers. Therefore it is useful to enter a [REF].pdb[ref] file name at",
1513 "the [TT]-o[tt] option when you want to convert a multi-chain [REF].pdb[ref] file.",
1516 "The option [TT]-vsite[tt] removes hydrogen and fast improper dihedral",
1517 "motions. Angular and out-of-plane motions can be removed by changing",
1518 "hydrogens into virtual sites and fixing angles, which fixes their",
1519 "position relative to neighboring atoms. Additionally, all atoms in the",
1520 "aromatic rings of the standard amino acids (i.e. PHE, TRP, TYR and HIS)",
1521 "can be converted into virtual sites, eliminating the fast improper dihedral",
1522 "fluctuations in these rings (but this feature is deprecated).",
1523 "[BB]Note[bb] that in this case all other hydrogen",
1524 "atoms are also converted to virtual sites. The mass of all atoms that are",
1525 "converted into virtual sites, is added to the heavy atoms.[PAR]",
1526 "Also slowing down of dihedral motion can be done with [TT]-heavyh[tt]",
1527 "done by increasing the hydrogen-mass by a factor of 4. This is also",
1528 "done for water hydrogens to slow down the rotational motion of water.",
1529 "The increase in mass of the hydrogens is subtracted from the bonded",
1530 "(heavy) atom so that the total mass of the system remains the same."
1533 settings->setHelpText(desc);
1535 options->addOption(BooleanOption("newrtp").store(&bNewRTP_).defaultValue(false).hidden().description(
1536 "Write the residue database in new format to [TT]new.rtp[tt]"));
1538 RealOption("lb").store(&long_bond_dist_).defaultValue(0.25).hidden().description("Long bond warning distance"));
1540 RealOption("sb").store(&short_bond_dist_).defaultValue(0.05).hidden().description("Short bond warning distance"));
1541 options->addOption(EnumOption<ChainSeparationType>("chainsep")
1542 .enumValue(c_chainSeparationTypeNames)
1543 .store(&chainSeparation_)
1544 .description("Condition in PDB files when a new chain should be "
1545 "started (adding termini)"));
1546 options->addOption(EnumOption<MergeType>("merge")
1547 .enumValue(c_mergeTypeNames)
1549 .description("Merge multiple chains into a single [moleculetype]"));
1550 options->addOption(StringOption("ff").store(&ff_).defaultValue("select").description(
1551 "Force field, interactive by default. Use [TT]-h[tt] for information."));
1553 EnumOption<WaterType>("water").store(&waterType_).enumValue(c_waterTypeNames).description("Water model to use"));
1554 options->addOption(BooleanOption("inter").store(&bInter_).defaultValue(false).description(
1555 "Set the next 8 options to interactive"));
1556 options->addOption(BooleanOption("ss").store(&bCysMan_).defaultValue(false).description(
1557 "Interactive SS bridge selection"));
1558 options->addOption(BooleanOption("ter").store(&bTerMan_).defaultValue(false).description(
1559 "Interactive termini selection, instead of charged (default)"));
1560 options->addOption(BooleanOption("lys").store(&bLysMan_).defaultValue(false).description(
1561 "Interactive lysine selection, instead of charged"));
1562 options->addOption(BooleanOption("arg").store(&bArgMan_).defaultValue(false).description(
1563 "Interactive arginine selection, instead of charged"));
1564 options->addOption(BooleanOption("asp").store(&bAspMan_).defaultValue(false).description(
1565 "Interactive aspartic acid selection, instead of charged"));
1566 options->addOption(BooleanOption("glu").store(&bGluMan_).defaultValue(false).description(
1567 "Interactive glutamic acid selection, instead of charged"));
1568 options->addOption(BooleanOption("gln").store(&bGlnMan_).defaultValue(false).description(
1569 "Interactive glutamine selection, instead of charged"));
1570 options->addOption(BooleanOption("his").store(&bHisMan_).defaultValue(false).description(
1571 "Interactive histidine selection, instead of checking H-bonds"));
1572 options->addOption(RealOption("angle").store(&angle_).defaultValue(135.0).description(
1573 "Minimum hydrogen-donor-acceptor angle for a H-bond (degrees)"));
1575 RealOption("dist").store(&distance_).defaultValue(0.3).description("Maximum donor-acceptor distance for a H-bond (nm)"));
1576 options->addOption(BooleanOption("una").store(&bUnA_).defaultValue(false).description(
1577 "Select aromatic rings with united CH atoms on phenylalanine, tryptophane and "
1579 options->addOption(BooleanOption("sort").store(&bSort_).defaultValue(true).hidden().description(
1580 "Sort the residues according to database, turning this off is dangerous as charge "
1581 "groups might be broken in parts"));
1583 BooleanOption("ignh").store(&bRemoveH_).defaultValue(false).description("Ignore hydrogen atoms that are in the coordinate file"));
1585 BooleanOption("missing")
1586 .store(&bAllowMissing_)
1587 .defaultValue(false)
1589 "Continue when atoms are missing and bonds cannot be made, dangerous"));
1591 BooleanOption("v").store(&bVerbose_).defaultValue(false).description("Be slightly more verbose in messages"));
1593 RealOption("posrefc").store(&posre_fc_).defaultValue(1000).description("Force constant for position restraints"));
1594 options->addOption(EnumOption<VSitesType>("vsite")
1595 .store(&vsitesType_)
1596 .enumValue(c_vsitesTypeNames)
1597 .description("Convert atoms to virtual sites"));
1598 options->addOption(BooleanOption("heavyh").store(&bHeavyH_).defaultValue(false).description(
1599 "Make hydrogen atoms heavy"));
1601 BooleanOption("deuterate").store(&bDeuterate_).defaultValue(false).description("Change the mass of hydrogens to 2 amu"));
1602 options->addOption(BooleanOption("chargegrp")
1603 .store(&bChargeGroups_)
1605 .description("Use charge groups in the [REF].rtp[ref] file"));
1606 options->addOption(BooleanOption("cmap").store(&bCmap_).defaultValue(true).description(
1607 "Use cmap torsions (if enabled in the [REF].rtp[ref] file)"));
1608 options->addOption(BooleanOption("renum")
1610 .defaultValue(false)
1611 .description("Renumber the residues consecutively in the output"));
1612 options->addOption(BooleanOption("rtpres")
1613 .store(&bRTPresname_)
1614 .defaultValue(false)
1615 .description("Use [REF].rtp[ref] entry names as residue names"));
1616 options->addOption(FileNameOption("f")
1619 .store(&inputConfFile_)
1621 .defaultBasename("protein")
1623 .description("Structure file"));
1624 options->addOption(FileNameOption("o")
1627 .store(&outputConfFile_)
1629 .defaultBasename("conf")
1630 .description("Structure file"));
1631 options->addOption(FileNameOption("p")
1634 .store(&topologyFile_)
1636 .defaultBasename("topol")
1637 .description("Topology file"));
1638 options->addOption(FileNameOption("i")
1641 .store(&includeTopologyFile_)
1643 .defaultBasename("posre")
1644 .description("Include file for topology"));
1645 options->addOption(FileNameOption("n")
1648 .store(&indexOutputFile_)
1649 .storeIsSet(&bIndexSet_)
1650 .defaultBasename("index")
1651 .description("Index file"));
1652 options->addOption(FileNameOption("q")
1656 .storeIsSet(&bOutputSet_)
1657 .defaultBasename("clean")
1659 .description("Structure file"));
1662 void pdb2gmx::optionsFinished()
1664 if (inputConfFile_.empty())
1666 GMX_THROW(InconsistentInputError("You must supply an input file"));
1670 /* if anything changes here, also change description of -inter */
1685 else if (bDeuterate_)
1694 /* Force field selection, interactive or direct */
1695 choose_ff(strcmp(ff_.c_str(), "select") == 0 ? nullptr : ff_.c_str(), forcefield_,
1696 sizeof(forcefield_), ffdir_, sizeof(ffdir_), *loggerPointer_);
1698 if (strlen(forcefield_) > 0)
1700 ffname_ = forcefield_;
1701 ffname_[0] = std::toupper(ffname_[0]);
1705 gmx_fatal(FARGS, "Empty forcefield string");
1711 char select[STRLEN];
1712 std::vector<DisulfideBond> ssbonds;
1716 const char* prev_atomname;
1717 const char* this_atomname;
1718 const char* prev_resname;
1719 const char* this_resname;
1724 int prev_chainnumber;
1725 int this_chainnumber;
1726 int this_chainstart;
1727 int prev_chainstart;
1729 gmx::LoggerBuilder builder;
1730 builder.addTargetStream(gmx::MDLogger::LogLevel::Info, &gmx::TextOutputFile::standardOutput());
1731 builder.addTargetStream(gmx::MDLogger::LogLevel::Warning, &gmx::TextOutputFile::standardError());
1732 gmx::LoggerOwner logOwner(builder.build());
1733 loggerPointer_ = std::make_unique<gmx::MDLogger>(logOwner.logger());
1734 const gmx::MDLogger& logger = *loggerPointer_;
1736 GMX_LOG(logger.info)
1738 .appendTextFormatted("Using the %s force field in directory %s", ffname_, ffdir_);
1740 choose_watermodel(c_waterTypeNames[waterType_], ffdir_, &watermodel_, logger);
1742 switch (vsitesType_)
1744 case VSitesType::None:
1746 bVsiteAromatics_ = false;
1748 case VSitesType::Hydrogens:
1750 bVsiteAromatics_ = false;
1752 case VSitesType::Aromatics:
1754 bVsiteAromatics_ = true;
1756 case VSitesType::Count:
1757 gmx_fatal(FARGS, "Internal inconsistency: VSitesType='%s'", c_vsitesTypeNames[vsitesType_]);
1760 /* Open the symbol table */
1762 open_symtab(&symtab);
1764 /* Residue type database */
1767 /* Read residue renaming database(s), if present */
1768 std::vector<std::string> rrn = fflib_search_file_end(ffdir_, ".r2b", FALSE);
1770 std::vector<RtpRename> rtprename;
1771 for (const auto& filename : rrn)
1773 GMX_LOG(logger.info).asParagraph().appendTextFormatted("going to rename %s", filename.c_str());
1774 FILE* fp = fflib_open(filename);
1775 read_rtprename(filename.c_str(), fp, &rtprename);
1779 /* Add all alternative names from the residue renaming database to the list
1780 of recognized amino/nucleic acids. */
1781 for (const auto& rename : rtprename)
1783 /* Only add names if the 'standard' gromacs/iupac base name was found */
1784 if (auto restype = rt.optionalTypeOfNamedDatabaseResidue(rename.gmx))
1786 rt.addResidue(rename.main, *restype);
1787 rt.addResidue(rename.nter, *restype);
1788 rt.addResidue(rename.cter, *restype);
1789 rt.addResidue(rename.bter, *restype);
1796 if (watermodel_ != nullptr && (strstr(watermodel_, "4p") || strstr(watermodel_, "4P")))
1800 else if (watermodel_ != nullptr && (strstr(watermodel_, "5p") || strstr(watermodel_, "5P")))
1810 char* title = nullptr;
1814 int natom = read_pdball(inputConfFile_.c_str(), bOutputSet_, outFile_.c_str(), &title, &pdba_all,
1815 &pdbx, &pbcType, box, bRemoveH_, &symtab, &rt, watres, &aps, bVerbose_);
1819 std::string message = formatString("No atoms found in pdb file %s\n", inputConfFile_.c_str());
1820 GMX_THROW(InconsistentInputError(message));
1823 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Analyzing pdb file");
1824 int nwaterchain = 0;
1826 modify_chain_numbers(&pdba_all, chainSeparation_, logger);
1828 int nchainmerges = 0;
1830 this_atomname = nullptr;
1832 this_resname = nullptr;
1835 this_chainnumber = -1;
1836 this_chainstart = 0;
1837 /* Keep the compiler happy */
1838 prev_chainstart = 0;
1841 std::vector<t_pdbchain> pdb_ch;
1844 bool bMerged = false;
1845 for (int i = 0; (i < natom); i++)
1847 ri = &pdba_all.resinfo[pdba_all.atom[i].resind];
1849 /* TODO this should live in a helper object, and consolidate
1850 that with code in modify_chain_numbers */
1851 prev_atomname = this_atomname;
1852 prev_atomnum = this_atomnum;
1853 prev_resname = this_resname;
1854 prev_resnum = this_resnum;
1855 prev_chainid = this_chainid;
1856 prev_chainnumber = this_chainnumber;
1859 prev_chainstart = this_chainstart;
1862 this_atomname = *pdba_all.atomname[i];
1863 this_atomnum = (pdba_all.pdbinfo != nullptr) ? pdba_all.pdbinfo[i].atomnr : i + 1;
1864 this_resname = *ri->name;
1865 this_resnum = ri->nr;
1866 this_chainid = ri->chainid;
1867 this_chainnumber = ri->chainnum;
1869 bWat_ = gmx::equalCaseInsensitive(*ri->name, watres);
1871 if ((i == 0) || (this_chainnumber != prev_chainnumber) || (bWat_ != bPrevWat_))
1875 "Must have pdbinfo from reading a PDB file if chain number is changing");
1876 this_chainstart = pdba_all.atom[i].resind;
1878 if (i > 0 && !bWat_)
1880 if (!strncmp(c_mergeTypeNames[mergeType_], "int", 3))
1882 GMX_LOG(logger.info)
1884 .appendTextFormatted(
1885 "Merge chain ending with residue %s%d (chain id '%c', atom %d "
1886 "%s) and chain starting with "
1887 "residue %s%d (chain id '%c', atom %d %s) into a single "
1888 "moleculetype (keeping termini)? [n/y]",
1889 prev_resname, prev_resnum, prev_chainid, prev_atomnum,
1890 prev_atomname, this_resname, this_resnum, this_chainid,
1891 this_atomnum, this_atomname);
1893 if (nullptr == fgets(select, STRLEN - 1, stdin))
1895 gmx_fatal(FARGS, "Error reading from stdin");
1897 bMerged = (select[0] == 'y');
1899 else if (!strncmp(c_mergeTypeNames[mergeType_], "all", 3))
1907 pdb_ch[numChains - 1].chainstart[pdb_ch[numChains - 1].nterpairs] =
1908 pdba_all.atom[i].resind - prev_chainstart;
1909 pdb_ch[numChains - 1].nterpairs++;
1910 pdb_ch[numChains - 1].chainstart.resize(pdb_ch[numChains - 1].nterpairs + 1);
1915 /* set natom for previous chain */
1918 pdb_ch[numChains - 1].natom = i - pdb_ch[numChains - 1].start;
1925 /* check if chain identifier was used before */
1926 for (int j = 0; (j < numChains); j++)
1928 if (pdb_ch[j].chainid != ' ' && pdb_ch[j].chainid == ri->chainid)
1930 GMX_LOG(logger.warning)
1932 .appendTextFormatted(
1933 "Chain identifier '%c' is used in two non-sequential "
1935 "They will be treated as separate chains unless you "
1936 "reorder your file.",
1940 t_pdbchain newChain;
1941 newChain.chainid = ri->chainid;
1942 newChain.chainnum = ri->chainnum;
1944 newChain.bAllWat = bWat_;
1947 newChain.nterpairs = 0;
1951 newChain.nterpairs = 1;
1953 newChain.chainstart.resize(newChain.nterpairs + 1);
1954 /* modified [numChains] to [0] below */
1955 newChain.chainstart[0] = 0;
1956 pdb_ch.push_back(newChain);
1962 pdb_ch.back().natom = natom - pdb_ch.back().start;
1964 /* set all the water blocks at the end of the chain */
1965 std::vector<int> swap_index(numChains);
1967 for (int i = 0; i < numChains; i++)
1969 if (!pdb_ch[i].bAllWat)
1975 for (int i = 0; i < numChains; i++)
1977 if (pdb_ch[i].bAllWat)
1983 if (nwaterchain > 1)
1985 GMX_LOG(logger.info)
1987 .appendTextFormatted("Moved all the water blocks to the end");
1991 std::vector<t_chain> chains(numChains);
1992 /* copy pdb data and x for all chains */
1993 for (int i = 0; (i < numChains); i++)
1995 int si = swap_index[i];
1996 chains[i].chainid = pdb_ch[si].chainid;
1997 chains[i].chainnum = pdb_ch[si].chainnum;
1998 chains[i].bAllWat = pdb_ch[si].bAllWat;
1999 chains[i].nterpairs = pdb_ch[si].nterpairs;
2000 chains[i].chainstart = pdb_ch[si].chainstart;
2001 chains[i].ntdb.clear();
2002 chains[i].ctdb.clear();
2003 chains[i].r_start.resize(pdb_ch[si].nterpairs);
2004 chains[i].r_end.resize(pdb_ch[si].nterpairs);
2006 snew(chains[i].pdba, 1);
2007 init_t_atoms(chains[i].pdba, pdb_ch[si].natom, true);
2008 for (j = 0; j < chains[i].pdba->nr; j++)
2010 chains[i].pdba->atom[j] = pdba_all.atom[pdb_ch[si].start + j];
2011 chains[i].pdba->atomname[j] = put_symtab(&symtab, *pdba_all.atomname[pdb_ch[si].start + j]);
2012 chains[i].pdba->pdbinfo[j] = pdba_all.pdbinfo[pdb_ch[si].start + j];
2013 chains[i].x.emplace_back(pdbx[pdb_ch[si].start + j]);
2015 /* Re-index the residues assuming that the indices are continuous */
2016 int k = chains[i].pdba->atom[0].resind;
2017 int nres = chains[i].pdba->atom[chains[i].pdba->nr - 1].resind - k + 1;
2018 chains[i].pdba->nres = nres;
2019 for (int j = 0; j < chains[i].pdba->nr; j++)
2021 chains[i].pdba->atom[j].resind -= k;
2023 srenew(chains[i].pdba->resinfo, nres);
2024 for (int j = 0; j < nres; j++)
2026 chains[i].pdba->resinfo[j] = pdba_all.resinfo[k + j];
2027 chains[i].pdba->resinfo[j].name = put_symtab(&symtab, *pdba_all.resinfo[k + j].name);
2028 /* make all chain identifiers equal to that of the chain */
2029 chains[i].pdba->resinfo[j].chainid = pdb_ch[si].chainid;
2033 if (nchainmerges > 0)
2035 GMX_LOG(logger.info)
2037 .appendTextFormatted("Merged chains into joint molecule definitions at %d places.",
2041 GMX_LOG(logger.info)
2043 .appendTextFormatted(
2044 "There are %d chains and %d blocks of water and "
2045 "%d residues with %d atoms",
2046 numChains - nwaterchain, nwaterchain, pdba_all.nres, natom);
2048 GMX_LOG(logger.info)
2050 .appendTextFormatted(" %5s %4s %6s", "chain", "#res", "#atoms");
2051 for (int i = 0; (i < numChains); i++)
2053 GMX_LOG(logger.info)
2055 .appendTextFormatted(" %d '%c' %5d %6d %s\n", i + 1,
2056 chains[i].chainid ? chains[i].chainid : '-', chains[i].pdba->nres,
2057 chains[i].pdba->nr, chains[i].bAllWat ? "(only water)" : "");
2060 check_occupancy(&pdba_all, inputConfFile_.c_str(), bVerbose_, logger);
2062 /* Read atomtypes... */
2063 PreprocessingAtomTypes atype = read_atype(ffdir_, &symtab);
2065 /* read residue database */
2066 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Reading residue database... (%s)", forcefield_);
2067 std::vector<std::string> rtpf = fflib_search_file_end(ffdir_, ".rtp", true);
2068 std::vector<PreprocessResidue> rtpFFDB;
2069 for (const auto& filename : rtpf)
2071 readResidueDatabase(filename, &rtpFFDB, &atype, &symtab, logger, false);
2075 /* Not correct with multiple rtp input files with different bonded types */
2076 FILE* fp = gmx_fio_fopen("new.rtp", "w");
2077 print_resall(fp, rtpFFDB, atype);
2081 /* read hydrogen database */
2082 std::vector<MoleculePatchDatabase> ah;
2083 read_h_db(ffdir_, &ah);
2085 /* Read Termini database... */
2086 std::vector<MoleculePatchDatabase> ntdb;
2087 std::vector<MoleculePatchDatabase> ctdb;
2088 std::vector<MoleculePatchDatabase*> tdblist;
2089 int nNtdb = read_ter_db(ffdir_, 'n', &ntdb, &atype);
2090 int nCtdb = read_ter_db(ffdir_, 'c', &ctdb, &atype);
2092 FILE* top_file = gmx_fio_fopen(topologyFile_.c_str(), "w");
2094 print_top_header(top_file, topologyFile_.c_str(), FALSE, ffdir_, mHmult_);
2097 std::vector<gmx::RVec> x;
2098 /* new pdb datastructure for sorting. */
2099 t_atoms** sortAtoms = nullptr;
2100 t_atoms** localAtoms = nullptr;
2101 snew(sortAtoms, numChains);
2102 snew(localAtoms, numChains);
2103 for (int chain = 0; (chain < numChains); chain++)
2105 cc = &(chains[chain]);
2107 /* set pdba, natom and nres to the current chain */
2110 natom = cc->pdba->nr;
2111 int nres = cc->pdba->nres;
2113 if (cc->chainid && (cc->chainid != ' '))
2115 GMX_LOG(logger.info)
2117 .appendTextFormatted("Processing chain %d '%c' (%d atoms, %d residues)",
2118 chain + 1, cc->chainid, natom, nres);
2122 GMX_LOG(logger.info)
2124 .appendTextFormatted("Processing chain %d (%d atoms, %d residues)", chain + 1,
2128 process_chain(pdba, x, bUnA_, bUnA_, bUnA_, bLysMan_, bAspMan_, bGluMan_, bHisMan_,
2129 bArgMan_, bGlnMan_, angle_, distance_, &symtab, rtprename);
2131 cc->chainstart[cc->nterpairs] = pdba->nres;
2133 for (int i = 0; i < cc->nterpairs; i++)
2135 find_nc_ter(pdba, cc->chainstart[i], cc->chainstart[i + 1], &(cc->r_start[j]),
2136 &(cc->r_end[j]), &rt, logger);
2138 if (cc->r_start[j] >= 0 && cc->r_end[j] >= 0)
2144 if (cc->nterpairs == 0)
2146 GMX_LOG(logger.info)
2148 .appendTextFormatted(
2149 "Problem with chain definition, or missing terminal residues. "
2150 "This chain does not appear to contain a recognized chain molecule. "
2151 "If this is incorrect, you can edit residuetypes.dat to modify the "
2155 /* Check for disulfides and other special bonds */
2156 ssbonds = makeDisulfideBonds(pdba, gmx::as_rvec_array(x.data()), bCysMan_, bVerbose_);
2158 if (!rtprename.empty())
2160 rename_resrtp(pdba, cc->nterpairs, cc->r_start, cc->r_end, rtprename, &symtab, bVerbose_, logger);
2163 for (int i = 0; i < cc->nterpairs; i++)
2166 * We first apply a filter so we only have the
2167 * termini that can be applied to the residue in question
2168 * (or a generic terminus if no-residue specific is available).
2170 /* First the N terminus */
2173 tdblist = filter_ter(ntdb, *pdba->resinfo[cc->r_start[i]].name);
2174 if (tdblist.empty())
2176 GMX_LOG(logger.info)
2178 .appendTextFormatted(
2179 "No suitable end (N or 5') terminus found in database - "
2180 "assuming this residue "
2181 "is already in a terminus-specific form and skipping terminus "
2183 cc->ntdb.push_back(nullptr);
2187 if (bTerMan_ && !tdblist.empty())
2189 sprintf(select, "Select start terminus type for %s-%d",
2190 *pdba->resinfo[cc->r_start[i]].name, pdba->resinfo[cc->r_start[i]].nr);
2191 cc->ntdb.push_back(choose_ter(tdblist, select));
2195 cc->ntdb.push_back(tdblist[0]);
2198 printf("Start terminus %s-%d: %s\n", *pdba->resinfo[cc->r_start[i]].name,
2199 pdba->resinfo[cc->r_start[i]].nr, (cc->ntdb[i])->name.c_str());
2205 cc->ntdb.push_back(nullptr);
2208 /* And the C terminus */
2211 tdblist = filter_ter(ctdb, *pdba->resinfo[cc->r_end[i]].name);
2212 if (tdblist.empty())
2214 GMX_LOG(logger.info)
2216 .appendTextFormatted(
2217 "No suitable end (C or 3') terminus found in database - "
2218 "assuming this residue"
2219 "is already in a terminus-specific form and skipping terminus "
2221 cc->ctdb.push_back(nullptr);
2225 if (bTerMan_ && !tdblist.empty())
2227 sprintf(select, "Select end terminus type for %s-%d",
2228 *pdba->resinfo[cc->r_end[i]].name, pdba->resinfo[cc->r_end[i]].nr);
2229 cc->ctdb.push_back(choose_ter(tdblist, select));
2233 cc->ctdb.push_back(tdblist[0]);
2235 printf("End terminus %s-%d: %s\n", *pdba->resinfo[cc->r_end[i]].name,
2236 pdba->resinfo[cc->r_end[i]].nr, (cc->ctdb[i])->name.c_str());
2242 cc->ctdb.push_back(nullptr);
2245 std::vector<MoleculePatchDatabase> hb_chain;
2246 /* lookup hackblocks and rtp for all residues */
2247 std::vector<PreprocessResidue> restp_chain;
2248 get_hackblocks_rtp(&hb_chain, &restp_chain, rtpFFDB, pdba->nres, pdba->resinfo, cc->nterpairs,
2249 &symtab, cc->ntdb, cc->ctdb, cc->r_start, cc->r_end, bAllowMissing_, logger);
2250 /* ideally, now we would not need the rtp itself anymore, but do
2251 everything using the hb and restp arrays. Unfortunately, that
2252 requires some re-thinking of code in gen_vsite.c, which I won't
2253 do now :( AF 26-7-99 */
2255 rename_atoms(nullptr, ffdir_, pdba, &symtab, restp_chain, false, &rt, false, bVerbose_);
2257 match_atomnames_with_rtp(restp_chain, hb_chain, pdba, &symtab, x, bVerbose_, logger);
2262 t_blocka* block = new_blocka();
2264 sort_pdbatoms(restp_chain, natom, &pdba, &sortAtoms[chain], &x, block, &gnames);
2265 remove_duplicate_atoms(pdba, x, bVerbose_, logger);
2270 GMX_LOG(logger.warning)
2272 .appendTextFormatted(
2273 "With the -remh option the generated "
2274 "index file (%s) might be useless "
2275 "(the index file is generated before hydrogens are added)",
2276 indexOutputFile_.c_str());
2278 write_index(indexOutputFile_.c_str(), block, gnames, false, 0);
2280 for (int i = 0; i < block->nr; i++)
2290 GMX_LOG(logger.warning)
2292 .appendTextFormatted(
2293 "Without sorting no check for duplicate atoms can be done");
2296 /* Generate Hydrogen atoms (and termini) in the sequence */
2297 GMX_LOG(logger.info)
2299 .appendTextFormatted(
2300 "Generating any missing hydrogen atoms and/or adding termini.");
2301 add_h(&pdba, &localAtoms[chain], &x, ah, &symtab, cc->nterpairs, cc->ntdb, cc->ctdb,
2302 cc->r_start, cc->r_end, bAllowMissing_);
2303 GMX_LOG(logger.info)
2305 .appendTextFormatted("Now there are %d residues with %d atoms", pdba->nres, pdba->nr);
2307 /* make up molecule name(s) */
2309 int k = (cc->nterpairs > 0 && cc->r_start[0] >= 0) ? cc->r_start[0] : 0;
2311 std::string restype = rt.typeOfNamedDatabaseResidue(*pdba->resinfo[k].name);
2313 std::string molname;
2321 this_chainid = cc->chainid;
2323 /* Add the chain id if we have one */
2324 if (this_chainid != ' ')
2326 suffix.append(formatString("_chain_%c", this_chainid));
2329 /* Check if there have been previous chains with the same id */
2331 for (int k = 0; k < chain; k++)
2333 if (cc->chainid == chains[k].chainid)
2338 /* Add the number for this chain identifier if there are multiple copies */
2341 suffix.append(formatString("%d", nid_used + 1));
2344 if (suffix.length() > 0)
2346 molname.append(restype);
2347 molname.append(suffix);
2354 std::string itp_fn = topologyFile_;
2356 std::string posre_fn = includeTopologyFile_;
2357 if ((numChains - nwaterchain > 1) && !cc->bAllWat)
2360 printf("Chain time...\n");
2361 // construct the itp file name
2362 itp_fn = stripSuffixIfPresent(itp_fn, ".top");
2364 itp_fn.append(molname);
2365 itp_fn.append(".itp");
2366 // now do the same for posre
2367 posre_fn = stripSuffixIfPresent(posre_fn, ".itp");
2368 posre_fn.append("_");
2369 posre_fn.append(molname);
2370 posre_fn.append(".itp");
2371 if (posre_fn == itp_fn)
2373 posre_fn = Path::concatenateBeforeExtension(posre_fn, "_pr");
2375 incls_.emplace_back();
2376 incls_.back() = itp_fn;
2377 itp_file_ = gmx_fio_fopen(itp_fn.c_str(), "w");
2384 mols_.emplace_back();
2387 mols_.back().name = "SOL";
2388 mols_.back().nr = pdba->nres;
2392 mols_.back().name = molname;
2393 mols_.back().nr = 1;
2398 print_top_comment(itp_file_, itp_fn.c_str(), ffdir_, true);
2404 top_file2 = nullptr;
2408 top_file2 = itp_file_;
2412 top_file2 = top_file;
2415 pdb2top(top_file2, posre_fn.c_str(), molname.c_str(), pdba, &x, &atype, &symtab, rtpFFDB,
2416 restp_chain, hb_chain, bAllowMissing_, bVsites_, bVsiteAromatics_, ffdir_, mHmult_,
2417 ssbonds, long_bond_dist_, short_bond_dist_, bDeuterate_, bChargeGroups_, bCmap_,
2418 bRenumRes_, bRTPresname_, logger);
2422 write_posres(posre_fn.c_str(), pdba, posre_fc_);
2427 gmx_fio_fclose(itp_file_);
2430 /* pdba and natom have been reassigned somewhere so: */
2435 if (watermodel_ == nullptr)
2437 for (int chain = 0; chain < numChains; chain++)
2439 if (chains[chain].bAllWat)
2441 auto message = formatString(
2442 "You have chosen not to include a water model, "
2443 "but there is water in the input file. Select a "
2444 "water model or remove the water from your input file.");
2445 GMX_THROW(InconsistentInputError(message));
2451 std::string waterFile = formatString("%s%c%s.itp", ffdir_, DIR_SEPARATOR, watermodel_);
2452 if (!fflib_fexist(waterFile))
2454 auto message = formatString(
2455 "The topology file '%s' for the selected water "
2456 "model '%s' can not be found in the force field "
2457 "directory. Select a different water model.",
2458 waterFile.c_str(), watermodel_);
2459 GMX_THROW(InconsistentInputError(message));
2463 print_top_mols(top_file, title, ffdir_, watermodel_, incls_, mols_);
2464 gmx_fio_fclose(top_file);
2466 /* now merge all chains back together */
2469 for (int i = 0; (i < numChains); i++)
2471 natom += chains[i].pdba->nr;
2472 nres += chains[i].pdba->nres;
2476 init_t_atoms(atoms, natom, false);
2477 for (int i = 0; i < atoms->nres; i++)
2479 sfree(atoms->resinfo[i].name);
2482 srenew(atoms->resinfo, nres);
2486 for (int i = 0; (i < numChains); i++)
2490 GMX_LOG(logger.info)
2492 .appendTextFormatted("Including chain %d in system: %d atoms %d residues",
2493 i + 1, chains[i].pdba->nr, chains[i].pdba->nres);
2495 for (int j = 0; (j < chains[i].pdba->nr); j++)
2497 atoms->atom[k] = chains[i].pdba->atom[j];
2498 atoms->atom[k].resind += l; /* l is processed nr of residues */
2499 atoms->atomname[k] = chains[i].pdba->atomname[j];
2500 atoms->resinfo[atoms->atom[k].resind].chainid = chains[i].chainid;
2501 x.push_back(chains[i].x[j]);
2504 for (int j = 0; (j < chains[i].pdba->nres); j++)
2506 atoms->resinfo[l] = chains[i].pdba->resinfo[j];
2509 atoms->resinfo[l].name = atoms->resinfo[l].rtp;
2513 done_atom(chains[i].pdba);
2518 GMX_LOG(logger.info)
2520 .appendTextFormatted("Now there are %d atoms and %d residues", k, l);
2521 print_sums(atoms, true, logger);
2525 GMX_LOG(logger.info).asParagraph().appendTextFormatted("Writing coordinate file...");
2526 clear_rvec(box_space);
2529 make_new_box(atoms->nr, gmx::as_rvec_array(x.data()), box, box_space, false);
2531 write_sto_conf(outputConfFile_.c_str(), title, atoms, gmx::as_rvec_array(x.data()), nullptr,
2534 done_symtab(&symtab);
2535 done_atom(&pdba_all);
2537 for (int chain = 0; chain < numChains; chain++)
2539 sfree(sortAtoms[chain]);
2540 sfree(localAtoms[chain]);
2548 GMX_LOG(logger.info)
2550 .appendTextFormatted("\t\t--------- PLEASE NOTE ------------");
2551 GMX_LOG(logger.info)
2553 .appendTextFormatted("You have successfully generated a topology from: %s.",
2554 inputConfFile_.c_str());
2555 if (watermodel_ != nullptr)
2557 GMX_LOG(logger.info)
2559 .appendTextFormatted("The %s force field and the %s water model are used.", ffname_,
2565 GMX_LOG(logger.info).asParagraph().appendTextFormatted("The %s force field is used.", ffname_);
2567 GMX_LOG(logger.info)
2569 .appendTextFormatted("\t\t--------- ETON ESAELP ------------");
2576 const char pdb2gmxInfo::name[] = "pdb2gmx";
2577 const char pdb2gmxInfo::shortDescription[] =
2578 "Convert coordinate files to topology and FF-compliant coordinate files";
2579 ICommandLineOptionsModulePointer pdb2gmxInfo::create()
2581 return std::make_unique<pdb2gmx>();