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, 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 "atoms.h"
  40
  41 #include <cstdio>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/topology/atomprop.h"
  47 #include "gromacs/topology/symtab.h"
  48 #include "gromacs/utility/compare.h"
  49 #include "gromacs/utility/fatalerror.h"
  50 #include "gromacs/utility/smalloc.h"
  51 #include "gromacs/utility/txtdump.h"
  52
  53 const char *ptype_str[eptNR+1] = {
  54     "Atom", "Nucleus", "Shell", "Bond", "VSite", nullptr
  55 };
  56
  57 void init_atom(t_atoms *at)
  58 {
  59     at->nr          = 0;
  60     at->nres        = 0;
  61     at->atom        = nullptr;
  62     at->resinfo     = nullptr;
  63     at->atomname    = nullptr;
  64     at->atomtype    = nullptr;
  65     at->atomtypeB   = nullptr;
  66     at->pdbinfo     = nullptr;
  67     at->haveMass    = FALSE;
  68     at->haveCharge  = FALSE;
  69     at->haveType    = FALSE;
  70     at->haveBState  = FALSE;
  71     at->havePdbInfo = FALSE;
  72 }
  73
  74 void init_atomtypes(t_atomtypes *at)
  75 {
  76     at->nr         = 0;
  77     at->radius     = nullptr;
  78     at->vol        = nullptr;
  79     at->atomnumber = nullptr;
  80     at->gb_radius  = nullptr;
  81     at->S_hct      = nullptr;
  82 }
  83
  84 void done_atom(t_atoms *at)
  85 {
  86     sfree(at->atom);
  87     sfree(at->resinfo);
  88     sfree(at->atomname);
  89     sfree(at->atomtype);
  90     sfree(at->atomtypeB);
  91     sfree(at->pdbinfo);
  92     init_atom(at);
  93 }
  94
  95 void done_atomtypes(t_atomtypes *atype)
  96 {
  97     atype->nr = 0;
  98     sfree(atype->radius);
  99     sfree(atype->vol);
 100     sfree(atype->surftens);
 101     sfree(atype->atomnumber);
 102     sfree(atype->gb_radius);
 103     sfree(atype->S_hct);
 104 }
 105
 106 void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
 107 {
 108     int i;
 109
 110     if (natom_extra > 0)
 111     {
 112         srenew(atoms->atomname, atoms->nr+natom_extra);
 113         srenew(atoms->atom, atoms->nr+natom_extra);
 114         if (nullptr != atoms->pdbinfo)
 115         {
 116             srenew(atoms->pdbinfo, atoms->nr+natom_extra);
 117         }
 118         if (nullptr != atoms->atomtype)
 119         {
 120             srenew(atoms->atomtype, atoms->nr+natom_extra);
 121         }
 122         if (nullptr != atoms->atomtypeB)
 123         {
 124             srenew(atoms->atomtypeB, atoms->nr+natom_extra);
 125         }
 126         for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
 127         {
 128             atoms->atomname[i] = nullptr;
 129             memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
 130             if (nullptr != atoms->pdbinfo)
 131             {
 132                 std::memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
 133             }
 134             if (nullptr != atoms->atomtype)
 135             {
 136                 atoms->atomtype[i] = nullptr;
 137             }
 138             if (nullptr != atoms->atomtypeB)
 139             {
 140                 atoms->atomtypeB[i] = nullptr;
 141             }
 142         }
 143         atoms->nr += natom_extra;
 144     }
 145     if (nres_extra > 0)
 146     {
 147         srenew(atoms->resinfo, atoms->nres+nres_extra);
 148         for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
 149         {
 150             std::memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
 151         }
 152         atoms->nres += nres_extra;
 153     }
 154 }
 155
 156 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
 157 {
 158     atoms->nr   = natoms;
 159     atoms->nres = 0;
 160     snew(atoms->atomname, natoms);
 161     atoms->atomtype  = nullptr;
 162     atoms->atomtypeB = nullptr;
 163     snew(atoms->resinfo, natoms);
 164     snew(atoms->atom, natoms);
 165     atoms->haveMass    = FALSE;
 166     atoms->haveCharge  = FALSE;
 167     atoms->haveType    = FALSE;
 168     atoms->haveBState  = FALSE;
 169     atoms->havePdbInfo = bPdbinfo;
 170     if (atoms->havePdbInfo)
 171     {
 172         snew(atoms->pdbinfo, natoms);
 173     }
 174     else
 175     {
 176         atoms->pdbinfo = nullptr;
 177     }
 178 }
 179
 180 void gmx_pdbinfo_init_default(t_pdbinfo *pdbinfo)
 181 {
 182     pdbinfo->type         = epdbATOM;
 183     pdbinfo->atomnr       = 0;
 184     pdbinfo->altloc       = ' ';
 185     pdbinfo->atomnm[0]    = '\0';
 186     pdbinfo->occup        = 1.0;
 187     pdbinfo->bfac         = 0.0;
 188     pdbinfo->bAnisotropic = FALSE;
 189     std::fill(pdbinfo->uij, pdbinfo->uij+6, 0.0);
 190 }
 191
 192 t_atoms *copy_t_atoms(const t_atoms *src)
 193 {
 194     t_atoms *dst;
 195     int      i;
 196
 197     snew(dst, 1);
 198     init_t_atoms(dst, src->nr, (nullptr != src->pdbinfo));
 199     dst->nr = src->nr;
 200     if (nullptr != src->atomname)
 201     {
 202         snew(dst->atomname, src->nr);
 203     }
 204     if (nullptr != src->atomtype)
 205     {
 206         snew(dst->atomtype, src->nr);
 207     }
 208     if (nullptr != src->atomtypeB)
 209     {
 210         snew(dst->atomtypeB, src->nr);
 211     }
 212     for (i = 0; (i < src->nr); i++)
 213     {
 214         dst->atom[i] = src->atom[i];
 215         if (nullptr != src->pdbinfo)
 216         {
 217             dst->pdbinfo[i] = src->pdbinfo[i];
 218         }
 219         if (nullptr != src->atomname)
 220         {
 221             dst->atomname[i]  = src->atomname[i];
 222         }
 223         if (nullptr != src->atomtype)
 224         {
 225             dst->atomtype[i] = src->atomtype[i];
 226         }
 227         if (nullptr != src->atomtypeB)
 228         {
 229             dst->atomtypeB[i] = src->atomtypeB[i];
 230         }
 231     }
 232     dst->nres = src->nres;
 233     for (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, int atom_ind, t_symtab *symtab,
 241                          const char *resname, int resnr, unsigned char ic,
 242                          int chainnum, char chainid)
 243 {
 244     t_resinfo *ri;
 245
 246     ri           = &atoms->resinfo[atoms->atom[atom_ind].resind];
 247     ri->name     = put_symtab(symtab, resname);
 248     ri->rtp      = nullptr;
 249     ri->nr       = resnr;
 250     ri->ic       = ic;
 251     ri->chainnum = chainnum;
 252     ri->chainid  = chainid;
 253 }
 254
 255 static void pr_atom(FILE *fp, int indent, const char *title, const t_atom *atom, int n)
 256 {
 257     int i;
 258
 259     if (available(fp, atom, indent, title))
 260     {
 261         indent = pr_title_n(fp, indent, title, n);
 262         for (i = 0; i < n; i++)
 263         {
 264             pr_indent(fp, indent);
 265             fprintf(fp, "%s[%6d]={type=%3hu, typeB=%3hu, ptype=%8s, m=%12.5e, "
 266                     "q=%12.5e, mB=%12.5e, qB=%12.5e, resind=%5d, atomnumber=%3d}\n",
 267                     title, i, atom[i].type, atom[i].typeB, ptype_str[atom[i].ptype],
 268                     atom[i].m, atom[i].q, atom[i].mB, atom[i].qB,
 269                     atom[i].resind, atom[i].atomnumber);
 270         }
 271     }
 272 }
 273
 274 static void pr_strings2(FILE *fp, int indent, const char *title,
 275                         char ***nm, char ***nmB, int n, gmx_bool bShowNumbers)
 276 {
 277     int i;
 278
 279     if (available(fp, nm, indent, title))
 280     {
 281         indent = pr_title_n(fp, indent, title, n);
 282         for (i = 0; i < n; i++)
 283         {
 284             pr_indent(fp, indent);
 285             fprintf(fp, "%s[%d]={name=\"%s\",nameB=\"%s\"}\n",
 286                     title, bShowNumbers ? i : -1, *(nm[i]), *(nmB[i]));
 287         }
 288     }
 289 }
 290
 291 static void pr_resinfo(FILE *fp, int indent, const char *title, const t_resinfo *resinfo, int n,
 292                        gmx_bool bShowNumbers)
 293 {
 294     int i;
 295
 296     if (available(fp, resinfo, indent, title))
 297     {
 298         indent = pr_title_n(fp, indent, title, n);
 299         for (i = 0; i < n; i++)
 300         {
 301             pr_indent(fp, indent);
 302             fprintf(fp, "%s[%d]={name=\"%s\", nr=%d, ic='%c'}\n",
 303                     title, bShowNumbers ? i : -1,
 304                     *(resinfo[i].name), resinfo[i].nr,
 305                     (resinfo[i].ic == '\0') ? ' ' : resinfo[i].ic);
 306         }
 307     }
 308 }
 309
 310 void pr_atoms(FILE *fp, int indent, const char *title, const t_atoms *atoms,
 311               gmx_bool bShownumbers)
 312 {
 313     if (available(fp, atoms, indent, title))
 314     {
 315         indent = pr_title(fp, indent, title);
 316         pr_atom(fp, indent, "atom", atoms->atom, atoms->nr);
 317         pr_strings(fp, indent, "atom", atoms->atomname, atoms->nr, bShownumbers);
 318         pr_strings2(fp, indent, "type", atoms->atomtype, atoms->atomtypeB, atoms->nr, bShownumbers);
 319         pr_resinfo(fp, indent, "residue", atoms->resinfo, atoms->nres, bShownumbers);
 320     }
 321 }
 322
 323
 324 void pr_atomtypes(FILE *fp, int indent, const char *title, const t_atomtypes *atomtypes,
 325                   gmx_bool bShowNumbers)
 326 {
 327     int i;
 328     if (available(fp, atomtypes, indent, title))
 329     {
 330         indent = pr_title(fp, indent, title);
 331         for (i = 0; i < atomtypes->nr; i++)
 332         {
 333             pr_indent(fp, indent);
 334             fprintf(fp,
 335                     "atomtype[%3d]={radius=%12.5e, volume=%12.5e, gb_radius=%12.5e, surftens=%12.5e, atomnumber=%4d, S_hct=%12.5e)}\n",
 336                     bShowNumbers ? i : -1, atomtypes->radius[i], atomtypes->vol[i],
 337                     atomtypes->gb_radius[i],
 338                     atomtypes->surftens[i], atomtypes->atomnumber[i], atomtypes->S_hct[i]);
 339         }
 340     }
 341 }
 342
 343 static void cmp_atom(FILE *fp, int index, const t_atom *a1, const t_atom *a2, real ftol, real abstol)
 344 {
 345     if (a2)
 346     {
 347         cmp_us(fp, "atom.type", index, a1->type, a2->type);
 348         cmp_us(fp, "atom.ptype", index, a1->ptype, a2->ptype);
 349         cmp_int(fp, "atom.resind", index, a1->resind, a2->resind);
 350         cmp_int(fp, "atom.atomnumber", index, a1->atomnumber, a2->atomnumber);
 351         cmp_real(fp, "atom.m", index, a1->m, a2->m, ftol, abstol);
 352         cmp_real(fp, "atom.q", index, a1->q, a2->q, ftol, abstol);
 353         cmp_us(fp, "atom.typeB", index, a1->typeB, a2->typeB);
 354         cmp_real(fp, "atom.mB", index, a1->mB, a2->mB, ftol, abstol);
 355         cmp_real(fp, "atom.qB", index, a1->qB, a2->qB, ftol, abstol);
 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, ftol, abstol);
 361         cmp_real(fp, "atom.q", index, a1->q, a1->qB, ftol, abstol);
 362     }
 363 }
 364
 365 void cmp_atoms(FILE *fp, const t_atoms *a1, const t_atoms *a2, real ftol, real abstol)
 366 {
 367     int i;
 368
 369     fprintf(fp, "comparing atoms\n");
 370
 371     if (a2)
 372     {
 373         cmp_int(fp, "atoms->nr", -1, a1->nr, a2->nr);
 374         for (i = 0; i < std::min(a1->nr, a2->nr); i++)
 375         {
 376             cmp_atom(fp, i, &(a1->atom[i]), &(a2->atom[i]), ftol, abstol);
 377         }
 378     }
 379     else
 380     {
 381         for (i = 0; (i < a1->nr); i++)
 382         {
 383             cmp_atom(fp, i, &(a1->atom[i]), nullptr, ftol, abstol);
 384         }
 385     }
 386 }
 387
 388 void atomsSetMassesBasedOnNames(t_atoms *atoms, gmx_bool printMissingMasses)
 389 {
 390     if (atoms->haveMass)
 391     {
 392         /* We could decide to anyhow assign then or generate a fatal error,
 393          * but it's probably most useful to keep the masses we have.
 394          */
 395         return;
 396     }
 397
 398     int            maxWarn  = (printMissingMasses ? 10 : 0);
 399     int            numWarn  = 0;
 400
 401     gmx_atomprop_t aps      = gmx_atomprop_init();
 402
 403     gmx_bool       haveMass = TRUE;
 404     for (int i = 0; i < atoms->nr; i++)
 405     {
 406         if (!gmx_atomprop_query(aps, epropMass,
 407                                 *atoms->resinfo[atoms->atom[i].resind].name,
 408                                 *atoms->atomname[i],
 409                                 &atoms->atom[i].m))
 410         {
 411             haveMass = FALSE;
 412
 413             if (numWarn < maxWarn)
 414             {
 415                 fprintf(stderr, "Can not find mass in database for atom %s in residue %d %s\n",
 416                         *atoms->atomname[i],
 417                         atoms->resinfo[atoms->atom[i].resind].nr,
 418                         *atoms->resinfo[atoms->atom[i].resind].name);
 419                 numWarn++;
 420             }
 421             else
 422             {
 423                 break;
 424             }
 425         }
 426     }
 427     atoms->haveMass = haveMass;
 428
 429     gmx_atomprop_destroy(aps);
 430 }