src/contrib/hexamer.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.3.99_development_20071104
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2006, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Groningen Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <string.h>
  41
  42 #include "gromacs/fileio/pdbio.h"
  43 #include "gromacs/fileio/confio.h"
  44 #include "symtab.h"
  45 #include "gromacs/utility/smalloc.h"
  46 #include "symtab.h"
  47 #include "macros.h"
  48 #include "copyrite.h"
  49 #include "gromacs/commandline/pargs.h"
  50 #include "index.h"
  51 #include "vec.h"
  52 #include "typedefs.h"
  53 #include "gbutil.h"
  54 #include "physics.h"
  55 #include "atomprop.h"
  56
  57 void copy_atom(t_symtab *tab,t_atoms *a1,int i1,t_atoms *a2,int i2,
  58                rvec xin[],rvec xout[],rvec vin[],rvec vout[])
  59 {
  60   a2->atom[i2]     = a1->atom[i1];
  61   a2->atomname[i2] = put_symtab(tab,*a1->atomname[i1]);
  62   a2->resname[a2->atom[i2].resnr] =
  63     put_symtab(tab,*a1->resname[a1->atom[i1].resnr]);
  64   copy_rvec(xin[i1],xout[i2]);
  65   copy_rvec(vin[i1],vout[i2]);
  66 }
  67
  68 static void rotate_x(int natom,rvec xin[],real angle,rvec xout[],
  69                      gmx_bool bZ,gmx_bool bUpsideDown,real dz)
  70 {
  71   int i;
  72   matrix mat;
  73
  74   angle *= DEG2RAD;
  75   clear_mat(mat);
  76   if (bZ) {
  77     mat[XX][XX] = cos(angle);
  78     mat[XX][YY] = sin(angle);
  79     mat[YY][XX] = -sin(angle);
  80     mat[YY][YY] = cos(angle);
  81     mat[ZZ][ZZ] = 1;
  82   }
  83   else {
  84     mat[XX][XX] = 1;
  85     mat[YY][YY] = cos(angle);
  86     mat[YY][ZZ] = sin(angle);
  87     mat[ZZ][YY] = -sin(angle);
  88     mat[ZZ][ZZ] = cos(angle);
  89   }
  90
  91   for(i=0; (i<natom); i++) {
  92     mvmul(mat,xin[i],xout[i]);
  93     if (bUpsideDown)
  94       xout[i][ZZ] *= -1;
  95     xout[i][ZZ] += dz;
  96   }
  97 }
  98
  99 static void prep_x(int natom,rvec x[],real rDist,real rAngleZ,real rAngleX)
 100 {
 101   int  i;
 102   rvec xcm;
 103   rvec *xx;
 104
 105   /* Center on Z-axis */
 106   clear_rvec(xcm);
 107   for(i=0; (i<natom); i++) {
 108     xcm[XX] += x[i][XX];
 109     xcm[YY] += x[i][YY];
 110     xcm[ZZ] += x[i][ZZ];
 111   }
 112   xcm[XX] /= natom;
 113   xcm[YY] /= natom;
 114   xcm[ZZ] /= natom;
 115   for(i=0; (i<natom); i++) {
 116     x[i][XX] -= xcm[XX];
 117     x[i][YY] -= xcm[YY];
 118     x[i][ZZ] -= xcm[ZZ];
 119   }
 120   if (rAngleZ != 0) {
 121     snew(xx,natom);
 122     rotate_x(natom,x,rAngleZ,xx,TRUE,FALSE,0);
 123     for(i=0; (i<natom); i++)
 124       copy_rvec(xx[i],x[i]);
 125     sfree(xx);
 126   }
 127   if (rAngleX != 0) {
 128     snew(xx,natom);
 129     rotate_x(natom,x,rAngleX,xx,FALSE,FALSE,0);
 130     for(i=0; (i<natom); i++)
 131       copy_rvec(xx[i],x[i]);
 132     sfree(xx);
 133   }
 134   if (rDist > 0) {
 135     for(i=0; (i<natom); i++)
 136       x[i][XX] += rDist;
 137   }
 138 }
 139
 140 int main(int argc, char *argv[])
 141 {
 142   t_symtab tab;
 143   static char *desc[] = {
 144     "[TT]hexamer[tt] takes a single input coordinate file and makes five symmetry",
 145     "related copies."
 146   };
 147 #define NPA asize(pa)
 148   t_filenm fnm[] = {
 149     { efSTX, "-f", NULL, ffREAD },
 150     { efPDB, "-o", NULL, ffWRITE }
 151   };
 152 #define NFILE asize(fnm)
 153   gmx_bool bCenter    = FALSE;
 154   gmx_bool bTrimer    = FALSE;
 155   gmx_bool bAlternate = FALSE;
 156   real rDist = 0,rAngleZ = 0,rAngleX = 0, alterz = 0;
 157   t_pargs pa[] = {
 158     { "-center",   FALSE, etBOOL,  {&bCenter},
 159       "Center molecule on Z-axis first" },
 160     { "-trimer",   FALSE, etBOOL,  {&bTrimer},
 161       "Make trimer rather than hexamer" },
 162     { "-alternate",FALSE, etBOOL,  {&bAlternate},
 163       "Turn every other molecule upside down" },
 164     { "-alterz",   FALSE, etREAL,  {&alterz},
 165       "Add this amount to Z-coordinate in every other molecule" },
 166     { "-radius",   FALSE, etREAL,  {&rDist},
 167       "Distance of protein axis from Z-axis (implies [TT]-center[tt])" },
 168     { "-anglez",   FALSE, etREAL,  {&rAngleZ},
 169       "Initial angle of rotation around Z-axis of protein" },
 170     { "-anglex",   FALSE, etREAL,  {&rAngleX},
 171       "Initial angle of rotation around X-axis of protein" }
 172   };
 173 #define NPA asize(pa)
 174   FILE    *fp;
 175   int     i,iout,now,natom;
 176   rvec    *xin,*vin,*xout;
 177   matrix  box;
 178   t_atoms atoms,aout;
 179   char    *infile,*outfile,title[256],buf[32];
 180
 181   CopyRight(stderr,argv[0]);
 182   parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa,
 183                     asize(desc),desc,0,NULL);
 184   bCenter = bCenter || (rDist > 0) || bAlternate;
 185
 186   infile  = ftp2fn(efSTX,NFILE,fnm);
 187   outfile = ftp2fn(efPDB,NFILE,fnm);
 188
 189   get_stx_coordnum(infile,&natom);
 190   init_t_atoms(&atoms,natom,TRUE);
 191   snew(xin,natom);
 192   snew(xout,natom);
 193   snew(vin,natom);
 194   read_stx_conf(infile,title,&atoms,xin,vin,box);
 195   printf("Read %d atoms\n",atoms.nr);
 196
 197   if (bCenter)
 198     prep_x(atoms.nr,xin,rDist,rAngleZ,rAngleX);
 199
 200   fp = gmx_ffopen(outfile,"w");
 201   for(i=0; (i<(bTrimer ? 3 : 6)); i++) {
 202     rotate_x(atoms.nr,xin,i*(bTrimer ? 120.0 : 60.0),xout,TRUE,
 203              bAlternate && ((i % 2) != 0),alterz*(((i % 2) == 0) ? 0 : 1));
 204     sprintf(buf,"Rotated %d degrees",i*(bTrimer ? 120 : 60));
 205     write_pdbfile(fp,buf,&atoms,xout,box,'A'+i,1+i);
 206   }
 207   gmx_ffclose(fp);
 208
 209   gmx_thanx(stderr);
 210
 211   return 0;
 212 }