2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, 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.
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.
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.
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.
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.
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.
44 #include "gromacs/fileio/pdbio.h"
45 #include "gromacs/fileio/confio.h"
47 #include "gromacs/utility/smalloc.h"
49 #include "gromacs/commandline/pargs.h"
50 #include "gromacs/fileio/strdb.h"
52 #include "gromacs/math/vec.h"
54 #include "gromacs/gmxlib/conformation-utilities.h"
55 #include "gromacs/math/units.h"
57 #include "gromacs/fileio/tpxio.h"
58 #include "gromacs/fileio/trxio.h"
66 #include "gromacs/utility/fatalerror.h"
84 real
calc_mass(t_atoms
*atoms
, gmx_bool bGetMass
, gmx_atomprop_t aps
)
90 for (i
= 0; (i
< atoms
->nr
); i
++)
94 gmx_atomprop_query(aps
, epropMass
,
95 *atoms
->resinfo
[atoms
->atom
[i
].resind
].name
,
96 *atoms
->atomname
[i
], &(atoms
->atom
[i
].m
));
98 tmass
+= atoms
->atom
[i
].m
;
104 real
calc_geom(int isize
, atom_id
*index
, rvec
*x
, rvec geom_center
, rvec min
,
105 rvec max
, gmx_bool bDiam
)
111 clear_rvec(geom_center
);
128 for (j
= 0; j
< DIM
; j
++)
130 min
[j
] = max
[j
] = x
[ii
][j
];
132 for (i
= 0; i
< isize
; i
++)
142 rvec_inc(geom_center
, x
[ii
]);
143 for (j
= 0; j
< DIM
; j
++)
145 if (x
[ii
][j
] < min
[j
])
149 if (x
[ii
][j
] > max
[j
])
158 for (j
= i
+ 1; j
< isize
; j
++)
160 d
= distance2(x
[ii
], x
[index
[j
]]);
161 diam2
= max(d
, diam2
);
166 for (j
= i
+ 1; j
< isize
; j
++)
168 d
= distance2(x
[i
], x
[j
]);
169 diam2
= max(d
, diam2
);
174 svmul(1.0 / isize
, geom_center
, geom_center
);
180 void center_conf(int natom
, rvec
*x
, rvec center
, rvec geom_cent
)
185 rvec_sub(center
, geom_cent
, shift
);
187 printf(" shift :%7.3f%7.3f%7.3f (nm)\n", shift
[XX
], shift
[YY
],
190 for (i
= 0; (i
< natom
); i
++)
192 rvec_inc(x
[i
], shift
);
196 void scale_conf(int natom
, rvec x
[], matrix box
, rvec scale
)
200 for (i
= 0; i
< natom
; i
++)
202 for (j
= 0; j
< DIM
; j
++)
207 for (i
= 0; i
< DIM
; i
++)
209 for (j
= 0; j
< DIM
; j
++)
211 box
[i
][j
] *= scale
[j
];
216 void read_bfac(const char *fn
, int *n_bfac
, double **bfac_val
, int **bfac_nr
)
221 *n_bfac
= get_lines(fn
, &bfac_lines
);
222 snew(*bfac_val
, *n_bfac
);
223 snew(*bfac_nr
, *n_bfac
);
224 fprintf(stderr
, "Reading %d B-factors from %s\n", *n_bfac
, fn
);
225 for (i
= 0; (i
< *n_bfac
); i
++)
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]);*/
234 void set_pdb_conf_bfac(int natoms
, int nres
, t_atoms
*atoms
, int n_bfac
,
235 double *bfac
, int *bfac_nr
, gmx_bool peratom
)
238 real bfac_min
, bfac_max
;
244 for (i
= 0; (i
< n_bfac
); i
++)
246 if (bfac_nr
[i
] - 1 >= atoms
->nres
)
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
)
257 if (bfac
[i
] < bfac_min
)
262 while ((bfac_max
> 99.99) || (bfac_min
< -99.99))
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
++)
274 while ((fabs(bfac_max
) < 0.5) && (fabs(bfac_min
) < 0.5))
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
++)
287 for (i
= 0; (i
< natoms
); i
++)
289 atoms
->pdbinfo
[i
].bfac
= 0;
294 fprintf(stderr
, "Will attach %d B-factors to %d residues\n", n_bfac
,
296 for (i
= 0; (i
< n_bfac
); i
++)
299 for (n
= 0; (n
< natoms
); n
++)
301 if (bfac_nr
[i
] == atoms
->resinfo
[atoms
->atom
[n
].resind
].nr
)
303 atoms
->pdbinfo
[n
].bfac
= bfac
[i
];
309 gmx_warning("Residue nr %d not found\n", bfac_nr
[i
]);
315 fprintf(stderr
, "Will attach %d B-factors to %d atoms\n", n_bfac
,
317 for (i
= 0; (i
< n_bfac
); i
++)
319 atoms
->pdbinfo
[bfac_nr
[i
] - 1].bfac
= bfac
[i
];
324 void pdb_legend(FILE *out
, int natoms
, int nres
, t_atoms
*atoms
, rvec x
[])
326 real bfac_min
, bfac_max
, xmin
, ymin
, zmin
;
335 for (i
= 0; (i
< natoms
); i
++)
337 xmin
= min(xmin
, x
[i
][XX
]);
338 ymin
= min(ymin
, x
[i
][YY
]);
339 zmin
= min(zmin
, x
[i
][ZZ
]);
340 bfac_min
= min(bfac_min
, atoms
->pdbinfo
[i
].bfac
);
341 bfac_max
= max(bfac_max
, atoms
->pdbinfo
[i
].bfac
);
343 fprintf(stderr
, "B-factors range from %g to %g\n", bfac_min
, bfac_max
);
344 for (i
= 1; (i
< 12); i
++)
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));
354 void visualize_images(const char *fn
, int ePBC
, matrix box
)
362 init_t_atoms(&atoms
, nat
, FALSE
);
365 /* FIXME: Constness should not be cast away */
367 ala
= (char *) "ALA";
368 for (i
= 0; i
< nat
; i
++)
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
/ NCUCVERT
;
376 calc_triclinic_images(box
, img
+ 1);
378 write_sto_conf(fn
, "Images", &atoms
, img
, NULL
, ePBC
, box
);
380 free_t_atoms(&atoms
, FALSE
);
384 void visualize_box(FILE *out
, int a0
, int r0
, matrix box
, rvec gridsize
)
388 int nx
, ny
, nz
, nbox
, nat
;
392 0, 1, 1, 3, 3, 2, 0, 2, 0, 4, 1, 5, 3, 7, 2, 6, 4, 5, 5, 7, 7, 6, 6,
399 nx
= (int) (gridsize
[XX
] + 0.5);
400 ny
= (int) (gridsize
[YY
] + 0.5);
401 nz
= (int) (gridsize
[ZZ
] + 0.5);
405 nat
= nbox
* NCUCVERT
;
407 calc_compact_unitcell_vertices(ecenterDEF
, box
, vert
);
409 for (z
= 0; z
< nz
; z
++)
411 for (y
= 0; y
< ny
; y
++)
413 for (x
= 0; x
< nx
; x
++)
415 for (i
= 0; i
< DIM
; i
++)
417 shift
[i
] = x
* box
[0][i
] + y
* box
[1][i
] + z
420 for (i
= 0; i
< NCUCVERT
; i
++)
422 rvec_add(vert
[i
], shift
, vert
[j
]);
429 for (i
= 0; i
< nat
; i
++)
431 fprintf(out
, get_pdbformat(), "ATOM", a0
+ i
, "C", "BOX", 'K' + i
432 / NCUCVERT
, r0
+ i
, 10 * vert
[i
][XX
], 10 * vert
[i
][YY
], 10
437 edge
= compact_unitcell_edges();
438 for (j
= 0; j
< nbox
; j
++)
440 for (i
= 0; i
< NCUCEDGE
; i
++)
442 fprintf(out
, "CONECT%5d%5d\n", a0
+ j
* NCUCVERT
+ edge
[2 * i
],
443 a0
+ j
* NCUCVERT
+ edge
[2 * i
+ 1]);
452 for (z
= 0; z
<= 1; z
++)
454 for (y
= 0; y
<= 1; y
++)
456 for (x
= 0; x
<= 1; x
++)
458 fprintf(out
, get_pdbformat(), "ATOM", a0
+ i
, "C", "BOX", 'K' + i
459 / 8, r0
+ i
, x
* 10 * box
[XX
][XX
],
460 y
* 10 * box
[YY
][YY
], z
* 10 * box
[ZZ
][ZZ
]);
466 for (i
= 0; i
< 24; i
+= 2)
468 fprintf(out
, "CONECT%5d%5d\n", a0
+ rectedge
[i
], a0
+ rectedge
[i
474 void calc_rotmatrix(rvec principal_axis
, rvec targetvec
, matrix rotmatrix
)
477 real ux
, uy
, uz
, costheta
, sintheta
;
479 costheta
= cos_angle(principal_axis
, targetvec
);
480 sintheta
= sqrt(1.0-costheta
*costheta
); /* sign is always positive since 0<theta<pi */
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
[XX
], principal_axis
[YY
], principal_axis
[ZZ
], targetvec
[XX
], targetvec
[YY
], targetvec
[ZZ
],
487 rotvec
[XX
], rotvec
[YY
], rotvec
[ZZ
]);
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
;
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]);
508 static void renum_resnr(t_atoms
*atoms
, int isize
, const int *index
,
511 int i
, resind_prev
, resind
;
514 for (i
= 0; i
< isize
; i
++)
516 resind
= atoms
->atom
[index
== NULL
? i
: index
[i
]].resind
;
517 if (resind
!= resind_prev
)
519 atoms
->resinfo
[resind
].nr
= resnr_start
;
522 resind_prev
= resind
;
526 int gmx_editconf(int argc
, char *argv
[])
530 "[THISMODULE] converts generic structure format to [TT].gro[tt], [TT].g96[tt]",
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.",
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.",
548 "Option [TT]-box[tt] requires only",
549 "one value for a cubic, rhombic dodecahedral, or truncated octahedral box.",
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.",
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].",
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.",
564 "[TT]-rotate[tt] rotates the coordinates and velocities.",
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.",
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.",
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.",
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.",
601 "The option [TT]-grasp[tt] is similar, but it puts the charges in the B-factor",
602 "and the radius in the occupancy.",
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].",
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.",
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."
618 "For complex molecules, the periodicity removal routine may break down, "
619 "in that case you can use [gmx-trjconv]."
621 static real dist
= 0.0, rbox
= 0.0, to_diam
= 0.0;
622 static gmx_bool bNDEF
= FALSE
, bRMPBC
= FALSE
, bCenter
= FALSE
, bReadVDW
=
623 FALSE
, bCONECT
= FALSE
;
624 static gmx_bool peratom
= FALSE
, bLegend
= FALSE
, bOrient
= FALSE
, bMead
=
625 FALSE
, bGrasp
= FALSE
, bSig56
= FALSE
;
627 { 1, 1, 1 }, newbox
=
628 { 0, 0, 0 }, newang
=
630 static real rho
= 1000.0, rvdw
= 0.12;
632 { 0, 0, 0 }, translation
=
633 { 0, 0, 0 }, rotangles
=
634 { 0, 0, 0 }, aligncenter
=
635 { 0, 0, 0 }, targetvec
=
637 static const char *btype
[] =
638 { NULL
, "triclinic", "cubic", "dodecahedron", "octahedron", NULL
},
642 static int resnr_start
= -1;
646 { "-ndef", FALSE
, etBOOL
,
647 { &bNDEF
}, "Choose output from default index groups" },
648 { "-visbox", FALSE
, etRVEC
,
650 "HIDDENVisualize a grid of boxes, -1 visualizes the 14 box images" },
651 { "-bt", FALSE
, etENUM
,
652 { btype
}, "Box type for [TT]-box[tt] and [TT]-d[tt]" },
653 { "-box", FALSE
, etRVEC
,
654 { newbox
}, "Box vector lengths (a,b,c)" },
655 { "-angles", FALSE
, etRVEC
,
656 { newang
}, "Angles between the box vectors (bc,ac,ab)" },
657 { "-d", FALSE
, etREAL
,
658 { &dist
}, "Distance between the solute and the box" },
659 { "-c", FALSE
, etBOOL
,
661 "Center molecule in box (implied by [TT]-box[tt] and [TT]-d[tt])" },
662 { "-center", FALSE
, etRVEC
,
663 { center
}, "Coordinates of geometrical center" },
664 { "-aligncenter", FALSE
, etRVEC
,
665 { aligncenter
}, "Center of rotation for alignment" },
666 { "-align", FALSE
, etRVEC
,
668 "Align to target vector" },
669 { "-translate", FALSE
, etRVEC
,
670 { translation
}, "Translation" },
671 { "-rotate", FALSE
, etRVEC
,
673 "Rotation around the X, Y and Z axes in degrees" },
674 { "-princ", FALSE
, etBOOL
,
676 "Orient molecule(s) along their principal axes" },
677 { "-scale", FALSE
, etRVEC
,
678 { scale
}, "Scaling factor" },
679 { "-density", FALSE
, etREAL
,
681 "Density (g/L) of the output box achieved by scaling" },
682 { "-pbc", FALSE
, etBOOL
,
684 "Remove the periodicity (make molecule whole again)" },
685 { "-resnr", FALSE
, etINT
,
687 " Renumber residues starting from resnr" },
688 { "-grasp", FALSE
, etBOOL
,
690 "Store the charge of the atom in the B-factor field and the radius of the atom in the occupancy field" },
692 "-rvdw", FALSE
, etREAL
,
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"
696 { "-sig56", FALSE
, etBOOL
,
698 "Use rmin/2 (minimum in the Van der Waals potential) rather than [GRK]sigma[grk]/2 " },
700 "-vdwread", FALSE
, etBOOL
,
702 "Read the Van der Waals radii from the file [TT]vdwradii.dat[tt] rather than computing the radii based on the force field"
704 { "-atom", FALSE
, etBOOL
,
705 { &peratom
}, "Force B-factor attachment per atom" },
706 { "-legend", FALSE
, etBOOL
,
707 { &bLegend
}, "Make B-factor legend" },
708 { "-label", FALSE
, etSTR
,
709 { &label
}, "Add chain label for all residues" },
711 "-conect", FALSE
, etBOOL
,
713 "Add CONECT records to a [TT].pdb[tt] file when written. Can only be done when a topology is present"
716 #define NPA asize(pa)
719 const char *infile
, *outfile
;
721 int outftp
, inftp
, natom
, i
, j
, n_bfac
, itype
, ntype
;
722 double *bfac
= NULL
, c6
, c12
;
724 t_topology
*top
= NULL
;
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
;
731 matrix box
, rotmatrix
, trans
;
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
;
741 { efSTX
, "-f", NULL
, ffREAD
},
742 { efNDX
, "-n", NULL
, ffOPTRD
},
743 { efSTO
, NULL
, NULL
, ffOPTWR
},
744 { efPQR
, "-mead", "mead", ffOPTWR
},
745 { efDAT
, "-bf", "bfact", ffOPTRD
}
747 #define NFILE asize(fnm)
749 if (!parse_common_args(&argc
, argv
, PCA_CAN_VIEW
, NFILE
, fnm
, NPA
, pa
,
750 asize(desc
), desc
, asize(bugs
), bugs
, &oenv
))
755 bIndex
= opt2bSet("-n", NFILE
, fnm
) || bNDEF
;
756 bMead
= opt2bSet("-mead", NFILE
, fnm
);
757 bSetSize
= opt2parg_bSet("-box", NPA
, pa
);
758 bSetAng
= opt2parg_bSet("-angles", NPA
, pa
);
759 bSetCenter
= opt2parg_bSet("-center", NPA
, pa
);
760 bDist
= opt2parg_bSet("-d", NPA
, pa
);
761 bAlign
= opt2parg_bSet("-align", NPA
, pa
);
762 /* Only automatically turn on centering without -noc */
763 if ((bDist
|| bSetSize
|| bSetCenter
) && !opt2parg_bSet("-c", NPA
, pa
))
767 bScale
= opt2parg_bSet("-scale", NPA
, pa
);
768 bRho
= opt2parg_bSet("-density", NPA
, pa
);
769 bTranslate
= opt2parg_bSet("-translate", NPA
, pa
);
770 bRotate
= opt2parg_bSet("-rotate", NPA
, pa
);
773 fprintf(stderr
, "WARNING: setting -density overrides -scale\n");
775 bScale
= bScale
|| bRho
;
776 bCalcGeom
= bCenter
|| bRotate
|| bOrient
|| bScale
;
777 bCalcDiam
= btype
[0][0] == 'c' || btype
[0][0] == 'd' || btype
[0][0] == 'o';
779 infile
= ftp2fn(efSTX
, NFILE
, fnm
);
782 outfile
= ftp2fn(efPQR
, NFILE
, fnm
);
786 outfile
= ftp2fn(efSTO
, NFILE
, fnm
);
788 outftp
= fn2ftp(outfile
);
789 inftp
= fn2ftp(infile
);
791 aps
= gmx_atomprop_init();
795 printf("Incompatible options -mead and -grasp. Turning off -grasp\n");
798 if (bGrasp
&& (outftp
!= efPDB
))
800 gmx_fatal(FARGS
, "Output file should be a .pdb file"
801 " when using the -grasp option\n");
803 if ((bMead
|| bGrasp
) && !((fn2ftp(infile
) == efTPR
) ||
804 (fn2ftp(infile
) == efTPA
) ||
805 (fn2ftp(infile
) == efTPB
)))
807 gmx_fatal(FARGS
, "Input file should be a .tp[abr] file"
808 " when using the -mead option\n");
811 get_stx_coordnum(infile
, &natom
);
812 init_t_atoms(&atoms
, natom
, TRUE
);
815 read_stx_conf(infile
, title
, &atoms
, x
, v
, &ePBC
, box
);
816 if (fn2ftp(infile
) == efPDB
)
818 get_pdb_atomnumber(&atoms
, aps
);
820 printf("Read %d atoms\n", atoms
.nr
);
822 /* Get the element numbers if available in a pdb file */
823 if (fn2ftp(infile
) == efPDB
)
825 get_pdb_atomnumber(&atoms
, aps
);
828 if (ePBC
!= epbcNONE
)
831 printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
832 vol
, 100*((int)(vol
*4.5)));
835 if (bMead
|| bGrasp
|| bCONECT
)
837 top
= read_top(infile
, NULL
);
842 if (atoms
.nr
!= top
->atoms
.nr
)
844 gmx_fatal(FARGS
, "Atom numbers don't match (%d vs. %d)", atoms
.nr
, top
->atoms
.nr
);
846 snew(atoms
.pdbinfo
, top
->atoms
.nr
);
847 ntype
= top
->idef
.atnr
;
848 for (i
= 0; (i
< atoms
.nr
); i
++)
850 /* Determine the Van der Waals radius from the force field */
853 if (!gmx_atomprop_query(aps
, epropVDW
,
854 *top
->atoms
.resinfo
[top
->atoms
.atom
[i
].resind
].name
,
855 *top
->atoms
.atomname
[i
], &vdw
))
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))
876 vdw
= 0.5*pow(sig6
, 1.0/6.0);
883 /* Factor of 10 for nm -> Angstroms */
888 atoms
.pdbinfo
[i
].occup
= top
->atoms
.atom
[i
].q
;
889 atoms
.pdbinfo
[i
].bfac
= vdw
;
893 atoms
.pdbinfo
[i
].occup
= vdw
;
894 atoms
.pdbinfo
[i
].bfac
= top
->atoms
.atom
[i
].q
;
899 for (i
= 0; (i
< natom
) && !bHaveV
; i
++)
901 for (j
= 0; (j
< DIM
) && !bHaveV
; j
++)
903 bHaveV
= bHaveV
|| (v
[i
][j
] != 0);
906 printf("%selocities found\n", bHaveV
? "V" : "No v");
912 gmx_fatal(FARGS
, "Sorry, can not visualize box with index groups");
916 gmx_fatal(FARGS
, "Sorry, can only visualize box with a pdb file");
919 else if (visbox
[0] == -1)
921 visualize_images("images.pdb", ePBC
, box
);
927 rm_gropbc(&atoms
, x
, box
);
934 fprintf(stderr
, "\nSelect a group for determining the system size:\n");
935 get_index(&atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
),
936 1, &ssize
, &sindex
, &sgrpname
);
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
[XX
], size
[YY
], size
[ZZ
]);
949 printf(" diameter :%7.3f (nm)\n", diam
);
951 printf(" center :%7.3f%7.3f%7.3f (nm)\n", gc
[XX
], gc
[YY
], gc
[ZZ
]);
952 printf(" box vectors :%7.3f%7.3f%7.3f (nm)\n",
953 norm(box
[XX
]), norm(box
[YY
]), norm(box
[ZZ
]));
954 printf(" box angles :%7.2f%7.2f%7.2f (degrees)\n",
955 norm2(box
[ZZ
]) == 0 ? 0 :
956 RAD2DEG
*acos(cos_angle_no_table(box
[YY
], box
[ZZ
])),
957 norm2(box
[ZZ
]) == 0 ? 0 :
958 RAD2DEG
*acos(cos_angle_no_table(box
[XX
], box
[ZZ
])),
959 norm2(box
[YY
]) == 0 ? 0 :
960 RAD2DEG
*acos(cos_angle_no_table(box
[XX
], box
[YY
])));
961 printf(" box volume :%7.2f (nm^3)\n", det(box
));
964 if (bRho
|| bOrient
|| bAlign
)
966 mass
= calc_mass(&atoms
, !fn2bTPX(infile
), aps
);
974 /* Get a group for principal component analysis */
975 fprintf(stderr
, "\nSelect group for the determining the orientation\n");
976 get_index(&atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &isize
, &index
, &grpnames
);
978 /* Orient the principal axes along the coordinate axes */
979 orient_princ(&atoms
, isize
, index
, natom
, x
, bHaveV
? v
: NULL
, NULL
);
986 /* scale coordinates and box */
989 /* Compute scaling constant */
993 dens
= (mass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
994 fprintf(stderr
, "Volume of input %g (nm^3)\n", vol
);
995 fprintf(stderr
, "Mass of input %g (a.m.u.)\n", mass
);
996 fprintf(stderr
, "Density of input %g (g/l)\n", dens
);
997 if (vol
== 0 || mass
== 0)
999 gmx_fatal(FARGS
, "Cannot scale density with "
1000 "zero mass (%g) or volume (%g)\n", mass
, vol
);
1003 scale
[XX
] = scale
[YY
] = scale
[ZZ
] = pow(dens
/rho
, 1.0/3.0);
1004 fprintf(stderr
, "Scaling all box vectors by %g\n", scale
[XX
]);
1006 scale_conf(atoms
.nr
, x
, box
, scale
);
1013 fprintf(stderr
, "\nSelect a group that you want to align:\n");
1014 get_index(&atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
),
1015 1, &asize
, &aindex
, &agrpname
);
1020 snew(aindex
, asize
);
1021 for (i
= 0; i
< asize
; i
++)
1026 printf("Aligning %d atoms (out of %d) to %g %g %g, center of rotation %g %g %g\n", asize
, natom
,
1027 targetvec
[XX
], targetvec
[YY
], targetvec
[ZZ
],
1028 aligncenter
[XX
], aligncenter
[YY
], aligncenter
[ZZ
]);
1029 /*subtract out pivot point*/
1030 for (i
= 0; i
< asize
; i
++)
1032 rvec_dec(x
[aindex
[i
]], aligncenter
);
1034 /*now determine transform and rotate*/
1036 principal_comp(asize
, aindex
, atoms
.atom
, x
, trans
, princd
);
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
[XX
] = trans
[0][2]; tmpvec
[YY
] = trans
[1][2]; tmpvec
[ZZ
] = trans
[2][2];
1041 calc_rotmatrix(tmpvec
, targetvec
, rotmatrix
);
1042 /* rotmatrix finished */
1044 for (i
= 0; i
< asize
; ++i
)
1046 mvmul(rotmatrix
, x
[aindex
[i
]], tmpvec
);
1047 copy_rvec(tmpvec
, x
[aindex
[i
]]);
1050 /*add pivot point back*/
1051 for (i
= 0; i
< asize
; i
++)
1053 rvec_inc(x
[aindex
[i
]], aligncenter
);
1065 fprintf(stderr
, "\nSelect a group that you want to translate:\n");
1066 get_index(&atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
),
1067 1, &ssize
, &sindex
, &sgrpname
);
1074 printf("Translating %d atoms (out of %d) by %g %g %g nm\n", ssize
, natom
,
1075 translation
[XX
], translation
[YY
], translation
[ZZ
]);
1078 for (i
= 0; i
< ssize
; i
++)
1080 rvec_inc(x
[sindex
[i
]], translation
);
1085 for (i
= 0; i
< natom
; i
++)
1087 rvec_inc(x
[i
], translation
);
1094 printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n", rotangles
[XX
], rotangles
[YY
], rotangles
[ZZ
]);
1095 for (i
= 0; i
< DIM
; i
++)
1097 rotangles
[i
] *= DEG2RAD
;
1099 rotate_conf(natom
, x
, v
, rotangles
[XX
], rotangles
[YY
], rotangles
[ZZ
]);
1104 /* recalc geometrical center and max and min coordinates and size */
1105 calc_geom(ssize
, sindex
, x
, gc
, min
, max
, FALSE
);
1106 rvec_sub(max
, min
, size
);
1107 if (bScale
|| bOrient
|| bRotate
)
1109 printf("new system size : %6.3f %6.3f %6.3f\n",
1110 size
[XX
], size
[YY
], size
[ZZ
]);
1114 if (bSetSize
|| bDist
|| (btype
[0][0] == 't' && bSetAng
))
1117 if (!(bSetSize
|| bDist
))
1119 for (i
= 0; i
< DIM
; i
++)
1121 newbox
[i
] = norm(box
[i
]);
1125 /* calculate new boxsize */
1126 switch (btype
[0][0])
1131 for (i
= 0; i
< DIM
; i
++)
1133 newbox
[i
] = size
[i
]+2*dist
;
1138 box
[XX
][XX
] = newbox
[XX
];
1139 box
[YY
][YY
] = newbox
[YY
];
1140 box
[ZZ
][ZZ
] = newbox
[ZZ
];
1144 matrix_convert(box
, newbox
, newang
);
1158 if (btype
[0][0] == 'c')
1160 for (i
= 0; i
< DIM
; i
++)
1165 else if (btype
[0][0] == 'd')
1171 box
[ZZ
][ZZ
] = d
*sqrt(2)/2;
1177 box
[YY
][YY
] = d
*sqrt(2)*2/3;
1179 box
[ZZ
][YY
] = d
*sqrt(2)/3;
1180 box
[ZZ
][ZZ
] = d
*sqrt(6)/3;
1186 /* calculate new coords for geometrical center */
1189 calc_box_center(ecenterDEF
, box
, center
);
1192 /* center molecule on 'center' */
1195 center_conf(natom
, x
, center
, gc
);
1201 calc_geom(ssize
, sindex
, x
, gc
, min
, max
, FALSE
);
1202 printf("new center :%7.3f%7.3f%7.3f (nm)\n", gc
[XX
], gc
[YY
], gc
[ZZ
]);
1204 if (bOrient
|| bScale
|| bDist
|| bSetSize
)
1206 printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n",
1207 norm(box
[XX
]), norm(box
[YY
]), norm(box
[ZZ
]));
1208 printf("new box angles :%7.2f%7.2f%7.2f (degrees)\n",
1209 norm2(box
[ZZ
]) == 0 ? 0 :
1210 RAD2DEG
*acos(cos_angle_no_table(box
[YY
], box
[ZZ
])),
1211 norm2(box
[ZZ
]) == 0 ? 0 :
1212 RAD2DEG
*acos(cos_angle_no_table(box
[XX
], box
[ZZ
])),
1213 norm2(box
[YY
]) == 0 ? 0 :
1214 RAD2DEG
*acos(cos_angle_no_table(box
[XX
], box
[YY
])));
1215 printf("new box volume :%7.2f (nm^3)\n", det(box
));
1218 if (check_box(epbcXYZ
, box
))
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
));
1226 if (bDist
&& btype
[0][0] == 't')
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
);
1235 else if (!opt2parg_bSet("-bt", NPA
, pa
))
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");
1242 if (bCONECT
&& (outftp
== efPDB
) && (inftp
== efTPR
))
1244 conect
= gmx_conect_generate(top
);
1253 fprintf(stderr
, "\nSelect a group for output:\n");
1254 get_index(&atoms
, opt2fn_null("-n", NFILE
, fnm
),
1255 1, &isize
, &index
, &grpname
);
1257 if (resnr_start
>= 0)
1259 renum_resnr(&atoms
, isize
, index
, resnr_start
);
1262 if (opt2parg_bSet("-label", NPA
, pa
))
1264 for (i
= 0; (i
< atoms
.nr
); i
++)
1266 atoms
.resinfo
[atoms
.atom
[i
].resind
].chainid
= label
[0];
1270 if (opt2bSet("-bf", NFILE
, fnm
) || bLegend
)
1272 gmx_fatal(FARGS
, "Sorry, cannot do bfactors with an index group.");
1275 if (outftp
== efPDB
)
1277 out
= gmx_ffopen(outfile
, "w");
1278 write_pdbfile_indexed(out
, title
, &atoms
, x
, ePBC
, box
, ' ', 1, isize
, index
, conect
, TRUE
);
1283 write_sto_conf_indexed(outfile
, title
, &atoms
, x
, bHaveV
? v
: NULL
, ePBC
, box
, isize
, index
);
1288 if (resnr_start
>= 0)
1290 renum_resnr(&atoms
, atoms
.nr
, NULL
, resnr_start
);
1293 if ((outftp
== efPDB
) || (outftp
== efPQR
))
1295 out
= gmx_ffopen(outfile
, "w");
1298 set_pdb_wide_format(TRUE
);
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");
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");
1312 else if (opt2bSet("-bf", NFILE
, fnm
))
1314 read_bfac(opt2fn("-bf", NFILE
, fnm
), &n_bfac
, &bfac
, &bfac_nr
);
1315 set_pdb_conf_bfac(atoms
.nr
, atoms
.nres
, &atoms
,
1316 n_bfac
, bfac
, bfac_nr
, peratom
);
1318 if (opt2parg_bSet("-label", NPA
, pa
))
1320 for (i
= 0; (i
< atoms
.nr
); i
++)
1322 atoms
.resinfo
[atoms
.atom
[i
].resind
].chainid
= label
[0];
1325 write_pdbfile(out
, title
, &atoms
, x
, ePBC
, box
, ' ', -1, conect
, TRUE
);
1328 pdb_legend(out
, atoms
.nr
, atoms
.nres
, &atoms
, x
);
1332 visualize_box(out
, bLegend
? atoms
.nr
+12 : atoms
.nr
,
1333 bLegend
? atoms
.nres
= 12 : atoms
.nres
, box
, visbox
);
1339 write_sto_conf(outfile
, title
, &atoms
, x
, bHaveV
? v
: NULL
, ePBC
, box
);
1342 gmx_atomprop_destroy(aps
);
1344 do_view(oenv
, outfile
, NULL
);