1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
82 static const char *pdbtp
[epdbNR
] =
83 { "ATOM ", "HETATM" };
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 int gmx_editconf(int argc
, char *argv
[])
456 "editconf converts generic structure format to [TT].gro[tt], [TT].g96[tt]",
459 "The box can be modified with options [TT]-box[tt], [TT]-d[tt] and",
460 "[TT]-angles[tt]. Both [TT]-box[tt] and [TT]-d[tt]",
461 "will center the system in the box, unless [TT]-noc[tt] is used.",
463 "Option [TT]-bt[tt] determines the box type: [TT]triclinic[tt] is a",
464 "triclinic box, [TT]cubic[tt] is a rectangular box with all sides equal",
465 "[TT]dodecahedron[tt] represents a rhombic dodecahedron and",
466 "[TT]octahedron[tt] is a truncated octahedron.",
467 "The last two are special cases of a triclinic box.",
468 "The length of the three box vectors of the truncated octahedron is the",
469 "shortest distance between two opposite hexagons.",
470 "The volume of a dodecahedron is 0.71 and that of a truncated octahedron",
471 "is 0.77 of that of a cubic box with the same periodic image distance.",
473 "Option [TT]-box[tt] requires only",
474 "one value for a cubic box, dodecahedron and a truncated octahedron.",
476 "With [TT]-d[tt] and a [TT]triclinic[tt] box the size of the system in the x, y",
477 "and z directions is used. With [TT]-d[tt] and [TT]cubic[tt],",
478 "[TT]dodecahedron[tt] or [TT]octahedron[tt] boxes, the dimensions are set",
479 "to the diameter of the system (largest distance between atoms) plus twice",
480 "the specified distance.",
482 "Option [TT]-angles[tt] is only meaningful with option [TT]-box[tt] and",
483 "a triclinic box and can not be used with option [TT]-d[tt].",
485 "When [TT]-n[tt] or [TT]-ndef[tt] is set, a group",
486 "can be selected for calculating the size and the geometric center,",
487 "otherwise the whole system is used.",
489 "[TT]-rotate[tt] rotates the coordinates and velocities.",
491 "[TT]-princ[tt] aligns the principal axes of the system along the",
492 "coordinate axes, this may allow you to decrease the box volume,",
493 "but beware that molecules can rotate significantly in a nanosecond.",
495 "Scaling is applied before any of the other operations are",
496 "performed. Boxes and coordinates can be scaled to give a certain density (option",
497 "[TT]-density[tt]). Note that this may be inaccurate in case a gro",
498 "file is given as input. A special feature of the scaling option, when the",
499 "factor -1 is given in one dimension, one obtains a mirror image,",
500 "mirrored in one of the plains, when one uses -1 in three dimensions",
501 "a point-mirror image is obtained.[PAR]",
502 "Groups are selected after all operations have been applied.[PAR]",
503 "Periodicity can be removed in a crude manner.",
504 "It is important that the box sizes at the bottom of your input file",
505 "are correct when the periodicity is to be removed.",
507 "When writing [TT].pdb[tt] files, B-factors can be",
508 "added with the [TT]-bf[tt] option. B-factors are read",
509 "from a file with with following format: first line states number of",
510 "entries in the file, next lines state an index",
511 "followed by a B-factor. The B-factors will be attached per residue",
512 "unless an index is larger than the number of residues or unless the",
513 "[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
514 "be added instead of B-factors. [TT]-legend[tt] will produce",
515 "a row of CA atoms with B-factors ranging from the minimum to the",
516 "maximum value found, effectively making a legend for viewing.",
518 "With the option -mead a special pdb (pqr) file for the MEAD electrostatics",
519 "program (Poisson-Boltzmann solver) can be made. A further prerequisite",
520 "is that the input file is a run input file.",
521 "The B-factor field is then filled with the Van der Waals radius",
522 "of the atoms while the occupancy field will hold the charge.",
524 "The option -grasp is similar, but it puts the charges in the B-factor",
525 "and the radius in the occupancy.",
527 "Option [TT]-align[tt] allows alignment",
528 "of the principal axis of a specified group against the given vector, ",
529 "with an optional center of rotation specified by [TT]-aligncenter[tt].",
531 "Finally with option [TT]-label[tt] editconf can add a chain identifier",
532 "to a pdb file, which can be useful for analysis with e.g. rasmol.",
534 "To convert a truncated octrahedron file produced by a package which uses",
535 "a cubic box with the corners cut off (such as Gromos) use:[BR]",
536 "[TT]editconf -f <in> -rotate 0 45 35.264 -bt o -box <veclen> -o <out>[tt][BR]",
537 "where [TT]veclen[tt] is the size of the cubic box times sqrt(3)/2." };
540 "For complex molecules, the periodicity removal routine may break down, ",
541 "in that case you can use trjconv." };
542 static real dist
= 0.0, rbox
= 0.0, to_diam
= 0.0;
543 static gmx_bool bNDEF
= FALSE
, bRMPBC
= FALSE
, bCenter
= FALSE
, bReadVDW
=
544 FALSE
, bCONECT
= FALSE
;
545 static gmx_bool peratom
= FALSE
, bLegend
= FALSE
, bOrient
= FALSE
, bMead
=
546 FALSE
, bGrasp
= FALSE
, bSig56
= FALSE
;
548 { 1, 1, 1 }, newbox
=
549 { 0, 0, 0 }, newang
=
551 static real rho
= 1000.0, rvdw
= 0.12;
553 { 0, 0, 0 }, translation
=
554 { 0, 0, 0 }, rotangles
=
555 { 0, 0, 0 }, aligncenter
=
556 { 0, 0, 0 }, targetvec
=
558 static const char *btype
[] =
559 { NULL
, "triclinic", "cubic", "dodecahedron", "octahedron", NULL
},
566 { "-ndef", FALSE
, etBOOL
,
567 { &bNDEF
}, "Choose output from default index groups" },
568 { "-visbox", FALSE
, etRVEC
,
570 "HIDDENVisualize a grid of boxes, -1 visualizes the 14 box images" },
571 { "-bt", FALSE
, etENUM
,
572 { btype
}, "Box type for -box and -d" },
573 { "-box", FALSE
, etRVEC
,
574 { newbox
}, "Box vector lengths (a,b,c)" },
575 { "-angles", FALSE
, etRVEC
,
576 { newang
}, "Angles between the box vectors (bc,ac,ab)" },
577 { "-d", FALSE
, etREAL
,
578 { &dist
}, "Distance between the solute and the box" },
579 { "-c", FALSE
, etBOOL
,
581 "Center molecule in box (implied by -box and -d)" },
582 { "-center", FALSE
, etRVEC
,
583 { center
}, "Coordinates of geometrical center" },
584 { "-aligncenter", FALSE
, etRVEC
,
585 { aligncenter
}, "Center of rotation for alignment" },
586 { "-align", FALSE
, etRVEC
,
588 "Align to target vector" },
589 { "-translate", FALSE
, etRVEC
,
590 { translation
}, "Translation" },
591 { "-rotate", FALSE
, etRVEC
,
593 "Rotation around the X, Y and Z axes in degrees" },
594 { "-princ", FALSE
, etBOOL
,
596 "Orient molecule(s) along their principal axes" },
597 { "-scale", FALSE
, etRVEC
,
598 { scale
}, "Scaling factor" },
599 { "-density", FALSE
, etREAL
,
601 "Density (g/l) of the output box achieved by scaling" },
602 { "-pbc", FALSE
, etBOOL
,
604 "Remove the periodicity (make molecule whole again)" },
605 { "-grasp", FALSE
, etBOOL
,
607 "Store the charge of the atom in the B-factor field and the radius of the atom in the occupancy field" },
609 "-rvdw", FALSE
, etREAL
,
611 "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" },
612 { "-sig56", FALSE
, etREAL
,
614 "Use rmin/2 (minimum in the Van der Waals potential) rather than sigma/2 " },
616 "-vdwread", FALSE
, etBOOL
,
618 "Read the Van der Waals radii from the file vdwradii.dat rather than computing the radii based on the force field" },
619 { "-atom", FALSE
, etBOOL
,
620 { &peratom
}, "Force B-factor attachment per atom" },
621 { "-legend", FALSE
, etBOOL
,
622 { &bLegend
}, "Make B-factor legend" },
623 { "-label", FALSE
, etSTR
,
624 { &label
}, "Add chain label for all residues" },
626 "-conect", FALSE
, etBOOL
,
628 "Add CONECT records to a pdb file when written. Can only be done when a topology is present" } };
629 #define NPA asize(pa)
632 const char *infile
, *outfile
;
634 int outftp
, inftp
, natom
, i
, j
, n_bfac
, itype
, ntype
;
635 double *bfac
= NULL
, c6
, c12
;
637 t_topology
*top
= NULL
;
639 char *grpname
, *sgrpname
, *agrpname
;
640 int isize
, ssize
, tsize
, asize
;
641 atom_id
*index
, *sindex
, *tindex
, *aindex
;
642 rvec
*x
, *v
, gc
, min
, max
, size
;
644 matrix box
,rotmatrix
,trans
;
646 gmx_bool bIndex
, bSetSize
, bSetAng
, bCubic
, bDist
, bSetCenter
, bAlign
;
647 gmx_bool bHaveV
, bScale
, bRho
, bTranslate
, bRotate
, bCalcGeom
, bCalcDiam
;
648 real xs
, ys
, zs
, xcent
, ycent
, zcent
, diam
= 0, mass
= 0, d
, vdw
;
654 { efSTX
, "-f", NULL
, ffREAD
},
655 { efNDX
, "-n", NULL
, ffOPTRD
},
656 { efSTO
, NULL
, NULL
, ffOPTWR
},
657 { efPQR
, "-mead", "mead", ffOPTWR
},
658 { efDAT
, "-bf", "bfact", ffOPTRD
} };
659 #define NFILE asize(fnm)
661 CopyRight(stderr
, argv
[0]);
662 parse_common_args(&argc
, argv
, PCA_CAN_VIEW
, NFILE
, fnm
, NPA
, pa
,
663 asize(desc
), desc
, asize(bugs
), bugs
, &oenv
);
665 bIndex
= opt2bSet("-n", NFILE
, fnm
) || bNDEF
;
666 bMead
= opt2bSet("-mead", NFILE
, fnm
);
667 bSetSize
= opt2parg_bSet("-box", NPA
, pa
);
668 bSetAng
= opt2parg_bSet("-angles", NPA
, pa
);
669 bSetCenter
= opt2parg_bSet("-center", NPA
, pa
);
670 bDist
= opt2parg_bSet("-d", NPA
, pa
);
671 bAlign
= opt2parg_bSet("-align", NPA
, pa
);
672 /* Only automatically turn on centering without -noc */
673 if ((bDist
|| bSetSize
|| bSetCenter
) && !opt2parg_bSet("-c", NPA
, pa
))
677 bScale
= opt2parg_bSet("-scale", NPA
, pa
);
678 bRho
= opt2parg_bSet("-density", NPA
, pa
);
679 bTranslate
= opt2parg_bSet("-translate", NPA
, pa
);
680 bRotate
= opt2parg_bSet("-rotate", NPA
, pa
);
682 fprintf(stderr
, "WARNING: setting -density overrides -scale\n");
683 bScale
= bScale
|| bRho
;
684 bCalcGeom
= bCenter
|| bRotate
|| bOrient
|| bScale
;
685 bCalcDiam
= btype
[0][0] == 'c' || btype
[0][0] == 'd' || btype
[0][0] == 'o';
687 infile
= ftp2fn(efSTX
, NFILE
, fnm
);
689 outfile
= ftp2fn(efPQR
, NFILE
, fnm
);
691 outfile
= ftp2fn(efSTO
, NFILE
, fnm
);
692 outftp
= fn2ftp(outfile
);
693 inftp
= fn2ftp(infile
);
695 aps
= gmx_atomprop_init();
699 printf("Incompatible options -mead and -grasp. Turning off -grasp\n");
702 if (bGrasp
&& (outftp
!= efPDB
))
703 gmx_fatal(FARGS
, "Output file should be a .pdb file"
704 " when using the -grasp option\n");
705 if ((bMead
|| bGrasp
) && !((fn2ftp(infile
) == efTPR
) ||
706 (fn2ftp(infile
) == efTPA
) ||
707 (fn2ftp(infile
) == efTPB
)))
708 gmx_fatal(FARGS
,"Input file should be a .tp[abr] file"
709 " when using the -mead option\n");
711 get_stx_coordnum(infile
,&natom
);
712 init_t_atoms(&atoms
,natom
,TRUE
);
715 read_stx_conf(infile
,title
,&atoms
,x
,v
,&ePBC
,box
);
716 if (fn2ftp(infile
) == efPDB
)
718 get_pdb_atomnumber(&atoms
,aps
);
720 printf("Read %d atoms\n",atoms
.nr
);
722 /* Get the element numbers if available in a pdb file */
723 if (fn2ftp(infile
) == efPDB
)
724 get_pdb_atomnumber(&atoms
,aps
);
726 if (ePBC
!= epbcNONE
)
729 printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
730 vol
,100*((int)(vol
*4.5)));
733 if (bMead
|| bGrasp
|| bCONECT
)
734 top
= read_top(infile
,NULL
);
738 if (atoms
.nr
!= top
->atoms
.nr
)
739 gmx_fatal(FARGS
,"Atom numbers don't match (%d vs. %d)",atoms
.nr
,top
->atoms
.nr
);
740 snew(atoms
.pdbinfo
,top
->atoms
.nr
);
741 ntype
= top
->idef
.atnr
;
742 for(i
=0; (i
<atoms
.nr
); i
++) {
743 /* Determine the Van der Waals radius from the force field */
745 if (!gmx_atomprop_query(aps
,epropVDW
,
746 *top
->atoms
.resinfo
[top
->atoms
.atom
[i
].resind
].name
,
747 *top
->atoms
.atomname
[i
],&vdw
))
751 itype
= top
->atoms
.atom
[i
].type
;
752 c12
= top
->idef
.iparams
[itype
*ntype
+itype
].lj
.c12
;
753 c6
= top
->idef
.iparams
[itype
*ntype
+itype
].lj
.c6
;
754 if ((c6
!= 0) && (c12
!= 0)) {
760 vdw
= 0.5*pow(sig6
,1.0/6.0);
765 /* Factor of 10 for nm -> Angstroms */
769 atoms
.pdbinfo
[i
].occup
= top
->atoms
.atom
[i
].q
;
770 atoms
.pdbinfo
[i
].bfac
= vdw
;
773 atoms
.pdbinfo
[i
].occup
= vdw
;
774 atoms
.pdbinfo
[i
].bfac
= top
->atoms
.atom
[i
].q
;
779 for (i
=0; (i
<natom
) && !bHaveV
; i
++)
780 for (j
=0; (j
<DIM
) && !bHaveV
; j
++)
781 bHaveV
=bHaveV
|| (v
[i
][j
]!=0);
782 printf("%selocities found\n",bHaveV
?"V":"No v");
786 gmx_fatal(FARGS
,"Sorry, can not visualize box with index groups");
788 gmx_fatal(FARGS
,"Sorry, can only visualize box with a pdb file");
789 } else if (visbox
[0] == -1)
790 visualize_images("images.pdb",ePBC
,box
);
794 rm_gropbc(&atoms
,x
,box
);
798 fprintf(stderr
,"\nSelect a group for determining the system size:\n");
799 get_index(&atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),
800 1,&ssize
,&sindex
,&sgrpname
);
805 diam
=calc_geom(ssize
,sindex
,x
,gc
,min
,max
,bCalcDiam
);
806 rvec_sub(max
, min
, size
);
807 printf(" system size :%7.3f%7.3f%7.3f (nm)\n",
808 size
[XX
], size
[YY
], size
[ZZ
]);
810 printf(" diameter :%7.3f (nm)\n",diam
);
811 printf(" center :%7.3f%7.3f%7.3f (nm)\n", gc
[XX
], gc
[YY
], gc
[ZZ
]);
812 printf(" box vectors :%7.3f%7.3f%7.3f (nm)\n",
813 norm(box
[XX
]), norm(box
[YY
]), norm(box
[ZZ
]));
814 printf(" box angles :%7.2f%7.2f%7.2f (degrees)\n",
815 norm2(box
[ZZ
])==0 ? 0 :
816 RAD2DEG
*acos(cos_angle_no_table(box
[YY
],box
[ZZ
])),
817 norm2(box
[ZZ
])==0 ? 0 :
818 RAD2DEG
*acos(cos_angle_no_table(box
[XX
],box
[ZZ
])),
819 norm2(box
[YY
])==0 ? 0 :
820 RAD2DEG
*acos(cos_angle_no_table(box
[XX
],box
[YY
])));
821 printf(" box volume :%7.2f (nm^3)\n",det(box
));
824 if (bRho
|| bOrient
|| bAlign
)
825 mass
= calc_mass(&atoms
,!fn2bTPX(infile
),aps
);
831 /* Get a group for principal component analysis */
832 fprintf(stderr
,"\nSelect group for the determining the orientation\n");
833 get_index(&atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),1,&isize
,&index
,&grpnames
);
835 /* Orient the principal axes along the coordinate axes */
836 orient_princ(&atoms
,isize
,index
,natom
,x
,bHaveV
? v
: NULL
, NULL
);
842 /* scale coordinates and box */
844 /* Compute scaling constant */
848 dens
= (mass
*AMU
)/(vol
*NANO
*NANO
*NANO
);
849 fprintf(stderr
,"Volume of input %g (nm^3)\n",vol
);
850 fprintf(stderr
,"Mass of input %g (a.m.u.)\n",mass
);
851 fprintf(stderr
,"Density of input %g (g/l)\n",dens
);
852 if (vol
==0 || mass
==0)
853 gmx_fatal(FARGS
,"Cannot scale density with "
854 "zero mass (%g) or volume (%g)\n",mass
,vol
);
856 scale
[XX
] = scale
[YY
] = scale
[ZZ
] = pow(dens
/rho
,1.0/3.0);
857 fprintf(stderr
,"Scaling all box vectors by %g\n",scale
[XX
]);
859 scale_conf(atoms
.nr
,x
,box
,scale
);
864 fprintf(stderr
,"\nSelect a group that you want to align:\n");
865 get_index(&atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),
866 1,&asize
,&aindex
,&agrpname
);
870 for (i
=0;i
<asize
;i
++)
873 printf("Aligning %d atoms (out of %d) to %g %g %g, center of rotation %g %g %g\n",asize
,natom
,
874 targetvec
[XX
],targetvec
[YY
],targetvec
[ZZ
],
875 aligncenter
[XX
],aligncenter
[YY
],aligncenter
[ZZ
]);
876 /*subtract out pivot point*/
877 for(i
=0; i
<asize
; i
++)
878 rvec_dec(x
[aindex
[i
]],aligncenter
);
879 /*now determine transform and rotate*/
881 principal_comp(asize
,aindex
,atoms
.atom
,x
, trans
,princd
);
883 unitv(targetvec
,targetvec
);
884 printf("Using %g %g %g as principal axis\n", trans
[0][2],trans
[1][2],trans
[2][2]);
885 tmpvec
[XX
]=trans
[0][2]; tmpvec
[YY
]=trans
[1][2]; tmpvec
[ZZ
]=trans
[2][2];
886 calc_rotmatrix(tmpvec
, targetvec
, rotmatrix
);
887 /* rotmatrix finished */
889 for (i
=0;i
<asize
;++i
)
891 mvmul(rotmatrix
,x
[aindex
[i
]],tmpvec
);
892 copy_rvec(tmpvec
,x
[aindex
[i
]]);
895 /*add pivot point back*/
896 for(i
=0; i
<asize
; i
++)
897 rvec_inc(x
[aindex
[i
]],aligncenter
);
904 fprintf(stderr
,"\nSelect a group that you want to translate:\n");
905 get_index(&atoms
,ftp2fn_null(efNDX
,NFILE
,fnm
),
906 1,&ssize
,&sindex
,&sgrpname
);
911 printf("Translating %d atoms (out of %d) by %g %g %g nm\n",ssize
,natom
,
912 translation
[XX
],translation
[YY
],translation
[ZZ
]);
914 for(i
=0; i
<ssize
; i
++)
915 rvec_inc(x
[sindex
[i
]],translation
);
918 for(i
=0; i
<natom
; i
++)
919 rvec_inc(x
[i
],translation
);
924 printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n",rotangles
[XX
],rotangles
[YY
],rotangles
[ZZ
]);
926 rotangles
[i
] *= DEG2RAD
;
927 rotate_conf(natom
,x
,v
,rotangles
[XX
],rotangles
[YY
],rotangles
[ZZ
]);
931 /* recalc geometrical center and max and min coordinates and size */
932 calc_geom(ssize
,sindex
,x
,gc
,min
,max
,FALSE
);
933 rvec_sub(max
, min
, size
);
934 if (bScale
|| bOrient
|| bRotate
)
935 printf("new system size : %6.3f %6.3f %6.3f\n",
936 size
[XX
],size
[YY
],size
[ZZ
]);
939 if (bSetSize
|| bDist
|| (btype
[0][0]=='t' && bSetAng
)) {
941 if (!(bSetSize
|| bDist
))
942 for (i
=0; i
<DIM
; i
++)
943 newbox
[i
] = norm(box
[i
]);
945 /* calculate new boxsize */
950 newbox
[i
] = size
[i
]+2*dist
;
952 box
[XX
][XX
] = newbox
[XX
];
953 box
[YY
][YY
] = newbox
[YY
];
954 box
[ZZ
][ZZ
] = newbox
[ZZ
];
956 matrix_convert(box
,newbox
,newang
);
966 if (btype
[0][0] == 'c')
969 else if (btype
[0][0] == 'd') {
974 box
[ZZ
][ZZ
] = d
*sqrt(2)/2;
978 box
[YY
][YY
] = d
*sqrt(2)*2/3;
980 box
[ZZ
][YY
] = d
*sqrt(2)/3;
981 box
[ZZ
][ZZ
] = d
*sqrt(6)/3;
987 /* calculate new coords for geometrical center */
989 calc_box_center(ecenterDEF
,box
,center
);
991 /* center molecule on 'center' */
993 center_conf(natom
,x
,center
,gc
);
997 calc_geom(ssize
,sindex
,x
, gc
, min
, max
, FALSE
);
998 printf("new center :%7.3f%7.3f%7.3f (nm)\n",gc
[XX
],gc
[YY
],gc
[ZZ
]);
1000 if (bOrient
|| bScale
|| bDist
|| bSetSize
) {
1001 printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n",
1002 norm(box
[XX
]), norm(box
[YY
]), norm(box
[ZZ
]));
1003 printf("new box angles :%7.2f%7.2f%7.2f (degrees)\n",
1004 norm2(box
[ZZ
])==0 ? 0 :
1005 RAD2DEG
*acos(cos_angle_no_table(box
[YY
],box
[ZZ
])),
1006 norm2(box
[ZZ
])==0 ? 0 :
1007 RAD2DEG
*acos(cos_angle_no_table(box
[XX
],box
[ZZ
])),
1008 norm2(box
[YY
])==0 ? 0 :
1009 RAD2DEG
*acos(cos_angle_no_table(box
[XX
],box
[YY
])));
1010 printf("new box volume :%7.2f (nm^3)\n",det(box
));
1013 if (check_box(epbcXYZ
,box
))
1014 printf("\nWARNING: %s\n",check_box(epbcXYZ
,box
));
1016 if (bDist
&& btype
[0][0]=='t')
1020 printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
1021 "distance from the solute to a box surface along the corresponding normal\n"
1022 "vector might be somewhat smaller than your specified value %f.\n"
1023 "You can check the actual value with g_mindist -pi\n",dist
);
1027 printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
1028 "If the molecule rotates the actual distance will be smaller. You might want\n"
1029 "to use a cubic box instead, or why not try a dodecahedron today?\n");
1032 if (bCONECT
&& (outftp
== efPDB
) && (inftp
== efTPR
))
1033 conect
= gmx_conect_generate(top
);
1038 fprintf(stderr
,"\nSelect a group for output:\n");
1039 get_index(&atoms
,opt2fn_null("-n",NFILE
,fnm
),
1040 1,&isize
,&index
,&grpname
);
1041 if (opt2bSet("-bf",NFILE
,fnm
))
1042 gmx_fatal(FARGS
,"combination not implemented: -bf -n or -bf -ndef");
1044 if (outftp
== efPDB
) {
1045 out
=ffopen(outfile
,"w");
1046 write_pdbfile_indexed(out
,title
,&atoms
,x
,ePBC
,box
,' ',1,isize
,index
,conect
,TRUE
);
1050 write_sto_conf_indexed(outfile
,title
,&atoms
,x
,bHaveV
?v
:NULL
,ePBC
,box
,
1055 if ((outftp
== efPDB
) || (outftp
== efPQR
)) {
1056 out
=ffopen(outfile
,"w");
1058 set_pdb_wide_format(TRUE
);
1059 fprintf(out
,"REMARK "
1060 "The B-factors in this file hold atomic radii\n");
1061 fprintf(out
,"REMARK "
1062 "The occupancy in this file hold atomic charges\n");
1065 fprintf(out
,"GRASP PDB FILE\nFORMAT NUMBER=1\n");
1066 fprintf(out
,"REMARK "
1067 "The B-factors in this file hold atomic charges\n");
1068 fprintf(out
,"REMARK "
1069 "The occupancy in this file hold atomic radii\n");
1071 else if (opt2bSet("-bf",NFILE
,fnm
)) {
1072 read_bfac(opt2fn("-bf",NFILE
,fnm
),&n_bfac
,&bfac
,&bfac_nr
);
1073 set_pdb_conf_bfac(atoms
.nr
,atoms
.nres
,&atoms
,
1074 n_bfac
,bfac
,bfac_nr
,peratom
);
1076 if (opt2parg_bSet("-label",NPA
,pa
)) {
1077 for(i
=0; (i
<atoms
.nr
); i
++)
1078 atoms
.resinfo
[atoms
.atom
[i
].resind
].chainid
=label
[0];
1080 write_pdbfile(out
,title
,&atoms
,x
,ePBC
,box
,' ',-1,conect
,TRUE
);
1082 pdb_legend(out
,atoms
.nr
,atoms
.nres
,&atoms
,x
);
1084 visualize_box(out
,bLegend
? atoms
.nr
+12 : atoms
.nr
,
1085 bLegend
? atoms
.nres
=12 : atoms
.nres
,box
,visbox
);
1089 write_sto_conf(outfile
,title
,&atoms
,x
,bHaveV
?v
:NULL
,ePBC
,box
);
1091 gmx_atomprop_destroy(aps
);
1093 do_view(oenv
,outfile
,NULL
);