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