?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Upped the version to 3.2.0
[gromacs.git] / src / tools / editconf.c
blob12aa6bff156ff2d18ca738629c5e393eebfbd2d2
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
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.
15
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.
20 *
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.
27 *
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.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * Green Red Orange Magenta Azure Cyan Skyblue
35 */
36 #include <math.h>
37 #include <string.h>
38 #include <ctype.h>
39 #include "pdbio.h"
40 #include "confio.h"
41 #include "symtab.h"
42 #include "smalloc.h"
43 #include "macros.h"
44 #include "copyrite.h"
45 #include "statutil.h"
46 #include "string2.h"
47 #include "strdb.h"
48 #include "index.h"
49 #include "vec.h"
50 #include "typedefs.h"
51 #include "gbutil.h"
52 #include "strdb.h"
53 #include "physics.h"
54 #include "atomprop.h"
55 #include "tpxio.h"
56 #include "pbc.h"
57 #include "princ.h"
58 #include "txtdump.h"
59 #include "viewit.h"
60
61 typedef struct {
62 char sanm[12];
63 int natm;
64 int nw;
65 char anm[6][12];
66 } t_simat;
67
68 typedef struct {
69 char reso[12];
70 char resn[12];
71 int nsatm;
72 t_simat sat[3];
73 } t_simlist;
74 static char *pdbtp[epdbNR]={"ATOM ","HETATM"};
75
76 real calc_mass(t_atoms *atoms,bool bGetMass,void *atomprop)
77 {
78 real tmass;
79 int i;
80
81 tmass = 0;
82 for(i=0; (i<atoms->nr); i++) {
83 if (bGetMass)
84 (void) query_atomprop(atomprop,epropMass,
85 *atoms->resname[atoms->atom[i].resnr],
86 *atoms->atomname[i],&(atoms->atom[i].m));
87 tmass += atoms->atom[i].m;
88 }
89
90 return tmass;
91 }
92
93 real calc_geom(int isize,atom_id *index,
94 rvec *x, rvec geom_center, rvec min, rvec max,bool bDiam)
95 {
96 real diam2,d;
97 char *grpnames;
98 int ii,i,j;
99
100 clear_rvec(geom_center);
101 diam2 = 0;
102 if (isize == 0) {
103 clear_rvec(min);
104 clear_rvec(max);
105 } else {
106 if (index)
107 ii=index[0];
108 else
109 ii=0;
110 for (j=0; j<DIM; j++)
111 min[j]=max[j]=x[ii][j];
112 for (i=0; i<isize; i++) {
113 if (index)
114 ii = index[i];
115 else
116 ii = i;
117 rvec_inc(geom_center,x[ii]);
118 for (j=0; j<DIM; j++) {
119 if (x[ii][j] < min[j]) min[j]=x[ii][j];
120 if (x[ii][j] > max[j]) max[j]=x[ii][j];
121 }
122 if (bDiam) {
123 if (index)
124 for (j=i+1; j<isize; j++) {
125 d = distance2(x[ii],x[index[j]]);
126 diam2 = max(d,diam2);
127 }
128 else
129 for (j=i+1; j<isize; j++) {
130 d = distance2(x[i],x[j]);
131 diam2 = max(d,diam2);
132 }
133 }
134 }
135 svmul(1.0/isize,geom_center,geom_center);
136 }
137
138 return sqrt(diam2);
139 }
140
141 void center_conf(int natom, rvec *x, rvec center, rvec geom_cent)
142 {
143 int i;
144 rvec shift;
145
146 rvec_sub(center,geom_cent,shift);
147
148 printf(" shift :%7.3f%7.3f%7.3f (nm)\n",
149 shift[XX],shift[YY],shift[ZZ]);
150
151 for (i=0; (i<natom); i++)
152 rvec_inc(x[i], shift);
153 }
154
155 void scale_conf(int natom,rvec x[],matrix box,rvec scale)
156 {
157 int i,j;
158
159 for(i=0; i<natom; i++) {
160 for (j=0; j<DIM; j++)
161 x[i][j] *= scale[j];
162 }
163 for (i=0; i<DIM; i++)
164 for (j=0; j<DIM; j++)
165 box[i][j] *= scale[j];
166 }
167
168 void rm_gropbc(t_atoms *atoms,rvec x[],matrix box)
169 {
170 real dist;
171 int n,m,d;
172
173 /* check periodic boundary */
174 for(n=1;(n<atoms->nr);n++) {
175 for(m=DIM-1; m>=0; m--) {
176 dist = x[n][m]-x[n-1][m];
177 if (fabs(dist) > 0.9*box[m][m]) {
178 if ( dist > 0 )
179 for(d=0; d<=m; d++)
180 x[n][d] -= box[m][d];
181 else
182 for(d=0; d<=m; d++)
183 x[n][d] += box[m][d];
184 }
185 }
186 }
187 }
188
189 void read_bfac(char *fn, int *n_bfac, double **bfac_val, int **bfac_nr)
190 {
191 int i;
192 char **bfac_lines;
193
194 *n_bfac = get_lines(fn, &bfac_lines);
195 snew(*bfac_val, *n_bfac);
196 snew(*bfac_nr, *n_bfac);
197 fprintf(stderr, "Reading %d B-factors from %s\n",*n_bfac,fn);
198 for(i=0; (i<*n_bfac); i++) {
199 /*fprintf(stderr, "Line %d: %s",i,bfac_lines[i]);*/
200 sscanf(bfac_lines[i],"%d %lf",&(*bfac_nr)[i],&(*bfac_val)[i]);
201 /*fprintf(stderr," nr %d val %g\n",(*bfac_nr)[i],(*bfac_val)[i]);*/
202 }
203
204 }
205
206 void set_pdb_conf_bfac(int natoms,int nres,t_atoms *atoms,
207 int n_bfac,double *bfac,int *bfac_nr,
208 bool peratom)
209 {
210 FILE *out;
211 real bfac_min,bfac_max;
212 int i,n;
213 bool found;
214 char buf[120];
215
216 bfac_max=-1e10;
217 bfac_min=1e10;
218 for(i=0; (i<n_bfac); i++) {
219 if (bfac_nr[i]-1>=atoms->nres)
220 peratom=TRUE;
221 if ((bfac_nr[i]-1<0) || (bfac_nr[i]-1>=atoms->nr))
222 fatal_error(0,"Index of B-Factor %d is out of range: %d (%g)",
223 i+1,bfac_nr[i],bfac[i]);
224 if (bfac[i] > bfac_max)
225 bfac_max = bfac[i];
226 if (bfac[i] < bfac_min)
227 bfac_min = bfac[i];
228 }
229 while ( (bfac_max > 99.99) || (bfac_min < -99.99) ) {
230 fprintf(stderr,"Range of values for B-factors too large (min %g, max %g) "
231 "will scale down a factor 10\n",bfac_min,bfac_max);
232 for(i=0; (i<n_bfac); i++)
233 bfac[i] /= 10;
234 bfac_max /= 10;
235 bfac_min /= 10;
236 }
237 while ( (fabs(bfac_max) < 0.5) && (fabs(bfac_min) < 0.5) ) {
238 fprintf(stderr,"Range of values for B-factors too small (min %g, max %g) "
239 "will scale up a factor 10\n",bfac_min,bfac_max);
240 for(i=0; (i<n_bfac); i++)
241 bfac[i] *= 10;
242 bfac_max *= 10;
243 bfac_min *= 10;
244 }
245
246 for(i=0; (i<natoms); i++)
247 atoms->pdbinfo[i].bfac=0;
248
249 if (!peratom) {
250 fprintf(stderr,"Will attach %d B-factors to %d residues\n",
251 n_bfac,nres);
252 for(i=0; (i<n_bfac); i++) {
253 found=FALSE;
254 for(n=0; (n<natoms); n++)
255 if ( bfac_nr[i] == (atoms->atom[n].resnr+1) ) {
256 atoms->pdbinfo[n].bfac=bfac[i];
257 found=TRUE;
258 }
259 if (!found) {
260 sprintf(buf,"Residue nr %d not found\n",bfac_nr[i]);
261 warning(buf);
262 }
263 }
264 } else {
265 fprintf(stderr,"Will attach %d B-factors to %d atoms\n",n_bfac,natoms);
266 for(i=0; (i<n_bfac); i++) {
267 atoms->pdbinfo[bfac_nr[i]-1].bfac=bfac[i];
268 }
269 }
270 }
271
272 void pdb_legend(FILE *out,int natoms,int nres,t_atoms *atoms,rvec x[])
273 {
274 real bfac_min,bfac_max,xmin,ymin,zmin;
275 int i;
276
277 bfac_max=-1e10;
278 bfac_min=1e10;
279 xmin = 1e10;
280 ymin = 1e10;
281 zmin = 1e10;
282 for (i=0; (i<natoms); i++) {
283 xmin = min(xmin,x[i][XX]);
284 ymin = min(ymin,x[i][YY]);
285 zmin = min(zmin,x[i][ZZ]);
286 bfac_min = min(bfac_min,atoms->pdbinfo[i].bfac);
287 bfac_max = max(bfac_max,atoms->pdbinfo[i].bfac);
288 }
289 fprintf(stderr,"B-factors range from %g to %g\n",bfac_min,bfac_max);
290 for (i=1; (i<12); i++) {
291 fprintf(out,pdbformat,
292 "ATOM ",natoms+1+i,"CA","LEG",' ',nres+1,
293 (xmin+(i*0.12))*10,ymin*10,zmin*10,1.0,
294 bfac_min+ ((i-1.0)*(bfac_max-bfac_min)/10) );
295 }
296 }
297 void visualize_images(char *fn,matrix box)
298 {
299 t_atoms atoms;
300 rvec *img;
301 char *c,*ala;
302 int nat,i;
303
304 nat = NTRICIMG+1;
305 init_t_atoms(&atoms,nat,FALSE);
306 atoms.nr = nat;
307 snew(img,nat);
308 c = "C";
309 ala = "ALA";
310 for(i=0; i<nat; i++) {
311 atoms.atomname[i] = &c;
312 atoms.atom[i].resnr = i;
313 atoms.resname[i] = &ala;
314 atoms.atom[i].chain = 'A'+i/NCUCVERT;
315 }
316 calc_triclinic_images(box,img+1);
317
318 write_sto_conf(fn,"Images",&atoms,img,NULL,box);
319
320 free_t_atoms(&atoms);
321 sfree(img);
322 }
323
324
325 void visualize_box(FILE *out,int a0,int r0,matrix box,rvec gridsize)
326 {
327 int *edge;
328 rvec *vert,shift;
329 int nx,ny,nz,nbox,nat;
330 int i,j,x,y,z;
331 int rectedge[24] = { 0,1, 1,3, 3,2, 0,2, 0,4, 1,5, 3,7, 2,6, 4,5, 5,7, 7,6, 6,4 };
332
333 a0++;
334 r0++;
335
336 nx = (int)(gridsize[XX]+0.5);
337 ny = (int)(gridsize[YY]+0.5);
338 nz = (int)(gridsize[ZZ]+0.5);
339 nbox = nx*ny*nz;
340 if (TRICLINIC(box)) {
341 nat = nbox*NCUCVERT;
342 snew(vert,nat);
343 calc_compact_unitcell_vertices(box,vert);
344 j = 0;
345 for(z=0; z<nz; z++)
346 for(y=0; y<ny; y++)
347 for(x=0; x<nx; x++) {
348 for(i=0; i<DIM; i++)
349 shift[i] = x*box[0][i]+y*box[1][i]+z*box[2][i];
350 for(i=0; i<NCUCVERT; i++) {
351 rvec_add(vert[i],shift,vert[j]);
352 j++;
353 }
354 }
355
356 for(i=0; i<nat; i++) {
357 fprintf(out,pdbformat,"ATOM",a0+i,"C","BOX",'K'+i/NCUCVERT,r0+i,
358 10*vert[i][XX],10*vert[i][YY],10*vert[i][ZZ]);
359 fprintf(out,"\n");
360 }
361
362 edge = compact_unitcell_edges();
363 for(j=0; j<nbox; j++)
364 for(i=0; i<NCUCEDGE; i++)
365 fprintf(out,"CONECT%5d%5d\n",
366 a0 + j*NCUCVERT + edge[2*i],
367 a0 + j*NCUCVERT + edge[2*i+1]);
368
369 sfree(vert);
370 } else {
371 i=0;
372 for(z=0; z<=1; z++)
373 for(y=0; y<=1; y++)
374 for(x=0; x<=1; x++) {
375 fprintf(out,pdbformat,"ATOM",a0+i,"C","BOX",'K'+i/8,r0+i,
376 x*10*box[XX][XX],y*10*box[YY][YY],z*10*box[ZZ][ZZ]);
377 fprintf(out,"\n");
378 i++;
379 }
380 for(i=0; i<24; i+=2)
381 fprintf(out,"CONECT%5d%5d\n",a0+rectedge[i],a0+rectedge[i+1]);
382 }
383 }
384
385 int main(int argc, char *argv[])
386 {
387 static char *desc[] = {
388 "editconf converts generic structure format to [TT].gro[tt], [TT].g96[tt]",
389 "or [TT].pdb[tt].",
390 "[PAR]",
391 "The box can be modified with options [TT]-box[tt], [TT]-d[tt] and",
392 "[TT]-angles[tt]. Both [TT]-box[tt] and [TT]-d[tt]",
393 "will center the system in the box.",
394 "[PAR]",
395 "Option [TT]-bt[tt] determines the box type: [TT]tric[tt] is a",
396 "triclinic box, [TT]cubic[tt] is a cubic box, [TT]dodecahedron[tt] is",
397 "a rhombic dodecahedron and [TT]octahedron[tt] is a truncated octahedron.",
398 "The last two are special cases of a triclinic box.",
399 "The length of the three box vectors of the truncated octahedron is the",
400 "shortest distance between two opposite hexagons.",
401 "The volume of a dodecahedron is 0.71 and that of a truncated octahedron",
402 "is 0.77 of that of a cubic box with the same periodic image distance.",
403 "[PAR]",
404 "Option [TT]-box[tt] requires only",
405 "one value for a cubic box, dodecahedron and a truncated octahedron.",
406 "With [TT]-d[tt] and [TT]tric[tt] the size of the system in the x, y",
407 "and z directions is used. With [TT]-d[tt] and [TT]cubic[tt],",
408 "[TT]dodecahedron[tt] or [TT]octahedron[tt] the diameter of the system",
409 "is used, which is the largest distance between two atoms.",
410 "[PAR]",
411 "Option [TT]-angles[tt] is only meaningful with option [TT]-box[tt] and",
412 "a triclinic box and can not be used with option [TT]-d[tt].",
413 "[PAR]",
414 "When [TT]-n[tt] or [TT]-ndef[tt] is set, a group",
415 "can be selected for calculating the size and the geometric center,",
416 "otherwise the whole system is used.",
417 "[PAR]",
418 "[TT]-rotate[tt] rotates the coordinates and velocities.",
419 "[TT]-princ[tt] aligns the principal axes of the system along the",
420 "coordinate axes, this may allow you to decrease the box volume,",
421 "but beware that molecules can rotate significantly in a nanosecond.[PAR]",
422 "Scaling is applied before any of the other operations are",
423 "performed. Boxes can be scaled to give a certain density (option",
424 "[TT]-density[tt]). A special feature of the scaling option, when the",
425 "factor -1 is given in one dimension, one obtains a mirror image,",
426 "mirrored in one of the plains, when one uses -1 in three dimensions",
427 "a point-mirror image is obtained.[PAR]",
428 "Groups are selected after all operations have been applied.[PAR]",
429 "Periodicity can be removed in a crude manner.",
430 "It is important that the box sizes at the bottom of your input file",
431 "are correct when the periodicity is to be removed.",
432 "[PAR]",
433 "The program can optionally rotate the solute molecule to align the",
434 "molecule along its principal axes ([TT]-rotate[tt])",
435 "[PAR]",
436 "When writing [TT].pdb[tt] files, B-factors can be",
437 "added with the [TT]-bf[tt] option. B-factors are read",
438 "from a file with with following format: first line states number of",
439 "entries in the file, next lines state an index",
440 "followed by a B-factor. The B-factors will be attached per residue",
441 "unless an index is larger than the number of residues or unless the",
442 "[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
443 "be added instead of B-factors. [TT]-legend[tt] will produce",
444 "a row of CA atoms with B-factors ranging from the minimum to the",
445 "maximum value found, effectively making a legend for viewing.",
446 "[PAR]",
447 "With the option -mead a special pdb file for the MEAD electrostatics",
448 "program (Poisson-Boltzmann solver) can be made. A further prerequisite",
449 "is that the input file is a run input file.",
450 "The B-factor field is then filled with the Van der Waals radius",
451 "of the atoms while the occupancy field will hold the charge.",
452 "[PAR]",
453 "The option -grasp is similar, but it puts the charges in the B-factor",
454 "and the radius in the occupancy.",
455 "[PAR]",
456 "Finally with option [TT]-label[tt] editconf can add a chain identifier",
457 "to a pdb file, which can be useful for analysis with e.g. rasmol."
458 "[PAR]",
459 "To convert a truncated octrahedron file produced by a package which uses",
460 "a cubic box with the corners cut off (such as Gromos) use:[BR]",
461 "[TT]editconf -f <in> -rotate 0 -45 -35.264 -bt o -box <veclen> -o <out>[tt][BR]",
462 "where [TT]veclen[tt] is the size of the cubic box times sqrt(3)/2."
463 };
464 static char *bugs[] = {
465 "For complex molecules, the periodicity removal routine may break down, "
466 "in that case you can use trjconv"
467 };
468 static real dist=0.0,rbox=0.0,to_diam=0.0;
469 static bool bNDEF=FALSE,bRMPBC=FALSE,bCenter=FALSE,bVOL=TRUE;
470 static bool peratom=FALSE,bLegend=FALSE,bOrient=FALSE,bMead=FALSE,bGrasp=FALSE;
471 static rvec scale={1,1,1},newbox={0,0,0},newang={90,90,90};
472 static real rho=1000.0,rvdw=0.12;
473 static rvec center={0,0,0},translation={0,0,0},rotangles={0,0,0};
474 static char *btype[]={ NULL, "tric", "cubic", "dodecahedron", "octahedron", NULL },*label="A";
475 static rvec visbox={0,0,0};
476 t_pargs pa[] = {
477 { "-ndef", FALSE, etBOOL, {&bNDEF},
478 "Choose output from default index groups" },
479 { "-visbox", FALSE, etRVEC, {visbox},
480 "HIDDENVisualize a grid of boxes, -1 visualizes the 14 box images" },
481 { "-bt", FALSE, etENUM, {btype},
482 "Box type for -box and -d" },
483 { "-box", FALSE, etRVEC, {newbox}, "Box vector lengths (a,b,c)" },
484 { "-angles", FALSE, etRVEC, {newang},
485 "Angles between the box vectors (bc,ac,ab)" },
486 { "-d", FALSE, etREAL, {&dist},
487 "Distance between the solute and the box" },
488 { "-c", FALSE, etBOOL, {&bCenter},
489 "Center molecule in box (implied by -box and -d)" },
490 { "-center", FALSE, etRVEC, {center}, "Coordinates of geometrical center"},
491 { "-translate", FALSE, etRVEC, {translation},
492 "Translation" },
493 { "-rotate", FALSE, etRVEC, {rotangles},
494 "Rotation around the X, Y and Z axes in degrees" },
495 { "-princ", FALSE, etBOOL, {&bOrient}, "Orient molecule(s) along their principal axes" },
496 { "-scale", FALSE, etRVEC, {scale}, "Scaling factor" },
497 { "-density",FALSE, etREAL, {&rho},
498 "Density (g/l) of the output box achieved by scaling" },
499 { "-vol", FALSE, etBOOL, {&bVOL},
500 "Compute and print volume of the box" },
501 { "-pbc", FALSE, etBOOL, {&bRMPBC},
502 "Remove the periodicity (make molecule whole again)" },
503 { "-mead", FALSE, etBOOL, {&bMead},
504 "Store the charge of the atom in the occupancy field and the radius of the atom in the B-factor field" },
505 { "-grasp", FALSE, etBOOL, {&bGrasp},
506 "Store the charge of the atom in the B-factor field and the radius of the atom in the occupancy field" },
507 { "-rvdw", FALSE, etREAL, {&rvdw},
508 "Default Van der Waals radius if one can not be found in the database" },
509 { "-atom", FALSE, etBOOL, {&peratom}, "Force B-factor attachment per atom" },
510 { "-legend", FALSE, etBOOL, {&bLegend}, "Make B-factor legend" },
511 { "-label", FALSE, etSTR, {&label}, "Add chain label for all residues" }
512 };
513 #define NPA asize(pa)
514
515 FILE *out;
516 void *atomprop;
517 char *infile,*outfile,title[STRLEN];
518 int outftp,natom,i,j,n_bfac;
519 double *bfac=NULL;
520 int *bfac_nr=NULL;
521 t_topology *top;
522 t_atoms atoms;
523 char *grpname,*sgrpname;
524 int isize,ssize;
525 atom_id *index,*sindex;
526 rvec *x,*v,gc,min,max,size;
527 matrix box;
528 bool bIndex,bSetSize,bSetAng,bCubic,bDist,bSetCenter;
529 bool bHaveV,bScale,bRho,bTranslate,bRotate,bCalcGeom,bCalcDiam;
530 real xs,ys,zs,xcent,ycent,zcent,diam=0,mass=0,d,vdw;
531 t_filenm fnm[] = {
532 { efSTX, "-f", NULL, ffREAD },
533 { efNDX, "-n", NULL, ffOPTRD },
534 { efSTO, NULL, NULL, ffWRITE },
535 { efDAT, "-bf", "bfact", ffOPTRD }
536 };
537 #define NFILE asize(fnm)
538
539 CopyRight(stderr,argv[0]);
540 parse_common_args(&argc,argv,PCA_CAN_VIEW,NFILE,fnm,NPA,pa,
541 asize(desc),desc,asize(bugs),bugs);
542
543 bIndex = opt2bSet("-n",NFILE,fnm) || bNDEF;
544 bSetSize = opt2parg_bSet("-box" ,NPA,pa);
545 bSetAng = opt2parg_bSet("-angles" ,NPA,pa);
546 bSetCenter= opt2parg_bSet("-center" ,NPA,pa);
547 bDist = opt2parg_bSet("-d" ,NPA,pa);
548 bCenter = bCenter || bDist || bSetCenter || bSetSize;
549 bScale = opt2parg_bSet("-scale" ,NPA,pa);
550 bRho = opt2parg_bSet("-density",NPA,pa);
551 bTranslate= opt2parg_bSet("-translate",NPA,pa);
552 bRotate = opt2parg_bSet("-rotate",NPA,pa);
553 if (bScale && bRho)
554 fprintf(stderr,"WARNING: setting -density overrides -scale\n");
555 bScale = bScale || bRho;
556 bCalcGeom = bCenter || bTranslate || bRotate || bOrient || bScale;
557 bCalcDiam = btype[0][0]=='c' || btype[0][0]=='d' || btype[0][0]=='o';
558
559 infile = ftp2fn(efSTX,NFILE,fnm);
560 outfile = ftp2fn(efSTO,NFILE,fnm);
561 outftp = fn2ftp(outfile);
562 atomprop = get_atomprop();
563 if (bMead && bGrasp) {
564 fprintf(stderr,"Incompatible options -mead and -grasp. Turning off -grasp\n");
565 bGrasp = FALSE;
566 }
567 if ((bMead || bGrasp) && (outftp != efPDB))
568 fatal_error(0,"Output file should be a .pdb file"
569 " when using the -mead option\n");
570 if ((bMead || bGrasp) && !((fn2ftp(infile) == efTPR) ||
571 (fn2ftp(infile) == efTPA) ||
572 (fn2ftp(infile) == efTPB)))
573 fatal_error(0,"Input file should be a .tp[abr] file"
574 " when using the -mead option\n");
575
576 get_stx_coordnum(infile,&natom);
577 init_t_atoms(&atoms,natom,TRUE);
578 snew(x,natom);
579 snew(v,natom);
580 read_stx_conf(infile,title,&atoms,x,v,box);
581 printf("Read %d atoms\n",atoms.nr);
582 if (bVOL) {
583 real vol = det(box);
584 printf("Volume: %g nm^3, corresponds to roughly %d electrons\n",
585 vol,100*((int)(vol*4.5)));
586 }
587
588 if (bMead || bGrasp) {
589 top = read_top(infile);
590 if (atoms.nr != top->atoms.nr)
591 fatal_error(0,"Atom numbers don't match (%d vs. %d)",
592 atoms.nr,top->atoms.nr);
593 for(i=0; (i<atoms.nr); i++) {
594 /* Factor of 10 for Angstroms */
595 if (query_atomprop(atomprop,epropVDW,
596 *top->atoms.resname[top->atoms.atom[i].resnr],
597 *top->atoms.atomname[i],&vdw))
598 vdw = 10*vdw;
599 else
600 vdw = 10*rvdw;
601 if (bMead) {
602 atoms.pdbinfo[i].occup = top->atoms.atom[i].q;
603 atoms.pdbinfo[i].bfac = vdw;
604 }
605 else {
606 /* Factor of 10 for Angstroms */
607 atoms.pdbinfo[i].occup = vdw;
608 atoms.pdbinfo[i].bfac = top->atoms.atom[i].q;
609 }
610 }
611 }
612 bHaveV=FALSE;
613 for (i=0; (i<natom) && !bHaveV; i++)
614 for (j=0; (j<DIM) && !bHaveV; j++)
615 bHaveV=bHaveV || (v[i][j]!=0);
616 printf("%selocities found\n",bHaveV?"V":"No v");
617
618 if (visbox[0] > 0) {
619 if (bIndex)
620 fatal_error(0,"Sorry, can not visualize box with index groups");
621 if (outftp != efPDB)
622 fatal_error(0,"Sorry, can only visualize box with a pdb file");
623 } else if (visbox[0] == -1)
624 visualize_images("images.pdb",box);
625
626 /* remove pbc */
627 if (bRMPBC)
628 rm_gropbc(&atoms,x,box);
629
630 if (bCalcGeom) {
631 if (bIndex) {
632 fprintf(stderr,"\nSelect a group for determining the system size:\n");
633 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),
634 1,&ssize,&sindex,&sgrpname);
635 } else {
636 ssize = atoms.nr;
637 sindex = NULL;
638 }
639 diam=calc_geom(ssize,sindex,x,gc,min,max,bCalcDiam);
640 rvec_sub(max, min, size);
641 printf(" system size :%7.3f%7.3f%7.3f (nm)\n",
642 size[XX], size[YY], size[ZZ]);
643 if (bCalcDiam)
644 printf(" diameter :%7.3f (nm)\n",diam);
645 printf(" center :%7.3f%7.3f%7.3f (nm)\n", gc[XX], gc[YY], gc[ZZ]);
646 printf(" box vectors :%7.3f%7.3f%7.3f (nm)\n",
647 norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
648 printf(" box angles :%7.2f%7.2f%7.2f (degrees)\n",
649 norm2(box[ZZ])==0 ? 0 :
650 RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
651 norm2(box[ZZ])==0 ? 0 :
652 RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
653 norm2(box[YY])==0 ? 0 :
654 RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
655 printf(" box volume :%7.2f (nm^3)\n",det(box));
656 }
657
658 if (bRho || bOrient)
659 mass = calc_mass(&atoms,!fn2bTPX(infile),atomprop);
660
661 if (bOrient) {
662 atom_id *index;
663 char *grpnames;
664
665 /* Get a group for principal component analysis */
666 fprintf(stderr,"\nSelect group for orientation of molecule:\n");
667 get_index(&atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpnames);
668
669 /* Orient the principal axes along the coordinate axes */
670 orient_princ(&atoms,isize,index,natom,x,bHaveV ? v : NULL, NULL);
671 sfree(index);
672 sfree(grpnames);
673 }
674
675 if ( bScale ) {
676 /* scale coordinates and box */
677 if (bRho) {
678 /* Compute scaling constant */
679 real vol,dens;
680
681 vol = det(box);
682 dens = (mass*AMU)/(vol*NANO*NANO*NANO);
683 fprintf(stderr,"Volume of input %g (nm^3)\n",vol);
684 fprintf(stderr,"Mass of input %g (a.m.u.)\n",mass);
685 fprintf(stderr,"Density of input %g (g/l)\n",dens);
686 if (vol==0 || mass==0)
687 fatal_error(0,"Cannot scale density with "
688 "zero mass (%g) or volume (%g)\n",mass,vol);
689
690 scale[XX] = scale[YY] = scale[ZZ] = pow(dens/rho,1.0/3.0);
691 fprintf(stderr,"Scaling all box vectors by %g\n",scale[XX]);
692 }
693 scale_conf(atoms.nr,x,box,scale);
694 }
695
696 if (bTranslate) {
697 printf("Translating by %g %g %g nm\n",translation[XX],translation[YY],translation[ZZ]);
698 for(i=0; i<natom; i++)
699 rvec_inc(x[i],translation);
700 }
701 if (bRotate) {
702 /* Rotate */
703 printf("Rotating %g, %g, %g degrees around the X, Y and Z axis respectively\n",rotangles[XX],rotangles[YY],rotangles[ZZ]);
704 for(i=0; i<DIM; i++)
705 rotangles[i] *= DEG2RAD;
706 rotate_conf(natom,x,v,rotangles[XX],rotangles[YY],rotangles[ZZ]);
707 }
708
709 if (bCalcGeom) {
710 /* recalc geometrical center and max and min coordinates and size */
711 calc_geom(ssize,sindex,x,gc,min,max,FALSE);
712 rvec_sub(max, min, size);
713 if (bScale || bOrient || bRotate)
714 printf("new system size : %6.3f %6.3f %6.3f\n",
715 size[XX],size[YY],size[ZZ]);
716 }
717
718 if (bSetSize || bDist || (btype[0][0]=='t' && bSetAng)) {
719 if (!(bSetSize || bDist))
720 for (i=0; i<DIM; i++)
721 newbox[i] = norm(box[i]);
722 clear_mat(box);
723 /* calculate new boxsize */
724 switch(btype[0][0]){
725 case 't':
726 if (bDist)
727 for(i=0; i<DIM; i++)
728 newbox[i] = size[i]+2*dist;
729 if (!bSetAng) {
730 box[XX][XX] = newbox[XX];
731 box[YY][YY] = newbox[YY];
732 box[ZZ][ZZ] = newbox[ZZ];
733 } else {
734 svmul(DEG2RAD,newang,newang);
735 box[XX][XX] = newbox[XX];
736 box[YY][XX] = newbox[YY]*cos(newang[ZZ]);
737 box[YY][YY] = newbox[YY]*sin(newang[ZZ]);
738 box[ZZ][XX] = newbox[ZZ]*cos(newang[YY]);
739 box[ZZ][YY] = newbox[ZZ]
740 *(cos(newang[XX])-cos(newang[YY])*cos(newang[ZZ]))/sin(newang[ZZ]);
741 box[ZZ][ZZ] = sqrt(sqr(newbox[ZZ])
742 -box[ZZ][XX]*box[ZZ][XX]-box[ZZ][YY]*box[ZZ][YY]);
743 }
744 if (bDist && TRICLINIC(box))
745 fprintf(stderr,"WARNING: the box is triclinic, the minimum distance between the solute and the box might be less than %f\nYou can check this with g_mindist -pi\n",dist);
746 break;
747 case 'c':
748 case 'd':
749 case 'o':
750 if (bSetSize)
751 d = newbox[0];
752 else
753 d = diam+2*dist;
754 if (btype[0][0] == 'c')
755 for(i=0; i<DIM; i++)
756 box[i][i] = d;
757 else if (btype[0][0] == 'd') {
758 box[XX][XX] = d;
759 box[YY][YY] = d;
760 box[ZZ][XX] = d/2;
761 box[ZZ][YY] = d/2;
762 box[ZZ][ZZ] = d*sqrt(2)/2;
763 } else {
764 box[XX][XX] = d;
765 box[YY][XX] = d/3;
766 box[YY][YY] = d*sqrt(2)*2/3;
767 box[ZZ][XX] = -d/3;
768 box[ZZ][YY] = d*sqrt(2)/3;
769 box[ZZ][ZZ] = d*sqrt(6)/3;
770 }
771 break;
772 }
773 }
774
775 /* calculate new coords for geometrical center */
776 if (!bSetCenter)
777 calc_box_center(box,center);
778
779 /* center molecule on 'center' */
780 if (bCenter)
781 center_conf(natom,x,center,gc);
782
783 /* print some */
784 if (bCalcGeom) {
785 calc_geom(ssize,sindex,x, gc, min, max, FALSE);
786 printf("new center :%7.3f%7.3f%7.3f (nm)\n",gc[XX],gc[YY],gc[ZZ]);
787 }
788 if (bOrient || bScale || bDist || bSetSize) {
789 printf("new box vectors :%7.3f%7.3f%7.3f (nm)\n",
790 norm(box[XX]), norm(box[YY]), norm(box[ZZ]));
791 printf("new box angles :%7.2f%7.2f%7.2f (degrees)\n",
792 norm2(box[ZZ])==0 ? 0 :
793 RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ])),
794 norm2(box[ZZ])==0 ? 0 :
795 RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
796 norm2(box[YY])==0 ? 0 :
797 RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
798 printf("new box volume :%7.2f (nm^3)\n",det(box));
799 }
800
801 if (check_box(box))
802 printf("\nWARNING: %s\n",check_box(box));
803
804 if (bIndex) {
805 fprintf(stderr,"\nSelect a group for output:\n");
806 get_index(&atoms,opt2fn_null("-n",NFILE,fnm),
807 1,&isize,&index,&grpname);
808 if (opt2bSet("-bf",NFILE,fnm))
809 fatal_error(0,"combination not implemented: -bf -n or -bf -ndef");
810 else
811 write_sto_conf_indexed(outfile,title,&atoms,x,bHaveV?v:NULL,box,
812 isize,index);
813 }
814 else {
815 if (outftp != efPDB) {
816 write_sto_conf(outfile,title,&atoms,x,bHaveV?v:NULL,box);
817 }
818 else {
819 out=ffopen(outfile,"w");
820 if (opt2bSet("-bf",NFILE,fnm)) {
821 read_bfac(opt2fn("-bf",NFILE,fnm),&n_bfac,&bfac,&bfac_nr);
822 set_pdb_conf_bfac(atoms.nr,atoms.nres,&atoms,
823 n_bfac,bfac,bfac_nr,peratom);
824 }
825 if (opt2parg_bSet("-label",NPA,pa)) {
826 for(i=0; (i<atoms.nr); i++)
827 atoms.atom[i].chain=label[0];
828 }
829 if (bMead || bGrasp) {
830 if (bMead) {
831 set_pdb_wide_format(TRUE);
832 fprintf(out,"REMARK "
833 "The b-factors in this file hold atomic radii\n");
834 fprintf(out,"REMARK "
835 "The occupancy in this file hold atomic charges\n");
836 }
837 else {
838 fprintf(out,"GRASP PDB FILE\nFORMAT NUMBER=1\n");
839 fprintf(out,"REMARK "
840 "The b-factors in this file hold atomic charges\n");
841 fprintf(out,"REMARK "
842 "The occupancy in this file hold atomic radii\n");
843 }
844 }
845 write_pdbfile(out,title,&atoms,x,box,0,-1);
846 if (bLegend)
847 pdb_legend(out,atoms.nr,atoms.nres,&atoms,x);
848 if (visbox[0] > 0)
849 visualize_box(out,bLegend ? atoms.nr+12 : atoms.nr,
850 bLegend? atoms.nres=12 : atoms.nres,box,visbox);
851 fclose(out);
852 }
853 }
854 done_atomprop(&atomprop);
855
856 do_view(outfile,NULL);
857
858 thanx(stderr);
859
860 return 0;
861 }
862
The ultimate molecular dynamics simulation package