src/gromacs/topology/atoms.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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "atoms.h"
  41
  42 #include <cstdio>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/topology/atomprop.h"
  48 #include "gromacs/topology/symtab.h"
  49 #include "gromacs/utility/compare.h"
  50 #include "gromacs/utility/enumerationhelpers.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/smalloc.h"
  53 #include "gromacs/utility/txtdump.h"
  54
  55 const char* enumValueToString(ParticleType enumValue)
  56 {
  57     static constexpr gmx::EnumerationArray<ParticleType, const char*> particleTypeNames = {
  58         "Atom", "Nucleus", "Shell", "Bond", "VSite"
  59     };
  60     return particleTypeNames[enumValue];
  61 }
  62
  63 void init_atom(t_atoms* at)
  64 {
  65     at->nr          = 0;
  66     at->nres        = 0;
  67     at->atom        = nullptr;
  68     at->resinfo     = nullptr;
  69     at->atomname    = nullptr;
  70     at->atomtype    = nullptr;
  71     at->atomtypeB   = nullptr;
  72     at->pdbinfo     = nullptr;
  73     at->haveMass    = FALSE;
  74     at->haveCharge  = FALSE;
  75     at->haveType    = FALSE;
  76     at->haveBState  = FALSE;
  77     at->havePdbInfo = FALSE;
  78 }
  79
  80 void init_atomtypes(t_atomtypes* at)
  81 {
  82     at->nr         = 0;
  83     at->atomnumber = nullptr;
  84 }
  85
  86 void done_atom(t_atoms* at)
  87 {
  88     sfree(at->atom);
  89     sfree(at->resinfo);
  90     sfree(at->atomname);
  91     sfree(at->atomtype);
  92     sfree(at->atomtypeB);
  93     sfree(at->pdbinfo);
  94     init_atom(at);
  95 }
  96
  97 void done_and_delete_atoms(t_atoms* atoms)
  98 {
  99     done_atom(atoms);
 100     delete atoms;
 101 }
 102
 103 void done_atomtypes(t_atomtypes* atype)
 104 {
 105     atype->nr = 0;
 106     sfree(atype->atomnumber);
 107 }
 108
 109 void add_t_atoms(t_atoms* atoms, int natom_extra, int nres_extra)
 110 {
 111     if (natom_extra > 0)
 112     {
 113         srenew(atoms->atomname, atoms->nr + natom_extra);
 114         srenew(atoms->atom, atoms->nr + natom_extra);
 115         if (nullptr != atoms->pdbinfo)
 116         {
 117             srenew(atoms->pdbinfo, atoms->nr + natom_extra);
 118         }
 119         if (nullptr != atoms->atomtype)
 120         {
 121             srenew(atoms->atomtype, atoms->nr + natom_extra);
 122         }
 123         if (nullptr != atoms->atomtypeB)
 124         {
 125             srenew(atoms->atomtypeB, atoms->nr + natom_extra);
 126         }
 127         for (int i = atoms->nr; (i < atoms->nr + natom_extra); i++)
 128         {
 129             atoms->atomname[i] = nullptr;
 130             memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
 131             if (nullptr != atoms->pdbinfo)
 132             {
 133                 std::memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
 134             }
 135             if (nullptr != atoms->atomtype)
 136             {
 137                 atoms->atomtype[i] = nullptr;
 138             }
 139             if (nullptr != atoms->atomtypeB)
 140             {
 141                 atoms->atomtypeB[i] = nullptr;
 142             }
 143         }
 144         atoms->nr += natom_extra;
 145     }
 146     if (nres_extra > 0)
 147     {
 148         srenew(atoms->resinfo, atoms->nres + nres_extra);
 149         for (int i = atoms->nres; (i < atoms->nres + nres_extra); i++)
 150         {
 151             std::memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
 152         }
 153         atoms->nres += nres_extra;
 154     }
 155 }
 156
 157 void init_t_atoms(t_atoms* atoms, int natoms, gmx_bool bPdbinfo)
 158 {
 159     atoms->nr   = natoms;
 160     atoms->nres = 0;
 161     snew(atoms->atomname, natoms);
 162     atoms->atomtype  = nullptr;
 163     atoms->atomtypeB = nullptr;
 164     snew(atoms->resinfo, natoms);
 165     snew(atoms->atom, natoms);
 166     atoms->haveMass    = FALSE;
 167     atoms->haveCharge  = FALSE;
 168     atoms->haveType    = FALSE;
 169     atoms->haveBState  = FALSE;
 170     atoms->havePdbInfo = bPdbinfo;
 171     if (atoms->havePdbInfo)
 172     {
 173         snew(atoms->pdbinfo, natoms);
 174     }
 175     else
 176     {
 177         atoms->pdbinfo = nullptr;
 178     }
 179 }
 180
 181 void gmx_pdbinfo_init_default(t_pdbinfo* pdbinfo)
 182 {
 183     pdbinfo->type         = PdbRecordType::Atom;
 184     pdbinfo->atomnr       = 0;
 185     pdbinfo->altloc       = ' ';
 186     pdbinfo->atomnm[0]    = '\0';
 187     pdbinfo->occup        = 1.0;
 188     pdbinfo->bfac         = 0.0;
 189     pdbinfo->bAnisotropic = FALSE;
 190     std::fill(pdbinfo->uij, pdbinfo->uij + 6, 0.0);
 191 }
 192
 193 t_atoms* copy_t_atoms(const t_atoms* src)
 194 {
 195     t_atoms* dst = nullptr;
 196
 197     snew(dst, 1);
 198     init_t_atoms(dst, src->nr, (nullptr != src->pdbinfo));
 199     if (nullptr != src->atomtype)
 200     {
 201         snew(dst->atomtype, src->nr);
 202     }
 203     if (nullptr != src->atomtypeB)
 204     {
 205         snew(dst->atomtypeB, src->nr);
 206     }
 207     for (int i = 0; (i < src->nr); i++)
 208     {
 209         dst->atom[i] = src->atom[i];
 210         if (nullptr != src->pdbinfo)
 211         {
 212             dst->pdbinfo[i] = src->pdbinfo[i];
 213         }
 214         if (nullptr != src->atomname)
 215         {
 216             dst->atomname[i] = src->atomname[i];
 217         }
 218         if (nullptr != src->atomtype)
 219         {
 220             dst->atomtype[i] = src->atomtype[i];
 221         }
 222         if (nullptr != src->atomtypeB)
 223         {
 224             dst->atomtypeB[i] = src->atomtypeB[i];
 225         }
 226     }
 227     dst->haveBState  = src->haveBState;
 228     dst->haveCharge  = src->haveCharge;
 229     dst->haveMass    = src->haveMass;
 230     dst->havePdbInfo = src->havePdbInfo;
 231     dst->haveType    = src->haveType;
 232     dst->nres        = src->nres;
 233     for (int i = 0; (i < src->nres); i++)
 234     {
 235         dst->resinfo[i] = src->resinfo[i];
 236     }
 237     return dst;
 238 }
 239
 240 void t_atoms_set_resinfo(t_atoms*      atoms,
 241                          int           atom_ind,
 242                          t_symtab*     symtab,
 243                          const char*   resname,
 244                          int           resnr,
 245                          unsigned char ic,
 246                          int           chainnum,
 247                          char          chainid)
 248 {
 249     t_resinfo* ri = &atoms->resinfo[atoms->atom[atom_ind].resind];
 250     ri->name      = put_symtab(symtab, resname);
 251     ri->rtp       = nullptr;
 252     ri->nr        = resnr;
 253     ri->ic        = ic;
 254     ri->chainnum  = chainnum;
 255     ri->chainid   = chainid;
 256 }
 257
 258 static void pr_atom(FILE* fp, int indent, const char* title, const t_atom* atom, int n)
 259 {
 260     if (available(fp, atom, indent, title))
 261     {
 262         indent = pr_title_n(fp, indent, title, n);
 263         for (int i = 0; i < n; i++)
 264         {
 265             pr_indent(fp, indent);
 266             fprintf(fp,
 267                     "%s[%6d]={type=%3hu, typeB=%3hu, ptype=%8s, m=%12.5e, "
 268                     "q=%12.5e, mB=%12.5e, qB=%12.5e, resind=%5d, atomnumber=%3d}\n",
 269                     title,
 270                     i,
 271                     atom[i].type,
 272                     atom[i].typeB,
 273                     enumValueToString(atom[i].ptype),
 274                     atom[i].m,
 275                     atom[i].q,
 276                     atom[i].mB,
 277                     atom[i].qB,
 278                     atom[i].resind,
 279                     atom[i].atomnumber);
 280         }
 281     }
 282 }
 283
 284 static void pr_strings2(FILE* fp, int indent, const char* title, char*** nm, char*** nmB, int n, gmx_bool bShowNumbers)
 285 {
 286     if (available(fp, nm, indent, title))
 287     {
 288         indent = pr_title_n(fp, indent, title, n);
 289         for (int i = 0; i < n; i++)
 290         {
 291             pr_indent(fp, indent);
 292             fprintf(fp, "%s[%d]={name=\"%s\",nameB=\"%s\"}\n", title, bShowNumbers ? i : -1, *(nm[i]), *(nmB[i]));
 293         }
 294     }
 295 }
 296
 297 static void pr_resinfo(FILE* fp, int indent, const char* title, const t_resinfo* resinfo, int n, gmx_bool bShowNumbers)
 298 {
 299     if (available(fp, resinfo, indent, title))
 300     {
 301         indent = pr_title_n(fp, indent, title, n);
 302         for (int i = 0; i < n; i++)
 303         {
 304             pr_indent(fp, indent);
 305             fprintf(fp,
 306                     "%s[%d]={name=\"%s\", nr=%d, ic='%c'}\n",
 307                     title,
 308                     bShowNumbers ? i : -1,
 309                     *(resinfo[i].name),
 310                     resinfo[i].nr,
 311                     (resinfo[i].ic == '\0') ? ' ' : resinfo[i].ic);
 312         }
 313     }
 314 }
 315
 316 void pr_atoms(FILE* fp, int indent, const char* title, const t_atoms* atoms, gmx_bool bShownumbers)
 317 {
 318     if (available(fp, atoms, indent, title))
 319     {
 320         indent = pr_title(fp, indent, title);
 321         pr_atom(fp, indent, "atom", atoms->atom, atoms->nr);
 322         pr_strings(fp, indent, "atom", atoms->atomname, atoms->nr, bShownumbers);
 323         pr_strings2(fp, indent, "type", atoms->atomtype, atoms->atomtypeB, atoms->nr, bShownumbers);
 324         pr_resinfo(fp, indent, "residue", atoms->resinfo, atoms->nres, bShownumbers);
 325     }
 326 }
 327
 328
 329 void pr_atomtypes(FILE* fp, int indent, const char* title, const t_atomtypes* atomtypes, gmx_bool bShowNumbers)
 330 {
 331     if (available(fp, atomtypes, indent, title))
 332     {
 333         indent = pr_title(fp, indent, title);
 334         for (int i = 0; i < atomtypes->nr; i++)
 335         {
 336             pr_indent(fp, indent);
 337             fprintf(fp, "atomtype[%3d]={atomnumber=%4d}\n", bShowNumbers ? i : -1, atomtypes->atomnumber[i]);
 338         }
 339     }
 340 }
 341
 342 static void compareAtom(FILE* fp, int index, const t_atom* a1, const t_atom* a2, real relativeTolerance, real absoluteTolerance)
 343 {
 344     if (a2)
 345     {
 346         cmp_us(fp, "atom.type", index, a1->type, a2->type);
 347         cmpEnum<ParticleType>(fp, "atom.ptype", a1->ptype, a2->ptype);
 348         cmp_int(fp, "atom.resind", index, a1->resind, a2->resind);
 349         cmp_int(fp, "atom.atomnumber", index, a1->atomnumber, a2->atomnumber);
 350         cmp_real(fp, "atom.m", index, a1->m, a2->m, relativeTolerance, absoluteTolerance);
 351         cmp_real(fp, "atom.q", index, a1->q, a2->q, relativeTolerance, absoluteTolerance);
 352         cmp_us(fp, "atom.typeB", index, a1->typeB, a2->typeB);
 353         cmp_real(fp, "atom.mB", index, a1->mB, a2->mB, relativeTolerance, absoluteTolerance);
 354         cmp_real(fp, "atom.qB", index, a1->qB, a2->qB, relativeTolerance, absoluteTolerance);
 355         cmp_str(fp, "elem", index, a1->elem, a2->elem);
 356     }
 357     else
 358     {
 359         cmp_us(fp, "atom.type", index, a1->type, a1->typeB);
 360         cmp_real(fp, "atom.m", index, a1->m, a1->mB, relativeTolerance, absoluteTolerance);
 361         cmp_real(fp, "atom.q", index, a1->q, a1->qB, relativeTolerance, absoluteTolerance);
 362     }
 363 }
 364
 365 static void compareResinfo(FILE* fp, int residue, const t_resinfo& r1, const t_resinfo& r2)
 366 {
 367     fprintf(fp, "comparing t_resinfo\n");
 368     cmp_str(fp, "name", residue, *r1.name, *r2.name);
 369     cmp_int(fp, "nr", residue, r1.nr, r2.nr);
 370     cmp_uc(fp, "ic", residue, r1.ic, r2.ic);
 371     cmp_int(fp, "chainnum", residue, r1.chainnum, r2.chainnum);
 372     cmp_uc(fp, "chainid", residue, r1.chainid, r2.chainid);
 373     if ((r1.rtp || r2.rtp) && (!r1.rtp || !r2.rtp))
 374     {
 375         fprintf(fp, "rtp info is present in topology %d but not in the other\n", r1.rtp ? 1 : 2);
 376     }
 377     if (r1.rtp && r2.rtp)
 378     {
 379         cmp_str(fp, "rtp", residue, *r1.rtp, *r2.rtp);
 380     }
 381 }
 382
 383 static void comparePdbinfo(FILE*            fp,
 384                            int              pdb,
 385                            const t_pdbinfo& pdb1,
 386                            const t_pdbinfo& pdb2,
 387                            real             relativeTolerance,
 388                            real             absoluteTolerance)
 389 {
 390     fprintf(fp, "comparing t_pdbinfo\n");
 391     cmpEnum<PdbRecordType>(fp, "type", pdb1.type, pdb2.type);
 392     cmp_int(fp, "atomnr", pdb, pdb1.atomnr, pdb2.atomnr);
 393     cmp_uc(fp, "altloc", pdb, pdb1.altloc, pdb2.altloc);
 394     cmp_str(fp, "atomnm", pdb, pdb1.atomnm, pdb2.atomnm);
 395     cmp_real(fp, "occup", pdb, pdb1.occup, pdb2.occup, relativeTolerance, absoluteTolerance);
 396     cmp_real(fp, "bfac", pdb, pdb1.bfac, pdb2.bfac, relativeTolerance, absoluteTolerance);
 397     cmp_bool(fp, "bAnistropic", pdb, pdb1.bAnisotropic, pdb2.bAnisotropic);
 398     for (int i = 0; i < 6; i++)
 399     {
 400         std::string buf = gmx::formatString("uij[%d]", i);
 401         cmp_int(fp, buf.c_str(), pdb, pdb1.uij[i], pdb2.uij[i]);
 402     }
 403 }
 404
 405
 406 void compareAtoms(FILE* fp, const t_atoms* a1, const t_atoms* a2, real relativeTolerance, real absoluteTolerance)
 407 {
 408     fprintf(fp, "comparing atoms\n");
 409
 410     if (a2)
 411     {
 412         cmp_int(fp, "atoms->nr", -1, a1->nr, a2->nr);
 413         cmp_int(fp, "atoms->nres", -1, a1->nres, a2->nres);
 414         cmp_bool(fp, "atoms->haveMass", -1, a1->haveMass, a2->haveMass);
 415         cmp_bool(fp, "atoms->haveCharge", -1, a1->haveCharge, a2->haveCharge);
 416         cmp_bool(fp, "atoms->haveType", -1, a1->haveType, a2->haveType);
 417         cmp_bool(fp, "atoms->haveBState", -1, a1->haveBState, a2->haveBState);
 418         cmp_bool(fp, "atoms->havePdbInfo", -1, a1->havePdbInfo, a2->havePdbInfo);
 419         for (int i = 0; i < std::min(a1->nr, a2->nr); i++)
 420         {
 421             compareAtom(fp, i, &(a1->atom[i]), &(a2->atom[i]), relativeTolerance, absoluteTolerance);
 422             if (a1->atomname && a2->atomname)
 423             {
 424                 cmp_str(fp, "atomname", i, *a1->atomname[i], *a2->atomname[i]);
 425             }
 426             if (a1->havePdbInfo && a2->havePdbInfo)
 427             {
 428                 comparePdbinfo(fp, i, a1->pdbinfo[i], a2->pdbinfo[i], relativeTolerance, absoluteTolerance);
 429             }
 430             if (a1->haveType && a2->haveType)
 431             {
 432                 cmp_str(fp, "atomtype", i, *a1->atomtype[i], *a2->atomtype[i]);
 433             }
 434             if (a1->haveBState && a2->haveBState)
 435             {
 436                 cmp_str(fp, "atomtypeB", i, *a1->atomtypeB[i], *a2->atomtypeB[i]);
 437             }
 438         }
 439         for (int i = 0; i < std::min(a1->nres, a2->nres); i++)
 440         {
 441             compareResinfo(fp, i, a1->resinfo[i], a2->resinfo[i]);
 442         }
 443     }
 444     else
 445     {
 446         for (int i = 0; (i < a1->nr); i++)
 447         {
 448             compareAtom(fp, i, &(a1->atom[i]), nullptr, relativeTolerance, absoluteTolerance);
 449         }
 450     }
 451 }
 452
 453 void atomsSetMassesBasedOnNames(t_atoms* atoms, gmx_bool printMissingMasses)
 454 {
 455     if (atoms->haveMass)
 456     {
 457         /* We could decide to anyhow assign then or generate a fatal error,
 458          * but it's probably most useful to keep the masses we have.
 459          */
 460         return;
 461     }
 462
 463     int maxWarn = (printMissingMasses ? 10 : 0);
 464     int numWarn = 0;
 465
 466     AtomProperties aps;
 467
 468     bool haveMass = true;
 469     for (int i = 0; i < atoms->nr; i++)
 470     {
 471         if (!aps.setAtomProperty(epropMass,
 472                                  *atoms->resinfo[atoms->atom[i].resind].name,
 473                                  *atoms->atomname[i],
 474                                  &atoms->atom[i].m))
 475         {
 476             haveMass = false;
 477
 478             if (numWarn < maxWarn)
 479             {
 480                 fprintf(stderr,
 481                         "Can not find mass in database for atom %s in residue %d %s\n",
 482                         *atoms->atomname[i],
 483                         atoms->resinfo[atoms->atom[i].resind].nr,
 484                         *atoms->resinfo[atoms->atom[i].resind].name);
 485                 numWarn++;
 486             }
 487             else
 488             {
 489                 break;
 490             }
 491         }
 492     }
 493     atoms->haveMass = haveMass;
 494 }