Bug Summary

File:gromacs/gmxana/gmx_editconf.c
Location:line 777, column 17
Description:Array access results in a null pointer dereference

Annotated Source Code

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, 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#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <string.h>
43
44#include "gromacs/fileio/pdbio.h"
45#include "gromacs/fileio/confio.h"
46#include "symtab.h"
47#include "gromacs/utility/smalloc.h"
48#include "macros.h"
49#include "gromacs/commandline/pargs.h"
50#include "gromacs/fileio/strdb.h"
51#include "index.h"
52#include "gromacs/math/vec.h"
53#include "typedefs.h"
54#include "gromacs/gmxlib/conformation-utilities.h"
55#include "physics.h"
56#include "atomprop.h"
57#include "gromacs/fileio/tpxio.h"
58#include "gromacs/fileio/trxio.h"
59#include "pbc.h"
60#include "princ.h"
61#include "txtdump.h"
62#include "viewit.h"
63#include "rmpbc.h"
64#include "gmx_ana.h"
65
66#include "gromacs/utility/fatalerror.h"
67
68typedef struct
69{
70 char sanm[12];
71 int natm;
72 int nw;
73 char anm[6][12];
74} t_simat;
75
76typedef struct
77{
78 char reso[12];
79 char resn[12];
80 int nsatm;
81 t_simat sat[3];
82} t_simlist;
83
84real calc_mass(t_atoms *atoms, gmx_bool bGetMass, gmx_atomprop_t aps)
85{
86 real tmass;
87 int i;
88
89 tmass = 0;
90 for (i = 0; (i < atoms->nr); i++)
91 {
92 if (bGetMass)
93 {
94 gmx_atomprop_query(aps, epropMass,
95 *atoms->resinfo[atoms->atom[i].resind].name,
96 *atoms->atomname[i], &(atoms->atom[i].m));
97 }
98 tmass += atoms->atom[i].m;
99 }
100
101 return tmass;
102}
103
104real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
105 rvec max, gmx_bool bDiam)
106{
107 real diam2, d;
108 char *grpnames;
109 int ii, i, j;
110
111 clear_rvec(geom_center);
112 diam2 = 0;
113 if (isize == 0)
114 {
115 clear_rvec(min);
116 clear_rvec(max);
117 }
118 else
119 {
120 if (index)
121 {
122 ii = index[0];
123 }
124 else
125 {
126 ii = 0;
127 }
128 for (j = 0; j < DIM3; j++)
129 {
130 min[j] = max[j] = x[ii][j];
131 }
132 for (i = 0; i < isize; i++)
133 {
134 if (index)
135 {
136 ii = index[i];
137 }
138 else
139 {
140 ii = i;
141 }
142 rvec_inc(geom_center, x[ii]);
143 for (j = 0; j < DIM3; j++)
144 {
145 if (x[ii][j] < min[j])
146 {
147 min[j] = x[ii][j];
148 }
149 if (x[ii][j] > max[j])
150 {
151 max[j] = x[ii][j];
152 }
153 }
154 if (bDiam)
155 {
156 if (index)
157 {
158 for (j = i + 1; j < isize; j++)
159 {
160 d = distance2(x[ii], x[index[j]]);
161 diam2 = max(d, diam2)(((d) > (diam2)) ? (d) : (diam2) );
162 }
163 }
164 else
165 {
166 for (j = i + 1; j < isize; j++)
167 {
168 d = distance2(x[i], x[j]);
169 diam2 = max(d, diam2)(((d) > (diam2)) ? (d) : (diam2) );
170 }
171 }
172 }
173 }
174 svmul(1.0 / isize, geom_center, geom_center);
175 }
176
177 return sqrt(diam2);
178}
179
180void center_conf(int natom, rvec *x, rvec center, rvec geom_cent)
181{
182 int i;
183 rvec shift;
184
185 rvec_sub(center, geom_cent, shift);
186
187 printf(" shift :%7.3f%7.3f%7.3f (nm)\n", shift[XX0], shift[YY1],
188 shift[ZZ2]);
189
190 for (i = 0; (i < natom); i++)
191 {
192 rvec_inc(x[i], shift);
193 }
194}
195
196void scale_conf(int natom, rvec x[], matrix box, rvec scale)
197{
198 int i, j;
199
200 for (i = 0; i < natom; i++)
201 {
202 for (j = 0; j < DIM3; j++)
203 {
204 x[i][j] *= scale[j];
205 }
206 }
207 for (i = 0; i < DIM3; i++)
208 {
209 for (j = 0; j < DIM3; j++)
210 {
211 box[i][j] *= scale[j];
212 }
213 }
214}
215
216void read_bfac(const char *fn, int *n_bfac, double **bfac_val, int **bfac_nr)
217{
218 int i;
219 char **bfac_lines;
220
221 *n_bfac = get_lines(fn, &bfac_lines);
222 snew(*bfac_val, *n_bfac)(*bfac_val) = save_calloc("*bfac_val", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 222, (*n_bfac), sizeof(*(*bfac_val)));
223 snew(*bfac_nr, *n_bfac)(*bfac_nr) = save_calloc("*bfac_nr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 223, (*n_bfac), sizeof(*(*bfac_nr)));
224 fprintf(stderrstderr, "Reading %d B-factors from %s\n", *n_bfac, fn);
225 for (i = 0; (i < *n_bfac); i++)
226 {
227 /*fprintf(stderr, "Line %d: %s",i,bfac_lines[i]);*/
228 sscanf(bfac_lines[i], "%d %lf", &(*bfac_nr)[i], &(*bfac_val)[i]);
229 /*fprintf(stderr," nr %d val %g\n",(*bfac_nr)[i],(*bfac_val)[i]);*/
230 }
231
232}
233
234void set_pdb_conf_bfac(int natoms, int nres, t_atoms *atoms, int n_bfac,
235 double *bfac, int *bfac_nr, gmx_bool peratom)
236{
237 FILE *out;
238 real bfac_min, bfac_max;
239 int i, n;
240 gmx_bool found;
241
242 bfac_max = -1e10;
243 bfac_min = 1e10;
244 for (i = 0; (i < n_bfac); i++)
245 {
246 if (bfac_nr[i] - 1 >= atoms->nres)
247 {
248 peratom = TRUE1;
249 }
250 /* if ((bfac_nr[i]-1<0) || (bfac_nr[i]-1>=atoms->nr))
251 gmx_fatal(FARGS,"Index of B-Factor %d is out of range: %d (%g)",
252 i+1,bfac_nr[i],bfac[i]); */
253 if (bfac[i] > bfac_max)
254 {
255 bfac_max = bfac[i];
256 }
257 if (bfac[i] < bfac_min)
258 {
259 bfac_min = bfac[i];
260 }
261 }
262 while ((bfac_max > 99.99) || (bfac_min < -99.99))
263 {
264 fprintf(stderrstderr,
265 "Range of values for B-factors too large (min %g, max %g) "
266 "will scale down a factor 10\n", bfac_min, bfac_max);
267 for (i = 0; (i < n_bfac); i++)
268 {
269 bfac[i] /= 10;
270 }
271 bfac_max /= 10;
272 bfac_min /= 10;
273 }
274 while ((fabs(bfac_max) < 0.5) && (fabs(bfac_min) < 0.5))
275 {
276 fprintf(stderrstderr,
277 "Range of values for B-factors too small (min %g, max %g) "
278 "will scale up a factor 10\n", bfac_min, bfac_max);
279 for (i = 0; (i < n_bfac); i++)
280 {
281 bfac[i] *= 10;
282 }
283 bfac_max *= 10;
284 bfac_min *= 10;
285 }
286
287 for (i = 0; (i < natoms); i++)
288 {
289 atoms->pdbinfo[i].bfac = 0;
290 }
291
292 if (!peratom)
293 {
294 fprintf(stderrstderr, "Will attach %d B-factors to %d residues\n", n_bfac,
295 nres);
296 for (i = 0; (i < n_bfac); i++)
297 {
298 found = FALSE0;
299 for (n = 0; (n < natoms); n++)
300 {
301 if (bfac_nr[i] == atoms->resinfo[atoms->atom[n].resind].nr)
302 {
303 atoms->pdbinfo[n].bfac = bfac[i];
304 found = TRUE1;
305 }
306 }
307 if (!found)
308 {
309 gmx_warning("Residue nr %d not found\n", bfac_nr[i]);
310 }
311 }
312 }
313 else
314 {
315 fprintf(stderrstderr, "Will attach %d B-factors to %d atoms\n", n_bfac,
316 natoms);
317 for (i = 0; (i < n_bfac); i++)
318 {
319 atoms->pdbinfo[bfac_nr[i] - 1].bfac = bfac[i];
320 }
321 }
322}
323
324void pdb_legend(FILE *out, int natoms, int nres, t_atoms *atoms, rvec x[])
325{
326 real bfac_min, bfac_max, xmin, ymin, zmin;
327 int i;
328 int space = ' ';
329
330 bfac_max = -1e10;
331 bfac_min = 1e10;
332 xmin = 1e10;
333 ymin = 1e10;
334 zmin = 1e10;
335 for (i = 0; (i < natoms); i++)
336 {
337 xmin = min(xmin, x[i][XX])(((xmin) < (x[i][0])) ? (xmin) : (x[i][0]) );
338 ymin = min(ymin, x[i][YY])(((ymin) < (x[i][1])) ? (ymin) : (x[i][1]) );
339 zmin = min(zmin, x[i][ZZ])(((zmin) < (x[i][2])) ? (zmin) : (x[i][2]) );
340 bfac_min = min(bfac_min, atoms->pdbinfo[i].bfac)(((bfac_min) < (atoms->pdbinfo[i].bfac)) ? (bfac_min) :
(atoms->pdbinfo[i].bfac) );
341 bfac_max = max(bfac_max, atoms->pdbinfo[i].bfac)(((bfac_max) > (atoms->pdbinfo[i].bfac)) ? (bfac_max) :
(atoms->pdbinfo[i].bfac) );
342 }
343 fprintf(stderrstderr, "B-factors range from %g to %g\n", bfac_min, bfac_max);
344 for (i = 1; (i < 12); i++)
345 {
346 fprintf(out,
347 "%-6s%5u %-4.4s%3.3s %c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f\n",
348 "ATOM ", natoms + 1 + i, "CA", "LEG", space, nres + 1, space,
349 (xmin + (i * 0.12)) * 10, ymin * 10, zmin * 10, 1.0, bfac_min
350 + ((i - 1.0) * (bfac_max - bfac_min) / 10));
351 }
352}
353
354void visualize_images(const char *fn, int ePBC, matrix box)
355{
356 t_atoms atoms;
357 rvec *img;
358 char *c, *ala;
359 int nat, i;
360
361 nat = NTRICIMG14 + 1;
362 init_t_atoms(&atoms, nat, FALSE0);
363 atoms.nr = nat;
364 snew(img, nat)(img) = save_calloc("img", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 364, (nat), sizeof(*(img)));
365 /* FIXME: Constness should not be cast away */
366 c = (char *) "C";
367 ala = (char *) "ALA";
368 for (i = 0; i < nat; i++)
369 {
370 atoms.atomname[i] = &c;
371 atoms.atom[i].resind = i;
372 atoms.resinfo[i].name = &ala;
373 atoms.resinfo[i].nr = i + 1;
374 atoms.resinfo[i].chainid = 'A' + i / NCUCVERT24;
375 }
376 calc_triclinic_images(box, img + 1);
377
378 write_sto_conf(fn, "Images", &atoms, img, NULL((void*)0), ePBC, box);
379
380 free_t_atoms(&atoms, FALSE0);
381 sfree(img)save_free("img", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 381, (img));
382}
383
384void visualize_box(FILE *out, int a0, int r0, matrix box, rvec gridsize)
385{
386 int *edge;
387 rvec *vert, shift;
388 int nx, ny, nz, nbox, nat;
389 int i, j, x, y, z;
390 int rectedge[24] =
391 {
392 0, 1, 1, 3, 3, 2, 0, 2, 0, 4, 1, 5, 3, 7, 2, 6, 4, 5, 5, 7, 7, 6, 6,
393 4
394 };
395
396 a0++;
397 r0++;
398
399 nx = (int) (gridsize[XX0] + 0.5);
400 ny = (int) (gridsize[YY1] + 0.5);
401 nz = (int) (gridsize[ZZ2] + 0.5);
402 nbox = nx * ny * nz;
403 if (TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0))
404 {
405 nat = nbox * NCUCVERT24;
406 snew(vert, nat)(vert) = save_calloc("vert", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 406, (nat), sizeof(*(vert)));
407 calc_compact_unitcell_vertices(ecenterDEF, box, vert);
408 j = 0;
409 for (z = 0; z < nz; z++)
410 {
411 for (y = 0; y < ny; y++)
412 {
413 for (x = 0; x < nx; x++)
414 {
415 for (i = 0; i < DIM3; i++)
416 {
417 shift[i] = x * box[0][i] + y * box[1][i] + z
418 * box[2][i];
419 }
420 for (i = 0; i < NCUCVERT24; i++)
421 {
422 rvec_add(vert[i], shift, vert[j]);
423 j++;
424 }
425 }
426 }
427 }
428
429 for (i = 0; i < nat; i++)
430 {
431 fprintf(out, get_pdbformat(), "ATOM", a0 + i, "C", "BOX", 'K' + i
432 / NCUCVERT24, r0 + i, 10 * vert[i][XX0], 10 * vert[i][YY1], 10
433 * vert[i][ZZ2]);
434 fprintf(out, "\n");
435 }
436
437 edge = compact_unitcell_edges();
438 for (j = 0; j < nbox; j++)
439 {
440 for (i = 0; i < NCUCEDGE36; i++)
441 {
442 fprintf(out, "CONECT%5d%5d\n", a0 + j * NCUCVERT24 + edge[2 * i],
443 a0 + j * NCUCVERT24 + edge[2 * i + 1]);
444 }
445 }
446
447 sfree(vert)save_free("vert", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 447, (vert));
448 }
449 else
450 {
451 i = 0;
452 for (z = 0; z <= 1; z++)
453 {
454 for (y = 0; y <= 1; y++)
455 {
456 for (x = 0; x <= 1; x++)
457 {
458 fprintf(out, get_pdbformat(), "ATOM", a0 + i, "C", "BOX", 'K' + i
459 / 8, r0 + i, x * 10 * box[XX0][XX0],
460 y * 10 * box[YY1][YY1], z * 10 * box[ZZ2][ZZ2]);
461 fprintf(out, "\n");
462 i++;
463 }
464 }
465 }
466 for (i = 0; i < 24; i += 2)
467 {
468 fprintf(out, "CONECT%5d%5d\n", a0 + rectedge[i], a0 + rectedge[i
469 + 1]);
470 }
471 }
472}
473
474void calc_rotmatrix(rvec principal_axis, rvec targetvec, matrix rotmatrix)
475{
476 rvec rotvec;
477 real ux, uy, uz, costheta, sintheta;
478
479 costheta = cos_angle(principal_axis, targetvec);
480 sintheta = sqrt(1.0-costheta*costheta); /* sign is always positive since 0<theta<pi */
481
482 /* Determine rotation from cross product with target vector */
483 cprod(principal_axis, targetvec, rotvec);
484 unitv(rotvec, rotvec);
485 printf("Aligning %g %g %g to %g %g %g : xprod %g %g %g\n",
486 principal_axis[XX0], principal_axis[YY1], principal_axis[ZZ2], targetvec[XX0], targetvec[YY1], targetvec[ZZ2],
487 rotvec[XX0], rotvec[YY1], rotvec[ZZ2]);
488
489 ux = rotvec[XX0];
490 uy = rotvec[YY1];
491 uz = rotvec[ZZ2];
492 rotmatrix[0][0] = ux*ux + (1.0-ux*ux)*costheta;
493 rotmatrix[0][1] = ux*uy*(1-costheta)-uz*sintheta;
494 rotmatrix[0][2] = ux*uz*(1-costheta)+uy*sintheta;
495 rotmatrix[1][0] = ux*uy*(1-costheta)+uz*sintheta;
496 rotmatrix[1][1] = uy*uy + (1.0-uy*uy)*costheta;
497 rotmatrix[1][2] = uy*uz*(1-costheta)-ux*sintheta;
498 rotmatrix[2][0] = ux*uz*(1-costheta)-uy*sintheta;
499 rotmatrix[2][1] = uy*uz*(1-costheta)+ux*sintheta;
500 rotmatrix[2][2] = uz*uz + (1.0-uz*uz)*costheta;
501
502 printf("Rotation matrix: \n%g %g %g\n%g %g %g\n%g %g %g\n",
503 rotmatrix[0][0], rotmatrix[0][1], rotmatrix[0][2],
504 rotmatrix[1][0], rotmatrix[1][1], rotmatrix[1][2],
505 rotmatrix[2][0], rotmatrix[2][1], rotmatrix[2][2]);
506}
507
508static void renum_resnr(t_atoms *atoms, int isize, const int *index,
509 int resnr_start)
510{
511 int i, resind_prev, resind;
512
513 resind_prev = -1;
514 for (i = 0; i < isize; i++)
515 {
516 resind = atoms->atom[index == NULL((void*)0) ? i : index[i]].resind;
517 if (resind != resind_prev)
518 {
519 atoms->resinfo[resind].nr = resnr_start;
520 resnr_start++;
521 }
522 resind_prev = resind;
523 }
524}
525
526int gmx_editconf(int argc, char *argv[])
527{
528 const char *desc[] =
529 {
530 "[THISMODULE] converts generic structure format to [TT].gro[tt], [TT].g96[tt]",
531 "or [TT].pdb[tt].",
532 "[PAR]",
533 "The box can be modified with options [TT]-box[tt], [TT]-d[tt] and",
534 "[TT]-angles[tt]. Both [TT]-box[tt] and [TT]-d[tt]",
535 "will center the system in the box, unless [TT]-noc[tt] is used.",
536 "[PAR]",
537 "Option [TT]-bt[tt] determines the box type: [TT]triclinic[tt] is a",
538 "triclinic box, [TT]cubic[tt] is a rectangular box with all sides equal",
539 "[TT]dodecahedron[tt] represents a rhombic dodecahedron and",
540 "[TT]octahedron[tt] is a truncated octahedron.",
541 "The last two are special cases of a triclinic box.",
542 "The length of the three box vectors of the truncated octahedron is the",
543 "shortest distance between two opposite hexagons.",
544 "Relative to a cubic box with some periodic image distance, the volume of a ",
545 "dodecahedron with this same periodic distance is 0.71 times that of the cube, ",
546 "and that of a truncated octahedron is 0.77 times.",
547 "[PAR]",
548 "Option [TT]-box[tt] requires only",
549 "one value for a cubic, rhombic dodecahedral, or truncated octahedral box.",
550 "[PAR]",
551 "With [TT]-d[tt] and a [TT]triclinic[tt] box the size of the system in the [IT]x[it]-, [IT]y[it]-,",
552 "and [IT]z[it]-directions is used. With [TT]-d[tt] and [TT]cubic[tt],",
553 "[TT]dodecahedron[tt] or [TT]octahedron[tt] boxes, the dimensions are set",
554 "to the diameter of the system (largest distance between atoms) plus twice",
555 "the specified distance.",
556 "[PAR]",
557 "Option [TT]-angles[tt] is only meaningful with option [TT]-box[tt] and",
558 "a triclinic box and cannot be used with option [TT]-d[tt].",
559 "[PAR]",
560 "When [TT]-n[tt] or [TT]-ndef[tt] is set, a group",
561 "can be selected for calculating the size and the geometric center,",
562 "otherwise the whole system is used.",
563 "[PAR]",
564 "[TT]-rotate[tt] rotates the coordinates and velocities.",
565 "[PAR]",
566 "[TT]-princ[tt] aligns the principal axes of the system along the",
567 "coordinate axes, with the longest axis aligned with the [IT]x[it]-axis. ",
568 "This may allow you to decrease the box volume,",
569 "but beware that molecules can rotate significantly in a nanosecond.",
570 "[PAR]",
571 "Scaling is applied before any of the other operations are",
572 "performed. Boxes and coordinates can be scaled to give a certain density (option",
573 "[TT]-density[tt]). Note that this may be inaccurate in case a [TT].gro[tt]",
574 "file is given as input. A special feature of the scaling option is that when the",
575 "factor -1 is given in one dimension, one obtains a mirror image,",
576 "mirrored in one of the planes. When one uses -1 in three dimensions, ",
577 "a point-mirror image is obtained.[PAR]",
578 "Groups are selected after all operations have been applied.[PAR]",
579 "Periodicity can be removed in a crude manner.",
580 "It is important that the box vectors at the bottom of your input file",
581 "are correct when the periodicity is to be removed.",
582 "[PAR]",
583 "When writing [TT].pdb[tt] files, B-factors can be",
584 "added with the [TT]-bf[tt] option. B-factors are read",
585 "from a file with with following format: first line states number of",
586 "entries in the file, next lines state an index",
587 "followed by a B-factor. The B-factors will be attached per residue",
588 "unless an index is larger than the number of residues or unless the",
589 "[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
590 "be added instead of B-factors. [TT]-legend[tt] will produce",
591 "a row of CA atoms with B-factors ranging from the minimum to the",
592 "maximum value found, effectively making a legend for viewing.",
593 "[PAR]",
594 "With the option [TT]-mead[tt] a special [TT].pdb[tt] ([TT].pqr[tt])",
595 "file for the MEAD electrostatics",
596 "program (Poisson-Boltzmann solver) can be made. A further prerequisite",
597 "is that the input file is a run input file.",
598 "The B-factor field is then filled with the Van der Waals radius",
599 "of the atoms while the occupancy field will hold the charge.",
600 "[PAR]",
601 "The option [TT]-grasp[tt] is similar, but it puts the charges in the B-factor",
602 "and the radius in the occupancy.",
603 "[PAR]",
604 "Option [TT]-align[tt] allows alignment",
605 "of the principal axis of a specified group against the given vector, ",
606 "with an optional center of rotation specified by [TT]-aligncenter[tt].",
607 "[PAR]",
608 "Finally, with option [TT]-label[tt], [TT]editconf[tt] can add a chain identifier",
609 "to a [TT].pdb[tt] file, which can be useful for analysis with e.g. Rasmol.",
610 "[PAR]",
611 "To convert a truncated octrahedron file produced by a package which uses",
612 "a cubic box with the corners cut off (such as GROMOS), use:[BR]",
613 "[TT]gmx editconf -f in -rotate 0 45 35.264 -bt o -box veclen -o out[tt][BR]",
614 "where [TT]veclen[tt] is the size of the cubic box times [SQRT]3[sqrt]/2."
615 };
616 const char *bugs[] =
617 {
618 "For complex molecules, the periodicity removal routine may break down, "
619 "in that case you can use [gmx-trjconv]."
620 };
621 static real dist = 0.0, rbox = 0.0, to_diam = 0.0;
622 static gmx_bool bNDEF = FALSE0, bRMPBC = FALSE0, bCenter = FALSE0, bReadVDW =
623 FALSE0, bCONECT = FALSE0;
624 static gmx_bool peratom = FALSE0, bLegend = FALSE0, bOrient = FALSE0, bMead =
625 FALSE0, bGrasp = FALSE0, bSig56 = FALSE0;
626 static rvec scale =
627 { 1, 1, 1 }, newbox =
628 { 0, 0, 0 }, newang =
629 { 90, 90, 90 };
630 static real rho = 1000.0, rvdw = 0.12;
631 static rvec center =
632 { 0, 0, 0 }, translation =
633 { 0, 0, 0 }, rotangles =
634 { 0, 0, 0 }, aligncenter =
635 { 0, 0, 0 }, targetvec =
636 { 0, 0, 0 };
637 static const char *btype[] =
638 { NULL((void*)0), "triclinic", "cubic", "dodecahedron", "octahedron", NULL((void*)0) },
639 *label = "A";
640 static rvec visbox =
641 { 0, 0, 0 };
642 static int resnr_start = -1;
643 t_pargs
644 pa[] =
645 {
646 { "-ndef", FALSE0, etBOOL,
647 { &bNDEF }, "Choose output from default index groups" },
648 { "-visbox", FALSE0, etRVEC,
649 { visbox },
650 "HIDDENVisualize a grid of boxes, -1 visualizes the 14 box images" },
651 { "-bt", FALSE0, etENUM,
652 { btype }, "Box type for [TT]-box[tt] and [TT]-d[tt]" },
653 { "-box", FALSE0, etRVEC,
654 { newbox }, "Box vector lengths (a,b,c)" },
655 { "-angles", FALSE0, etRVEC,
656 { newang }, "Angles between the box vectors (bc,ac,ab)" },
657 { "-d", FALSE0, etREAL,
658 { &dist }, "Distance between the solute and the box" },
659 { "-c", FALSE0, etBOOL,
660 { &bCenter },
661 "Center molecule in box (implied by [TT]-box[tt] and [TT]-d[tt])" },
662 { "-center", FALSE0, etRVEC,
663 { center }, "Coordinates of geometrical center" },
664 { "-aligncenter", FALSE0, etRVEC,
665 { aligncenter }, "Center of rotation for alignment" },
666 { "-align", FALSE0, etRVEC,
667 { targetvec },
668 "Align to target vector" },
669 { "-translate", FALSE0, etRVEC,
670 { translation }, "Translation" },
671 { "-rotate", FALSE0, etRVEC,
672 { rotangles },
673 "Rotation around the X, Y and Z axes in degrees" },
674 { "-princ", FALSE0, etBOOL,
675 { &bOrient },
676 "Orient molecule(s) along their principal axes" },
677 { "-scale", FALSE0, etRVEC,
678 { scale }, "Scaling factor" },
679 { "-density", FALSE0, etREAL,
680 { &rho },
681 "Density (g/L) of the output box achieved by scaling" },
682 { "-pbc", FALSE0, etBOOL,
683 { &bRMPBC },
684 "Remove the periodicity (make molecule whole again)" },
685 { "-resnr", FALSE0, etINT,
686 { &resnr_start },
687 " Renumber residues starting from resnr" },
688 { "-grasp", FALSE0, etBOOL,
689 { &bGrasp },
690 "Store the charge of the atom in the B-factor field and the radius of the atom in the occupancy field" },
691 {
692 "-rvdw", FALSE0, etREAL,
693 { &rvdw },
694 "Default Van der Waals radius (in nm) if one can not be found in the database or if no parameters are present in the topology file"
695 },
696 { "-sig56", FALSE0, etBOOL,
697 { &bSig56 },
698 "Use rmin/2 (minimum in the Van der Waals potential) rather than [GRK]sigma[grk]/2 " },
699 {
700 "-vdwread", FALSE0, etBOOL,
701 { &bReadVDW },
702 "Read the Van der Waals radii from the file [TT]vdwradii.dat[tt] rather than computing the radii based on the force field"
703 },
704 { "-atom", FALSE0, etBOOL,
705 { &peratom }, "Force B-factor attachment per atom" },
706 { "-legend", FALSE0, etBOOL,
707 { &bLegend }, "Make B-factor legend" },
708 { "-label", FALSE0, etSTR,
709 { &label }, "Add chain label for all residues" },
710 {
711 "-conect", FALSE0, etBOOL,
712 { &bCONECT },
713 "Add CONECT records to a [TT].pdb[tt] file when written. Can only be done when a topology is present"
714 }
715 };
716#define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])))
717
718 FILE *out;
719 const char *infile, *outfile;
720 char title[STRLEN4096];
721 int outftp, inftp, natom, i, j, n_bfac, itype, ntype;
722 double *bfac = NULL((void*)0), c6, c12;
723 int *bfac_nr = NULL((void*)0);
724 t_topology *top = NULL((void*)0);
725 t_atoms atoms;
726 char *grpname, *sgrpname, *agrpname;
727 int isize, ssize, tsize, asize;
728 atom_id *index, *sindex, *tindex, *aindex;
729 rvec *x, *v, gc, min, max, size;
730 int ePBC;
731 matrix box, rotmatrix, trans;
732 rvec princd, tmpvec;
733 gmx_bool bIndex, bSetSize, bSetAng, bCubic, bDist, bSetCenter, bAlign;
734 gmx_bool bHaveV, bScale, bRho, bTranslate, bRotate, bCalcGeom, bCalcDiam;
735 real xs, ys, zs, xcent, ycent, zcent, diam = 0, mass = 0, d, vdw;
736 gmx_atomprop_t aps;
737 gmx_conect conect;
738 output_env_t oenv;
739 t_filenm fnm[] =
740 {
741 { efSTX, "-f", NULL((void*)0), ffREAD1<<1 },
742 { efNDX, "-n", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
743 { efSTO, NULL((void*)0), NULL((void*)0), ffOPTWR(1<<2| 1<<3) },
744 { efPQR, "-mead", "mead", ffOPTWR(1<<2| 1<<3) },
745 { efDAT, "-bf", "bfact", ffOPTRD(1<<1 | 1<<3) }
746 };
747#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
748
749 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW(1<<5), NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa,
1
Taking false branch
→
750 asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, asize(bugs)((int)(sizeof(bugs)/sizeof((bugs)[0]))), bugs, &oenv))
751 {
752 return 0;
753 }
754
755 bIndex = opt2bSet("-n", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || bNDEF;
756 bMead = opt2bSet("-mead", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
757 bSetSize = opt2parg_bSet("-box", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
758 bSetAng = opt2parg_bSet("-angles", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
759 bSetCenter = opt2parg_bSet("-center", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
760 bDist = opt2parg_bSet("-d", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
761 bAlign = opt2parg_bSet("-align", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
762 /* Only automatically turn on centering without -noc */
763 if ((bDist || bSetSize || bSetCenter) && !opt2parg_bSet("-c", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa))
764 {
765 bCenter = TRUE1;
766 }
767 bScale = opt2parg_bSet("-scale", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
768 bRho = opt2parg_bSet("-density", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
769 bTranslate = opt2parg_bSet("-translate", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
770 bRotate = opt2parg_bSet("-rotate", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
771 if (bScale && bRho)
772 {
773 fprintf(stderrstderr, "WARNING: setting -density overrides -scale\n");
774 }
775 bScale = bScale || bRho;
776 bCalcGeom = bCenter || bRotate || bOrient || bScale;
777 bCalcDiam = btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o';
2
←
Array access results in a null pointer dereference
778
779 infile = ftp2fn(efSTX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
780 if (bMead)
781 {
782 outfile = ftp2fn(efPQR, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
783 }
784 else
785 {
786 outfile = ftp2fn(efSTO, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
787 }
788 outftp = fn2ftp(outfile);
789 inftp = fn2ftp(infile);
790
791 aps = gmx_atomprop_init();
792
793 if (bMead && bGrasp)
794 {
795 printf("Incompatible options -mead and -grasp. Turning off -grasp\n");
796 bGrasp = FALSE0;
797 }
798 if (bGrasp && (outftp != efPDB))
799 {
800 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 800, "Output file should be a .pdb file"
801 " when using the -grasp option\n");
802 }
803 if ((bMead || bGrasp) && !((fn2ftp(infile) == efTPR) ||
804 (fn2ftp(infile) == efTPA) ||
805 (fn2ftp(infile) == efTPB)))
806 {
807 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 807, "Input file should be a .tp[abr] file"
808 " when using the -mead option\n");
809 }
810
811 get_stx_coordnum(infile, &natom);
812 init_t_atoms(&atoms, natom, TRUE1);
813 snew(x, natom)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 813, (natom), sizeof(*(x)));
814 snew(v, natom)(v) = save_calloc("v", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 814, (natom), sizeof(*(v)));
815 read_stx_conf(infile, title, &atoms, x, v, &ePBC, box);
816 if (fn2ftp(infile) == efPDB)
817 {
818 get_pdb_atomnumber(&atoms, aps);
819 }
820 printf("Read %d atoms\n", atoms.nr);
821
822 /* Get the element numbers if available in a pdb file */
823 if (fn2ftp(infile) == efPDB)
824 {
825 get_pdb_atomnumber(&atoms, aps);
826 }
827
828 if (ePBC != epbcNONE)
829 {
830 real vol = det(box);
831 printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
832 vol, 100*((int)(vol*4.5)));
833 }
834
835 if (bMead || bGrasp || bCONECT)
836 {
837 top = read_top(infile, NULL((void*)0));
838 }
839
840 if (bMead || bGrasp)
841 {
842 if (atoms.nr != top->atoms.nr)
843 {
844 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 844, "Atom numbers don't match (%d vs. %d)", atoms.nr, top->atoms.nr);
845 }
846 snew(atoms.pdbinfo, top->atoms.nr)(atoms.pdbinfo) = save_calloc("atoms.pdbinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 846, (top->atoms.nr), sizeof(*(atoms.pdbinfo)));
847 ntype = top->idef.atnr;
848 for (i = 0; (i < atoms.nr); i++)
849 {
850 /* Determine the Van der Waals radius from the force field */
851 if (bReadVDW)
852 {
853 if (!gmx_atomprop_query(aps, epropVDW,
854 *top->atoms.resinfo[top->atoms.atom[i].resind].name,
855 *top->atoms.atomname[i], &vdw))
856 {
857 vdw = rvdw;
858 }
859 }
860 else
861 {
862 itype = top->atoms.atom[i].type;
863 c12 = top->idef.iparams[itype*ntype+itype].lj.c12;
864 c6 = top->idef.iparams[itype*ntype+itype].lj.c6;
865 if ((c6 != 0) && (c12 != 0))
866 {
867 real sig6;
868 if (bSig56)
869 {
870 sig6 = 2*c12/c6;
871 }
872 else
873 {
874 sig6 = c12/c6;
875 }
876 vdw = 0.5*pow(sig6, 1.0/6.0);
877 }
878 else
879 {
880 vdw = rvdw;
881 }
882 }
883 /* Factor of 10 for nm -> Angstroms */
884 vdw *= 10;
885
886 if (bMead)
887 {
888 atoms.pdbinfo[i].occup = top->atoms.atom[i].q;
889 atoms.pdbinfo[i].bfac = vdw;
890 }
891 else
892 {
893 atoms.pdbinfo[i].occup = vdw;
894 atoms.pdbinfo[i].bfac = top->atoms.atom[i].q;
895 }
896 }
897 }
898 bHaveV = FALSE0;
899 for (i = 0; (i < natom) && !bHaveV; i++)
900 {
901 for (j = 0; (j < DIM3) && !bHaveV; j++)
902 {
903 bHaveV = bHaveV || (v[i][j] != 0);
904 }
905 }
906 printf("%selocities found\n", bHaveV ? "V" : "No v");
907
908 if (visbox[0] > 0)
909 {
910 if (bIndex)
911 {
912 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 912, "Sorry, can not visualize box with index groups");
913 }
914 if (outftp != efPDB)
915 {
916 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 916, "Sorry, can only visualize box with a pdb file");
917 }
918 }
919 else if (visbox[0] == -1)
920 {
921 visualize_images("images.pdb", ePBC, box);
922 }
923
924 /* remove pbc */
925 if (bRMPBC)
926 {
927 rm_gropbc(&atoms, x, box);
928 }
929
930 if (bCalcGeom)
931 {
932 if (bIndex)
933 {
934 fprintf(stderrstderr, "\nSelect a group for determining the system size:\n");
935 get_index(&atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
936 1, &ssize, &sindex, &sgrpname);
937 }
938 else
939 {
940 ssize = atoms.nr;
941 sindex = NULL((void*)0);
942 }
943 diam = calc_geom(ssize, sindex, x, gc, min, max, bCalcDiam);
944 rvec_sub(max, min, size);
945 printf(" system size :%7.3f%7.3f%7.3f (nm)\n",
946 size[XX0], size[YY1], size[ZZ2]);
947 if (bCalcDiam)
948 {
949 printf(" diameter :%7.3f (nm)\n", diam);
950 }
951 printf(" center :%7.3f%7.3f%7.3f (nm)\n", gc[XX0], gc[YY1], gc[ZZ2]);
952 printf(" box vectors :%7.3f%7.3f%7.3f (nm)\n",
953 norm(box[XX0]), norm(box[YY1]), norm(box[ZZ2]));
954 printf(" box angles :%7.2f%7.2f%7.2f (degrees)\n",
955 norm2(box[ZZ2]) == 0 ? 0 :
956 RAD2DEG(180.0/3.14159265358979323846)*acos(cos_angle_no_table(box[YY1], box[ZZ2])),
957 norm2(box[ZZ2]) == 0 ? 0 :
958 RAD2DEG(180.0/3.14159265358979323846)*acos(cos_angle_no_table(box[XX0], box[ZZ2])),
959 norm2(box[YY1]) == 0 ? 0 :
960 RAD2DEG(180.0/3.14159265358979323846)*acos(cos_angle_no_table(box[XX0], box[YY1])));
961 printf(" box volume :%7.2f (nm^3)\n", det(box));
962 }
963
964 if (bRho || bOrient || bAlign)
965 {
966 mass = calc_mass(&atoms, !fn2bTPX(infile), aps);
967 }
968
969 if (bOrient)
970 {
971 atom_id *index;
972 char *grpnames;
973
974 /* Get a group for principal component analysis */
975 fprintf(stderrstderr, "\nSelect group for the determining the orientation\n");
976 get_index(&atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &isize, &index, &grpnames);
977
978 /* Orient the principal axes along the coordinate axes */
979 orient_princ(&atoms, isize, index, natom, x, bHaveV ? v : NULL((void*)0), NULL((void*)0));
980 sfree(index)save_free("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 980, (index));
981 sfree(grpnames)save_free("grpnames", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 981, (grpnames));
982 }
983
984 if (bScale)
985 {
986 /* scale coordinates and box */
987 if (bRho)
988 {
989 /* Compute scaling constant */
990 real vol, dens;
991
992 vol = det(box);
993 dens = (mass*AMU(1.6605402e-27))/(vol*NANO(1e-9)*NANO(1e-9)*NANO(1e-9));
994 fprintf(stderrstderr, "Volume of input %g (nm^3)\n", vol);
995 fprintf(stderrstderr, "Mass of input %g (a.m.u.)\n", mass);
996 fprintf(stderrstderr, "Density of input %g (g/l)\n", dens);
997 if (vol == 0 || mass == 0)
998 {
999 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 999, "Cannot scale density with "
1000 "zero mass (%g) or volume (%g)\n", mass, vol);
1001 }
1002
1003 scale[XX0] = scale[YY1] = scale[ZZ2] = pow(dens/rho, 1.0/3.0);
1004 fprintf(stderrstderr, "Scaling all box vectors by %g\n", scale[XX0]);
1005 }
1006 scale_conf(atoms.nr, x, box, scale);
1007 }
1008
1009 if (bAlign)
1010 {
1011 if (bIndex)
1012 {
1013 fprintf(stderrstderr, "\nSelect a group that you want to align:\n");
1014 get_index(&atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1015 1, &asize, &aindex, &agrpname);
1016 }
1017 else
1018 {
1019 asize = atoms.nr;
1020 snew(aindex, asize)(aindex) = save_calloc("aindex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 1020, (asize), sizeof(*(aindex)));
1021 for (i = 0; i < asize; i++)
1022 {
1023 aindex[i] = i;
1024 }
1025 }
1026 printf("Aligning %d atoms (out of %d) to %g %g %g, center of rotation %g %g %g\n", asize, natom,
1027 targetvec[XX0], targetvec[YY1], targetvec[ZZ2],
1028 aligncenter[XX0], aligncenter[YY1], aligncenter[ZZ2]);
1029 /*subtract out pivot point*/
1030 for (i = 0; i < asize; i++)
1031 {
1032 rvec_dec(x[aindex[i]], aligncenter);
1033 }
1034 /*now determine transform and rotate*/
1035 /*will this work?*/
1036 principal_comp(asize, aindex, atoms.atom, x, trans, princd);
1037
1038 unitv(targetvec, targetvec);
1039 printf("Using %g %g %g as principal axis\n", trans[0][2], trans[1][2], trans[2][2]);
1040 tmpvec[XX0] = trans[0][2]; tmpvec[YY1] = trans[1][2]; tmpvec[ZZ2] = trans[2][2];
1041 calc_rotmatrix(tmpvec, targetvec, rotmatrix);
1042 /* rotmatrix finished */
1043
1044 for (i = 0; i < asize; ++i)
1045 {
1046 mvmul(rotmatrix, x[aindex[i]], tmpvec);
1047 copy_rvec(tmpvec, x[aindex[i]]);
1048 }
1049
1050 /*add pivot point back*/
1051 for (i = 0; i < asize; i++)
1052 {
1053 rvec_inc(x[aindex[i]], aligncenter);
1054 }
1055 if (!bIndex)
1056 {
1057 sfree(aindex)save_free("aindex", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 1057, (aindex));
1058 }
1059 }
1060
1061 if (bTranslate)
1062 {
1063 if (bIndex)
1064 {
1065 fprintf(stderrstderr, "\nSelect a group that you want to translate:\n");
1066 get_index(&atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1067 1, &ssize, &sindex, &sgrpname);
1068 }
1069 else
1070 {
1071 ssize = atoms.nr;
1072 sindex = NULL((void*)0);
1073 }
1074 printf("Translating %d atoms (out of %d) by %g %g %g nm\n", ssize, natom,
1075 translation[XX0], translation[YY1], translation[ZZ2]);
1076 if (sindex)
1077 {
1078 for (i = 0; i < ssize; i++)
1079 {
1080 rvec_inc(x[sindex[i]], translation);
1081 }
1082 }
1083 else
1084 {
1085 for (i = 0; i < natom; i++)
1086 {
1087 rvec_inc(x[i], translation);
1088 }
1089 }
1090 }
1091 if (bRotate)
1092 {
1093 /* Rotate */
1094 printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n", rotangles[XX0], rotangles[YY1], rotangles[ZZ2]);
1095 for (i = 0; i < DIM3; i++)
1096 {
1097 rotangles[i] *= DEG2RAD(3.14159265358979323846/180.0);
1098 }
1099 rotate_conf(natom, x, v, rotangles[XX0], rotangles[YY1], rotangles[ZZ2]);
1100 }
1101
1102 if (bCalcGeom)
1103 {
1104 /* recalc geometrical center and max and min coordinates and size */
1105 calc_geom(ssize, sindex, x, gc, min, max, FALSE0);
1106 rvec_sub(max, min, size);
1107 if (bScale || bOrient || bRotate)
1108 {
1109 printf("new system size : %6.3f %6.3f %6.3f\n",
1110 size[XX0], size[YY1], size[ZZ2]);
1111 }
1112 }
1113
1114 if (bSetSize || bDist || (btype[0][0] == 't' && bSetAng))
1115 {
1116 ePBC = epbcXYZ;
1117 if (!(bSetSize || bDist))
1118 {
1119 for (i = 0; i < DIM3; i++)
1120 {
1121 newbox[i] = norm(box[i]);
1122 }
1123 }
1124 clear_mat(box);
1125 /* calculate new boxsize */
1126 switch (btype[0][0])
1127 {
1128 case 't':
1129 if (bDist)
1130 {
1131 for (i = 0; i < DIM3; i++)
1132 {
1133 newbox[i] = size[i]+2*dist;
1134 }
1135 }
1136 if (!bSetAng)
1137 {
1138 box[XX0][XX0] = newbox[XX0];
1139 box[YY1][YY1] = newbox[YY1];
1140 box[ZZ2][ZZ2] = newbox[ZZ2];
1141 }
1142 else
1143 {
1144 matrix_convert(box, newbox, newang);
1145 }
1146 break;
1147 case 'c':
1148 case 'd':
1149 case 'o':
1150 if (bSetSize)
1151 {
1152 d = newbox[0];
1153 }
1154 else
1155 {
1156 d = diam+2*dist;
1157 }
1158 if (btype[0][0] == 'c')
1159 {
1160 for (i = 0; i < DIM3; i++)
1161 {
1162 box[i][i] = d;
1163 }
1164 }
1165 else if (btype[0][0] == 'd')
1166 {
1167 box[XX0][XX0] = d;
1168 box[YY1][YY1] = d;
1169 box[ZZ2][XX0] = d/2;
1170 box[ZZ2][YY1] = d/2;
1171 box[ZZ2][ZZ2] = d*sqrt(2)/2;
1172 }
1173 else
1174 {
1175 box[XX0][XX0] = d;
1176 box[YY1][XX0] = d/3;
1177 box[YY1][YY1] = d*sqrt(2)*2/3;
1178 box[ZZ2][XX0] = -d/3;
1179 box[ZZ2][YY1] = d*sqrt(2)/3;
1180 box[ZZ2][ZZ2] = d*sqrt(6)/3;
1181 }
1182 break;
1183 }
1184 }
1185
1186 /* calculate new coords for geometrical center */
1187 if (!bSetCenter)
1188 {
1189 calc_box_center(ecenterDEF, box, center);
1190 }
1191
1192 /* center molecule on 'center' */
1193 if (bCenter)
1194 {
1195 center_conf(natom, x, center, gc);
1196 }
1197
1198 /* print some */
1199 if (bCalcGeom)
1200 {
1201 calc_geom(ssize, sindex, x, gc, min, max, FALSE0);
1202 printf("new center :%7.3f%7.3f%7.3f (nm)\n", gc[XX0], gc[YY1], gc[ZZ2]);
1203 }
1204 if (bOrient || bScale || bDist || bSetSize)
1205 {
1206 printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n",
1207 norm(box[XX0]), norm(box[YY1]), norm(box[ZZ2]));
1208 printf("new box angles :%7.2f%7.2f%7.2f (degrees)\n",
1209 norm2(box[ZZ2]) == 0 ? 0 :
1210 RAD2DEG(180.0/3.14159265358979323846)*acos(cos_angle_no_table(box[YY1], box[ZZ2])),
1211 norm2(box[ZZ2]) == 0 ? 0 :
1212 RAD2DEG(180.0/3.14159265358979323846)*acos(cos_angle_no_table(box[XX0], box[ZZ2])),
1213 norm2(box[YY1]) == 0 ? 0 :
1214 RAD2DEG(180.0/3.14159265358979323846)*acos(cos_angle_no_table(box[XX0], box[YY1])));
1215 printf("new box volume :%7.2f (nm^3)\n", det(box));
1216 }
1217
1218 if (check_box(epbcXYZ, box))
1219 {
1220 printf("\nWARNING: %s\n"
1221 "See the GROMACS manual for a description of the requirements that\n"
1222 "must be satisfied by descriptions of simulation cells.\n",
1223 check_box(epbcXYZ, box));
1224 }
1225
1226 if (bDist && btype[0][0] == 't')
1227 {
1228 if (TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0))
1229 {
1230 printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
1231 "distance from the solute to a box surface along the corresponding normal\n"
1232 "vector might be somewhat smaller than your specified value %f.\n"
1233 "You can check the actual value with g_mindist -pi\n", dist);
1234 }
1235 else if (!opt2parg_bSet("-bt", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa))
1236 {
1237 printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
1238 "If the molecule rotates the actual distance will be smaller. You might want\n"
1239 "to use a cubic box instead, or why not try a dodecahedron today?\n");
1240 }
1241 }
1242 if (bCONECT && (outftp == efPDB) && (inftp == efTPR))
1243 {
1244 conect = gmx_conect_generate(top);
1245 }
1246 else
1247 {
1248 conect = NULL((void*)0);
1249 }
1250
1251 if (bIndex)
1252 {
1253 fprintf(stderrstderr, "\nSelect a group for output:\n");
1254 get_index(&atoms, opt2fn_null("-n", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1255 1, &isize, &index, &grpname);
1256
1257 if (resnr_start >= 0)
1258 {
1259 renum_resnr(&atoms, isize, index, resnr_start);
1260 }
1261
1262 if (opt2parg_bSet("-label", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa))
1263 {
1264 for (i = 0; (i < atoms.nr); i++)
1265 {
1266 atoms.resinfo[atoms.atom[i].resind].chainid = label[0];
1267 }
1268 }
1269
1270 if (opt2bSet("-bf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || bLegend)
1271 {
1272 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_editconf.c"
, 1272, "Sorry, cannot do bfactors with an index group.");
1273 }
1274
1275 if (outftp == efPDB)
1276 {
1277 out = gmx_ffopen(outfile, "w");
1278 write_pdbfile_indexed(out, title, &atoms, x, ePBC, box, ' ', 1, isize, index, conect, TRUE1);
1279 gmx_ffclose(out);
1280 }
1281 else
1282 {
1283 write_sto_conf_indexed(outfile, title, &atoms, x, bHaveV ? v : NULL((void*)0), ePBC, box, isize, index);
1284 }
1285 }
1286 else
1287 {
1288 if (resnr_start >= 0)
1289 {
1290 renum_resnr(&atoms, atoms.nr, NULL((void*)0), resnr_start);
1291 }
1292
1293 if ((outftp == efPDB) || (outftp == efPQR))
1294 {
1295 out = gmx_ffopen(outfile, "w");
1296 if (bMead)
1297 {
1298 set_pdb_wide_format(TRUE1);
1299 fprintf(out, "REMARK "
1300 "The B-factors in this file hold atomic radii\n");
1301 fprintf(out, "REMARK "
1302 "The occupancy in this file hold atomic charges\n");
1303 }
1304 else if (bGrasp)
1305 {
1306 fprintf(out, "GRASP PDB FILE\nFORMAT NUMBER=1\n");
1307 fprintf(out, "REMARK "
1308 "The B-factors in this file hold atomic charges\n");
1309 fprintf(out, "REMARK "
1310 "The occupancy in this file hold atomic radii\n");
1311 }
1312 else if (opt2bSet("-bf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
1313 {
1314 read_bfac(opt2fn("-bf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &n_bfac, &bfac, &bfac_nr);
1315 set_pdb_conf_bfac(atoms.nr, atoms.nres, &atoms,
1316 n_bfac, bfac, bfac_nr, peratom);
1317 }
1318 if (opt2parg_bSet("-label", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa))
1319 {
1320 for (i = 0; (i < atoms.nr); i++)
1321 {
1322 atoms.resinfo[atoms.atom[i].resind].chainid = label[0];
1323 }
1324 }
1325 write_pdbfile(out, title, &atoms, x, ePBC, box, ' ', -1, conect, TRUE1);
1326 if (bLegend)
1327 {
1328 pdb_legend(out, atoms.nr, atoms.nres, &atoms, x);
1329 }
1330 if (visbox[0] > 0)
1331 {
1332 visualize_box(out, bLegend ? atoms.nr+12 : atoms.nr,
1333 bLegend ? atoms.nres = 12 : atoms.nres, box, visbox);
1334 }
1335 gmx_ffclose(out);
1336 }
1337 else
1338 {
1339 write_sto_conf(outfile, title, &atoms, x, bHaveV ? v : NULL((void*)0), ePBC, box);
1340 }
1341 }
1342 gmx_atomprop_destroy(aps);
1343
1344 do_view(oenv, outfile, NULL((void*)0));
1345
1346 return 0;
1347}