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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
85 real calc_mass(t_atoms *atoms, gmx_bool bGetMass, gmx_atomprop_t aps)
91 for (i = 0; (i < atoms->nr); i++)
95 gmx_atomprop_query(aps, epropMass,
96 *atoms->resinfo[atoms->atom[i].resind].name,
97 *atoms->atomname[i], &(atoms->atom[i].m));
99 tmass += atoms->atom[i].m;
105 real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,
106 rvec max, gmx_bool bDiam)
112 clear_rvec(geom_center);
125 for (j = 0; j < DIM; j++)
126 min[j] = max[j] = x[ii][j];
127 for (i = 0; i < isize; i++)
133 rvec_inc(geom_center, x[ii]);
134 for (j = 0; j < DIM; j++)
136 if (x[ii][j] < min[j])
138 if (x[ii][j] > max[j])
144 for (j = i + 1; j < isize; j++)
146 d = distance2(x[ii], x[index[j]]);
147 diam2 = max(d,diam2);
150 for (j = i + 1; j < isize; j++)
152 d = distance2(x[i], x[j]);
153 diam2 = max(d,diam2);
157 svmul(1.0 / isize, geom_center, geom_center);
163 void center_conf(int natom, rvec *x, rvec center, rvec geom_cent)
168 rvec_sub(center, geom_cent, shift);
170 printf(" shift :%7.3f%7.3f%7.3f (nm)\n", shift[XX], shift[YY],
173 for (i = 0; (i < natom); i++)
174 rvec_inc(x[i], shift);
177 void scale_conf(int natom, rvec x[], matrix box, rvec scale)
181 for (i = 0; i < natom; i++)
183 for (j = 0; j < DIM; j++)
186 for (i = 0; i < DIM; i++)
187 for (j = 0; j < DIM; j++)
188 box[i][j] *= scale[j];
191 void read_bfac(const char *fn, int *n_bfac, double **bfac_val, int **bfac_nr)
196 *n_bfac = get_lines(fn, &bfac_lines);
197 snew(*bfac_val, *n_bfac);
198 snew(*bfac_nr, *n_bfac);
199 fprintf(stderr, "Reading %d B-factors from %s\n", *n_bfac, fn);
200 for (i = 0; (i < *n_bfac); i++)
202 /*fprintf(stderr, "Line %d: %s",i,bfac_lines[i]);*/
203 sscanf(bfac_lines[i], "%d %lf", &(*bfac_nr)[i], &(*bfac_val)[i]);
204 /*fprintf(stderr," nr %d val %g\n",(*bfac_nr)[i],(*bfac_val)[i]);*/
209 void set_pdb_conf_bfac(int natoms, int nres, t_atoms *atoms, int n_bfac,
210 double *bfac, int *bfac_nr, gmx_bool peratom)
213 real bfac_min, bfac_max;
219 for (i = 0; (i < n_bfac); i++)
221 if (bfac_nr[i] - 1 >= atoms->nres)
223 /* if ((bfac_nr[i]-1<0) || (bfac_nr[i]-1>=atoms->nr))
224 gmx_fatal(FARGS,"Index of B-Factor %d is out of range: %d (%g)",
225 i+1,bfac_nr[i],bfac[i]); */
226 if (bfac[i] > bfac_max)
228 if (bfac[i] < bfac_min)
231 while ((bfac_max > 99.99) || (bfac_min < -99.99))
234 "Range of values for B-factors too large (min %g, max %g) "
235 "will scale down a factor 10\n", bfac_min, bfac_max);
236 for (i = 0; (i < n_bfac); i++)
241 while ((fabs(bfac_max) < 0.5) && (fabs(bfac_min) < 0.5))
244 "Range of values for B-factors too small (min %g, max %g) "
245 "will scale up a factor 10\n", bfac_min, bfac_max);
246 for (i = 0; (i < n_bfac); i++)
252 for (i = 0; (i < natoms); i++)
253 atoms->pdbinfo[i].bfac = 0;
257 fprintf(stderr, "Will attach %d B-factors to %d residues\n", n_bfac,
259 for (i = 0; (i < n_bfac); i++)
262 for (n = 0; (n < natoms); n++)
263 if (bfac_nr[i] == atoms->resinfo[atoms->atom[n].resind].nr)
265 atoms->pdbinfo[n].bfac = bfac[i];
270 gmx_warning("Residue nr %d not found\n", bfac_nr[i]);
276 fprintf(stderr, "Will attach %d B-factors to %d atoms\n", n_bfac,
278 for (i = 0; (i < n_bfac); i++)
280 atoms->pdbinfo[bfac_nr[i] - 1].bfac = bfac[i];
285 void pdb_legend(FILE *out, int natoms, int nres, t_atoms *atoms, rvec x[])
287 real bfac_min, bfac_max, xmin, ymin, zmin;
296 for (i = 0; (i < natoms); i++)
298 xmin = min(xmin,x[i][XX]);
299 ymin = min(ymin,x[i][YY]);
300 zmin = min(zmin,x[i][ZZ]);
301 bfac_min = min(bfac_min,atoms->pdbinfo[i].bfac);
302 bfac_max = max(bfac_max,atoms->pdbinfo[i].bfac);
304 fprintf(stderr, "B-factors range from %g to %g\n", bfac_min, bfac_max);
305 for (i = 1; (i < 12); i++)
308 "%-6s%5u %-4.4s%3.3s %c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f\n",
309 "ATOM ", natoms + 1 + i, "CA", "LEG", space, nres + 1, space,
310 (xmin + (i * 0.12)) * 10, ymin * 10, zmin * 10, 1.0, bfac_min
311 + ((i - 1.0) * (bfac_max - bfac_min) / 10));
315 void visualize_images(const char *fn, int ePBC, matrix box)
323 init_t_atoms(&atoms, nat, FALSE);
326 /* FIXME: Constness should not be cast away */
328 ala = (char *) "ALA";
329 for (i = 0; i < nat; i++)
331 atoms.atomname[i] = &c;
332 atoms.atom[i].resind = i;
333 atoms.resinfo[i].name = &ala;
334 atoms.resinfo[i].nr = i + 1;
335 atoms.resinfo[i].chainid = 'A' + i / NCUCVERT;
337 calc_triclinic_images(box, img + 1);
339 write_sto_conf(fn, "Images", &atoms, img, NULL, ePBC, box);
341 free_t_atoms(&atoms, FALSE);
345 void visualize_box(FILE *out, int a0, int r0, matrix box, rvec gridsize)
349 int nx, ny, nz, nbox, nat;
352 { 0, 1, 1, 3, 3, 2, 0, 2, 0, 4, 1, 5, 3, 7, 2, 6, 4, 5, 5, 7, 7, 6, 6,
358 nx = (int) (gridsize[XX] + 0.5);
359 ny = (int) (gridsize[YY] + 0.5);
360 nz = (int) (gridsize[ZZ] + 0.5);
364 nat = nbox * NCUCVERT;
366 calc_compact_unitcell_vertices(ecenterDEF, box, vert);
368 for (z = 0; z < nz; z++)
369 for (y = 0; y < ny; y++)
370 for (x = 0; x < nx; x++)
372 for (i = 0; i < DIM; i++)
373 shift[i] = x * box[0][i] + y * box[1][i] + z
375 for (i = 0; i < NCUCVERT; i++)
377 rvec_add(vert[i], shift, vert[j]);
382 for (i = 0; i < nat; i++)
384 fprintf(out, pdbformat, "ATOM", a0 + i, "C", "BOX", 'K' + i
385 / NCUCVERT, r0 + i, 10 * vert[i][XX], 10 * vert[i][YY], 10
390 edge = compact_unitcell_edges();
391 for (j = 0; j < nbox; j++)
392 for (i = 0; i < NCUCEDGE; i++)
393 fprintf(out, "CONECT%5d%5d\n", a0 + j * NCUCVERT + edge[2 * i],
394 a0 + j * NCUCVERT + edge[2 * i + 1]);
401 for (z = 0; z <= 1; z++)
402 for (y = 0; y <= 1; y++)
403 for (x = 0; x <= 1; x++)
405 fprintf(out, pdbformat, "ATOM", a0 + i, "C", "BOX", 'K' + i
406 / 8, r0 + i, x * 10 * box[XX][XX],
407 y * 10 * box[YY][YY], z * 10 * box[ZZ][ZZ]);
411 for (i = 0; i < 24; i += 2)
412 fprintf(out, "CONECT%5d%5d\n", a0 + rectedge[i], a0 + rectedge[i
417 void calc_rotmatrix(rvec principal_axis, rvec targetvec, matrix rotmatrix)
420 real ux,uy,uz,costheta,sintheta;
422 costheta = cos_angle(principal_axis,targetvec);
423 sintheta=sqrt(1.0-costheta*costheta); /* sign is always positive since 0<theta<pi */
425 /* Determine rotation from cross product with target vector */
426 cprod(principal_axis,targetvec,rotvec);
427 unitv(rotvec,rotvec);
428 printf("Aligning %g %g %g to %g %g %g : xprod %g %g %g\n",
429 principal_axis[XX],principal_axis[YY],principal_axis[ZZ],targetvec[XX],targetvec[YY],targetvec[ZZ],
430 rotvec[XX],rotvec[YY],rotvec[ZZ]);
435 rotmatrix[0][0]=ux*ux + (1.0-ux*ux)*costheta;
436 rotmatrix[0][1]=ux*uy*(1-costheta)-uz*sintheta;
437 rotmatrix[0][2]=ux*uz*(1-costheta)+uy*sintheta;
438 rotmatrix[1][0]=ux*uy*(1-costheta)+uz*sintheta;
439 rotmatrix[1][1]=uy*uy + (1.0-uy*uy)*costheta;
440 rotmatrix[1][2]=uy*uz*(1-costheta)-ux*sintheta;
441 rotmatrix[2][0]=ux*uz*(1-costheta)-uy*sintheta;
442 rotmatrix[2][1]=uy*uz*(1-costheta)+ux*sintheta;
443 rotmatrix[2][2]=uz*uz + (1.0-uz*uz)*costheta;
445 printf("Rotation matrix: \n%g %g %g\n%g %g %g\n%g %g %g\n",
446 rotmatrix[0][0],rotmatrix[0][1],rotmatrix[0][2],
447 rotmatrix[1][0],rotmatrix[1][1],rotmatrix[1][2],
448 rotmatrix[2][0],rotmatrix[2][1],rotmatrix[2][2]);
451 static void renum_resnr(t_atoms *atoms,int isize,const int *index,
454 int i,resind_prev,resind;
457 for(i=0; i<isize; i++)
459 resind = atoms->atom[index == NULL ? i : index[i]].resind;
460 if (resind != resind_prev)
462 atoms->resinfo[resind].nr = resnr_start;
465 resind_prev = resind;
469 int gmx_editconf(int argc, char *argv[])
474 "[TT]editconf[tt] converts generic structure format to [TT].gro[tt], [TT].g96[tt]",
477 "The box can be modified with options [TT]-box[tt], [TT]-d[tt] and",
478 "[TT]-angles[tt]. Both [TT]-box[tt] and [TT]-d[tt]",
479 "will center the system in the box, unless [TT]-noc[tt] is used.",
481 "Option [TT]-bt[tt] determines the box type: [TT]triclinic[tt] is a",
482 "triclinic box, [TT]cubic[tt] is a rectangular box with all sides equal",
483 "[TT]dodecahedron[tt] represents a rhombic dodecahedron and",
484 "[TT]octahedron[tt] is a truncated octahedron.",
485 "The last two are special cases of a triclinic box.",
486 "The length of the three box vectors of the truncated octahedron is the",
487 "shortest distance between two opposite hexagons.",
488 "Relative to a cubic box with some periodic image distance, the volume of a ",
489 "dodecahedron with this same periodic distance is 0.71 times that of the cube, ",
490 "and that of a truncated octahedron is 0.77 times.",
492 "Option [TT]-box[tt] requires only",
493 "one value for a cubic, rhombic dodecahedral, or truncated octahedral box.",
495 "With [TT]-d[tt] and a [TT]triclinic[tt] box the size of the system in the [IT]x[it]-, [IT]y[it]-,",
496 "and [IT]z[it]-directions is used. With [TT]-d[tt] and [TT]cubic[tt],",
497 "[TT]dodecahedron[tt] or [TT]octahedron[tt] boxes, the dimensions are set",
498 "to the diameter of the system (largest distance between atoms) plus twice",
499 "the specified distance.",
501 "Option [TT]-angles[tt] is only meaningful with option [TT]-box[tt] and",
502 "a triclinic box and cannot be used with option [TT]-d[tt].",
504 "When [TT]-n[tt] or [TT]-ndef[tt] is set, a group",
505 "can be selected for calculating the size and the geometric center,",
506 "otherwise the whole system is used.",
508 "[TT]-rotate[tt] rotates the coordinates and velocities.",
510 "[TT]-princ[tt] aligns the principal axes of the system along the",
511 "coordinate axes, with the longest axis aligned with the [IT]x[it]-axis. ",
512 "This may allow you to decrease the box volume,",
513 "but beware that molecules can rotate significantly in a nanosecond.",
515 "Scaling is applied before any of the other operations are",
516 "performed. Boxes and coordinates can be scaled to give a certain density (option",
517 "[TT]-density[tt]). Note that this may be inaccurate in case a [TT].gro[tt]",
518 "file is given as input. A special feature of the scaling option is that when the",
519 "factor -1 is given in one dimension, one obtains a mirror image,",
520 "mirrored in one of the planes. When one uses -1 in three dimensions, ",
521 "a point-mirror image is obtained.[PAR]",
522 "Groups are selected after all operations have been applied.[PAR]",
523 "Periodicity can be removed in a crude manner.",
524 "It is important that the box vectors at the bottom of your input file",
525 "are correct when the periodicity is to be removed.",
527 "When writing [TT].pdb[tt] files, B-factors can be",
528 "added with the [TT]-bf[tt] option. B-factors are read",
529 "from a file with with following format: first line states number of",
530 "entries in the file, next lines state an index",
531 "followed by a B-factor. The B-factors will be attached per residue",
532 "unless an index is larger than the number of residues or unless the",
533 "[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
534 "be added instead of B-factors. [TT]-legend[tt] will produce",
535 "a row of CA atoms with B-factors ranging from the minimum to the",
536 "maximum value found, effectively making a legend for viewing.",
538 "With the option [TT]-mead[tt] a special [TT].pdb[tt] ([TT].pqr[tt])",
539 "file for the MEAD electrostatics",
540 "program (Poisson-Boltzmann solver) can be made. A further prerequisite",
541 "is that the input file is a run input file.",
542 "The B-factor field is then filled with the Van der Waals radius",
543 "of the atoms while the occupancy field will hold the charge.",
545 "The option [TT]-grasp[tt] is similar, but it puts the charges in the B-factor",
546 "and the radius in the occupancy.",
548 "Option [TT]-align[tt] allows alignment",
549 "of the principal axis of a specified group against the given vector, ",
550 "with an optional center of rotation specified by [TT]-aligncenter[tt].",
552 "Finally, with option [TT]-label[tt], [TT]editconf[tt] can add a chain identifier",
553 "to a [TT].pdb[tt] file, which can be useful for analysis with e.g. Rasmol.",
555 "To convert a truncated octrahedron file produced by a package which uses",
556 "a cubic box with the corners cut off (such as GROMOS), use:[BR]",
557 "[TT]editconf -f in -rotate 0 45 35.264 -bt o -box veclen -o out[tt][BR]",
558 "where [TT]veclen[tt] is the size of the cubic box times [SQRT]3[sqrt]/2." };
561 "For complex molecules, the periodicity removal routine may break down, "
562 "in that case you can use [TT]trjconv[tt]." };
563 static real dist = 0.0, rbox = 0.0, to_diam = 0.0;
564 static gmx_bool bNDEF = FALSE, bRMPBC = FALSE, bCenter = FALSE, bReadVDW =
565 FALSE, bCONECT = FALSE;
566 static gmx_bool peratom = FALSE, bLegend = FALSE, bOrient = FALSE, bMead =
567 FALSE, bGrasp = FALSE, bSig56 = FALSE;
569 { 1, 1, 1 }, newbox =
570 { 0, 0, 0 }, newang =
572 static real rho = 1000.0, rvdw = 0.12;
574 { 0, 0, 0 }, translation =
575 { 0, 0, 0 }, rotangles =
576 { 0, 0, 0 }, aligncenter =
577 { 0, 0, 0 }, targetvec =
579 static const char *btype[] =
580 { NULL, "triclinic", "cubic", "dodecahedron", "octahedron", NULL },
584 static int resnr_start = -1;
588 { "-ndef", FALSE, etBOOL,
589 { &bNDEF }, "Choose output from default index groups" },
590 { "-visbox", FALSE, etRVEC,
592 "HIDDENVisualize a grid of boxes, -1 visualizes the 14 box images" },
593 { "-bt", FALSE, etENUM,
594 { btype }, "Box type for [TT]-box[tt] and [TT]-d[tt]" },
595 { "-box", FALSE, etRVEC,
596 { newbox }, "Box vector lengths (a,b,c)" },
597 { "-angles", FALSE, etRVEC,
598 { newang }, "Angles between the box vectors (bc,ac,ab)" },
599 { "-d", FALSE, etREAL,
600 { &dist }, "Distance between the solute and the box" },
601 { "-c", FALSE, etBOOL,
603 "Center molecule in box (implied by [TT]-box[tt] and [TT]-d[tt])" },
604 { "-center", FALSE, etRVEC,
605 { center }, "Coordinates of geometrical center" },
606 { "-aligncenter", FALSE, etRVEC,
607 { aligncenter }, "Center of rotation for alignment" },
608 { "-align", FALSE, etRVEC,
610 "Align to target vector" },
611 { "-translate", FALSE, etRVEC,
612 { translation }, "Translation" },
613 { "-rotate", FALSE, etRVEC,
615 "Rotation around the X, Y and Z axes in degrees" },
616 { "-princ", FALSE, etBOOL,
618 "Orient molecule(s) along their principal axes" },
619 { "-scale", FALSE, etRVEC,
620 { scale }, "Scaling factor" },
621 { "-density", FALSE, etREAL,
623 "Density (g/L) of the output box achieved by scaling" },
624 { "-pbc", FALSE, etBOOL,
626 "Remove the periodicity (make molecule whole again)" },
627 { "-resnr", FALSE, etINT,
629 " Renumber residues starting from resnr" },
630 { "-grasp", FALSE, etBOOL,
632 "Store the charge of the atom in the B-factor field and the radius of the atom in the occupancy field" },
634 "-rvdw", FALSE, etREAL,
636 "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" },
637 { "-sig56", FALSE, etBOOL,
639 "Use rmin/2 (minimum in the Van der Waals potential) rather than [GRK]sigma[grk]/2 " },
641 "-vdwread", FALSE, etBOOL,
643 "Read the Van der Waals radii from the file [TT]vdwradii.dat[tt] rather than computing the radii based on the force field" },
644 { "-atom", FALSE, etBOOL,
645 { &peratom }, "Force B-factor attachment per atom" },
646 { "-legend", FALSE, etBOOL,
647 { &bLegend }, "Make B-factor legend" },
648 { "-label", FALSE, etSTR,
649 { &label }, "Add chain label for all residues" },
651 "-conect", FALSE, etBOOL,
653 "Add CONECT records to a [TT].pdb[tt] file when written. Can only be done when a topology is present" } };
654 #define NPA asize(pa)
657 const char *infile, *outfile;
659 int outftp, inftp, natom, i, j, n_bfac, itype, ntype;
660 double *bfac = NULL, c6, c12;
662 t_topology *top = NULL;
664 char *grpname, *sgrpname, *agrpname;
665 int isize, ssize, tsize, asize;
666 atom_id *index, *sindex, *tindex, *aindex;
667 rvec *x, *v, gc, min, max, size;
669 matrix box,rotmatrix,trans;
671 gmx_bool bIndex, bSetSize, bSetAng, bCubic, bDist, bSetCenter, bAlign;
672 gmx_bool bHaveV, bScale, bRho, bTranslate, bRotate, bCalcGeom, bCalcDiam;
673 real xs, ys, zs, xcent, ycent, zcent, diam = 0, mass = 0, d, vdw;
679 { efSTX, "-f", NULL, ffREAD },
680 { efNDX, "-n", NULL, ffOPTRD },
681 { efSTO, NULL, NULL, ffOPTWR },
682 { efPQR, "-mead", "mead", ffOPTWR },
683 { efDAT, "-bf", "bfact", ffOPTRD } };
684 #define NFILE asize(fnm)
686 CopyRight(stderr, argv[0]);
687 parse_common_args(&argc, argv, PCA_CAN_VIEW, NFILE, fnm, NPA, pa,
688 asize(desc), desc, asize(bugs), bugs, &oenv);
690 bIndex = opt2bSet("-n", NFILE, fnm) || bNDEF;
691 bMead = opt2bSet("-mead", NFILE, fnm);
692 bSetSize = opt2parg_bSet("-box", NPA, pa);
693 bSetAng = opt2parg_bSet("-angles", NPA, pa);
694 bSetCenter = opt2parg_bSet("-center", NPA, pa);
695 bDist = opt2parg_bSet("-d", NPA, pa);
696 bAlign = opt2parg_bSet("-align", NPA, pa);
697 /* Only automatically turn on centering without -noc */
698 if ((bDist || bSetSize || bSetCenter) && !opt2parg_bSet("-c", NPA, pa))
702 bScale = opt2parg_bSet("-scale", NPA, pa);
703 bRho = opt2parg_bSet("-density", NPA, pa);
704 bTranslate = opt2parg_bSet("-translate", NPA, pa);
705 bRotate = opt2parg_bSet("-rotate", NPA, pa);
707 fprintf(stderr, "WARNING: setting -density overrides -scale\n");
708 bScale = bScale || bRho;
709 bCalcGeom = bCenter || bRotate || bOrient || bScale;
710 bCalcDiam = btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o';
712 infile = ftp2fn(efSTX, NFILE, fnm);
714 outfile = ftp2fn(efPQR, NFILE, fnm);
716 outfile = ftp2fn(efSTO, NFILE, fnm);
717 outftp = fn2ftp(outfile);
718 inftp = fn2ftp(infile);
720 aps = gmx_atomprop_init();
724 printf("Incompatible options -mead and -grasp. Turning off -grasp\n");
727 if (bGrasp && (outftp != efPDB))
728 gmx_fatal(FARGS, "Output file should be a .pdb file"
729 " when using the -grasp option\n");
730 if ((bMead || bGrasp) && !((fn2ftp(infile) == efTPR) ||
731 (fn2ftp(infile) == efTPA) ||
732 (fn2ftp(infile) == efTPB)))
733 gmx_fatal(FARGS,"Input file should be a .tp[abr] file"
734 " when using the -mead option\n");
736 get_stx_coordnum(infile,&natom);
737 init_t_atoms(&atoms,natom,TRUE);
740 read_stx_conf(infile,title,&atoms,x,v,&ePBC,box);
741 if (fn2ftp(infile) == efPDB)
743 get_pdb_atomnumber(&atoms,aps);
745 printf("Read %d atoms\n",atoms.nr);
747 /* Get the element numbers if available in a pdb file */
748 if (fn2ftp(infile) == efPDB)
749 get_pdb_atomnumber(&atoms,aps);
751 if (ePBC != epbcNONE)
754 printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
755 vol,100*((int)(vol*4.5)));
758 if (bMead || bGrasp || bCONECT)
759 top = read_top(infile,NULL);
763 if (atoms.nr != top->atoms.nr)
764 gmx_fatal(FARGS,"Atom numbers don't match (%d vs. %d)",atoms.nr,top->atoms.nr);
765 snew(atoms.pdbinfo,top->atoms.nr);
766 ntype = top->idef.atnr;
767 for(i=0; (i<atoms.nr); i++) {
768 /* Determine the Van der Waals radius from the force field */
770 if (!gmx_atomprop_query(aps,epropVDW,
771 *top->atoms.resinfo[top->atoms.atom[i].resind].name,
772 *top->atoms.atomname[i],&vdw))
776 itype = top->atoms.atom[i].type;
777 c12 = top->idef.iparams[itype*ntype+itype].lj.c12;
778 c6 = top->idef.iparams[itype*ntype+itype].lj.c6;
779 if ((c6 != 0) && (c12 != 0)) {
785 vdw = 0.5*pow(sig6,1.0/6.0);
790 /* Factor of 10 for nm -> Angstroms */
794 atoms.pdbinfo[i].occup = top->atoms.atom[i].q;
795 atoms.pdbinfo[i].bfac = vdw;
798 atoms.pdbinfo[i].occup = vdw;
799 atoms.pdbinfo[i].bfac = top->atoms.atom[i].q;
804 for (i=0; (i<natom) && !bHaveV; i++)
805 for (j=0; (j<DIM) && !bHaveV; j++)
806 bHaveV=bHaveV || (v[i][j]!=0);
807 printf("%selocities found\n",bHaveV?"V":"No v");
811 gmx_fatal(FARGS,"Sorry, can not visualize box with index groups");
813 gmx_fatal(FARGS,"Sorry, can only visualize box with a pdb file");
814 } else if (visbox[0] == -1)
815 visualize_images("images.pdb",ePBC,box);
819 rm_gropbc(&atoms,x,box);
823 fprintf(stderr,"\nSelect a group for determining the system size:\n");
824 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
825 1,&ssize,&sindex,&sgrpname);
830 diam=calc_geom(ssize,sindex,x,gc,min,max,bCalcDiam);
831 rvec_sub(max, min, size);
832 printf(" system size :%7.3f%7.3f%7.3f (nm)\n",
833 size[XX], size[YY], size[ZZ]);
835 printf(" diameter :%7.3f (nm)\n",diam);
836 printf(" center :%7.3f%7.3f%7.3f (nm)\n", gc[XX], gc[YY], gc[ZZ]);
837 printf(" box vectors :%7.3f%7.3f%7.3f (nm)\n",
838 norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
839 printf(" box angles :%7.2f%7.2f%7.2f (degrees)\n",
840 norm2(box[ZZ])==0 ? 0 :
841 RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
842 norm2(box[ZZ])==0 ? 0 :
843 RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
844 norm2(box[YY])==0 ? 0 :
845 RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
846 printf(" box volume :%7.2f (nm^3)\n",det(box));
849 if (bRho || bOrient || bAlign)
850 mass = calc_mass(&atoms,!fn2bTPX(infile),aps);
856 /* Get a group for principal component analysis */
857 fprintf(stderr,"\nSelect group for the determining the orientation\n");
858 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpnames);
860 /* Orient the principal axes along the coordinate axes */
861 orient_princ(&atoms,isize,index,natom,x,bHaveV ? v : NULL, NULL);
867 /* scale coordinates and box */
869 /* Compute scaling constant */
873 dens = (mass*AMU)/(vol*NANO*NANO*NANO);
874 fprintf(stderr,"Volume of input %g (nm^3)\n",vol);
875 fprintf(stderr,"Mass of input %g (a.m.u.)\n",mass);
876 fprintf(stderr,"Density of input %g (g/l)\n",dens);
877 if (vol==0 || mass==0)
878 gmx_fatal(FARGS,"Cannot scale density with "
879 "zero mass (%g) or volume (%g)\n",mass,vol);
881 scale[XX] = scale[YY] = scale[ZZ] = pow(dens/rho,1.0/3.0);
882 fprintf(stderr,"Scaling all box vectors by %g\n",scale[XX]);
884 scale_conf(atoms.nr,x,box,scale);
889 fprintf(stderr,"\nSelect a group that you want to align:\n");
890 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
891 1,&asize,&aindex,&agrpname);
895 for (i=0;i<asize;i++)
898 printf("Aligning %d atoms (out of %d) to %g %g %g, center of rotation %g %g %g\n",asize,natom,
899 targetvec[XX],targetvec[YY],targetvec[ZZ],
900 aligncenter[XX],aligncenter[YY],aligncenter[ZZ]);
901 /*subtract out pivot point*/
902 for(i=0; i<asize; i++)
903 rvec_dec(x[aindex[i]],aligncenter);
904 /*now determine transform and rotate*/
906 principal_comp(asize,aindex,atoms.atom,x, trans,princd);
908 unitv(targetvec,targetvec);
909 printf("Using %g %g %g as principal axis\n", trans[0][2],trans[1][2],trans[2][2]);
910 tmpvec[XX]=trans[0][2]; tmpvec[YY]=trans[1][2]; tmpvec[ZZ]=trans[2][2];
911 calc_rotmatrix(tmpvec, targetvec, rotmatrix);
912 /* rotmatrix finished */
914 for (i=0;i<asize;++i)
916 mvmul(rotmatrix,x[aindex[i]],tmpvec);
917 copy_rvec(tmpvec,x[aindex[i]]);
920 /*add pivot point back*/
921 for(i=0; i<asize; i++)
922 rvec_inc(x[aindex[i]],aligncenter);
929 fprintf(stderr,"\nSelect a group that you want to translate:\n");
930 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
931 1,&ssize,&sindex,&sgrpname);
936 printf("Translating %d atoms (out of %d) by %g %g %g nm\n",ssize,natom,
937 translation[XX],translation[YY],translation[ZZ]);
939 for(i=0; i<ssize; i++)
940 rvec_inc(x[sindex[i]],translation);
943 for(i=0; i<natom; i++)
944 rvec_inc(x[i],translation);
949 printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n",rotangles[XX],rotangles[YY],rotangles[ZZ]);
951 rotangles[i] *= DEG2RAD;
952 rotate_conf(natom,x,v,rotangles[XX],rotangles[YY],rotangles[ZZ]);
956 /* recalc geometrical center and max and min coordinates and size */
957 calc_geom(ssize,sindex,x,gc,min,max,FALSE);
958 rvec_sub(max, min, size);
959 if (bScale || bOrient || bRotate)
960 printf("new system size : %6.3f %6.3f %6.3f\n",
961 size[XX],size[YY],size[ZZ]);
964 if (bSetSize || bDist || (btype[0][0]=='t' && bSetAng)) {
966 if (!(bSetSize || bDist))
967 for (i=0; i<DIM; i++)
968 newbox[i] = norm(box[i]);
970 /* calculate new boxsize */
975 newbox[i] = size[i]+2*dist;
977 box[XX][XX] = newbox[XX];
978 box[YY][YY] = newbox[YY];
979 box[ZZ][ZZ] = newbox[ZZ];
981 matrix_convert(box,newbox,newang);
991 if (btype[0][0] == 'c')
994 else if (btype[0][0] == 'd') {
999 box[ZZ][ZZ] = d*sqrt(2)/2;
1003 box[YY][YY] = d*sqrt(2)*2/3;
1005 box[ZZ][YY] = d*sqrt(2)/3;
1006 box[ZZ][ZZ] = d*sqrt(6)/3;
1012 /* calculate new coords for geometrical center */
1014 calc_box_center(ecenterDEF,box,center);
1016 /* center molecule on 'center' */
1018 center_conf(natom,x,center,gc);
1022 calc_geom(ssize,sindex,x, gc, min, max, FALSE);
1023 printf("new center :%7.3f%7.3f%7.3f (nm)\n",gc[XX],gc[YY],gc[ZZ]);
1025 if (bOrient || bScale || bDist || bSetSize) {
1026 printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n",
1027 norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
1028 printf("new box angles :%7.2f%7.2f%7.2f (degrees)\n",
1029 norm2(box[ZZ])==0 ? 0 :
1030 RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
1031 norm2(box[ZZ])==0 ? 0 :
1032 RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
1033 norm2(box[YY])==0 ? 0 :
1034 RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
1035 printf("new box volume :%7.2f (nm^3)\n",det(box));
1038 if (check_box(epbcXYZ,box))
1039 printf("\nWARNING: %s\n"
1040 "See the GROMACS manual for a description of the requirements that\n"
1041 "must be satisfied by descriptions of simulation cells.\n",
1042 check_box(epbcXYZ,box));
1044 if (bDist && btype[0][0]=='t')
1048 printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
1049 "distance from the solute to a box surface along the corresponding normal\n"
1050 "vector might be somewhat smaller than your specified value %f.\n"
1051 "You can check the actual value with g_mindist -pi\n",dist);
1053 else if (!opt2parg_bSet("-bt", NPA, pa))
1055 printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
1056 "If the molecule rotates the actual distance will be smaller. You might want\n"
1057 "to use a cubic box instead, or why not try a dodecahedron today?\n");
1060 if (bCONECT && (outftp == efPDB) && (inftp == efTPR))
1061 conect = gmx_conect_generate(top);
1066 fprintf(stderr,"\nSelect a group for output:\n");
1067 get_index(&atoms,opt2fn_null("-n",NFILE,fnm),
1068 1,&isize,&index,&grpname);
1070 if (resnr_start >= 0)
1072 renum_resnr(&atoms,isize,index,resnr_start);
1075 if (opt2parg_bSet("-label",NPA,pa)) {
1076 for(i=0; (i<atoms.nr); i++)
1077 atoms.resinfo[atoms.atom[i].resind].chainid=label[0];
1080 if (opt2bSet("-bf",NFILE,fnm) || bLegend)
1082 gmx_fatal(FARGS,"Sorry, cannot do bfactors with an index group.");
1085 if (outftp == efPDB)
1087 out=ffopen(outfile,"w");
1088 write_pdbfile_indexed(out,title,&atoms,x,ePBC,box,' ',1,isize,index,conect,TRUE);
1093 write_sto_conf_indexed(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box,isize,index);
1098 if (resnr_start >= 0)
1100 renum_resnr(&atoms,atoms.nr,NULL,resnr_start);
1103 if ((outftp == efPDB) || (outftp == efPQR)) {
1104 out=ffopen(outfile,"w");
1106 set_pdb_wide_format(TRUE);
1107 fprintf(out,"REMARK "
1108 "The B-factors in this file hold atomic radii\n");
1109 fprintf(out,"REMARK "
1110 "The occupancy in this file hold atomic charges\n");
1113 fprintf(out,"GRASP PDB FILE\nFORMAT NUMBER=1\n");
1114 fprintf(out,"REMARK "
1115 "The B-factors in this file hold atomic charges\n");
1116 fprintf(out,"REMARK "
1117 "The occupancy in this file hold atomic radii\n");
1119 else if (opt2bSet("-bf",NFILE,fnm)) {
1120 read_bfac(opt2fn("-bf",NFILE,fnm),&n_bfac,&bfac,&bfac_nr);
1121 set_pdb_conf_bfac(atoms.nr,atoms.nres,&atoms,
1122 n_bfac,bfac,bfac_nr,peratom);
1124 if (opt2parg_bSet("-label",NPA,pa)) {
1125 for(i=0; (i<atoms.nr); i++)
1126 atoms.resinfo[atoms.atom[i].resind].chainid=label[0];
1128 write_pdbfile(out,title,&atoms,x,ePBC,box,' ',-1,conect,TRUE);
1130 pdb_legend(out,atoms.nr,atoms.nres,&atoms,x);
1132 visualize_box(out,bLegend ? atoms.nr+12 : atoms.nr,
1133 bLegend? atoms.nres=12 : atoms.nres,box,visbox);
1137 write_sto_conf(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box);
1139 gmx_atomprop_destroy(aps);
1141 do_view(oenv,outfile,NULL);