src/contrib/hexamer.c

   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_mkyaw_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include <ctype.h>
  34 #include "assert.h"
  35 #include "pdbio.h"
  36 #include "confio.h"
  37 #include "symtab.h"
  38 #include "smalloc.h"
  39 #include "symtab.h"
  40 #include "macros.h"
  41 #include "smalloc.h"
  42 #include "copyrite.h"
  43 #include "statutil.h"
  44 #include "string2.h"
  45 #include "strdb.h"
  46 #include "rdgroup.h"
  47 #include "vec.h"
  48 #include "typedefs.h"
  49 #include "gbutil.h"
  50 #include "strdb.h"
  51 #include "rdgroup.h"
  52 #include "physics.h"
  53 #include "atomprop.h"
  54
  55 void copy_atom(t_symtab *tab,t_atoms *a1,int i1,t_atoms *a2,int i2,
  56                rvec xin[],rvec xout[],rvec vin[],rvec vout[])
  57 {
  58   a2->atom[i2]     = a1->atom[i1];
  59   a2->atomname[i2] = put_symtab(tab,*a1->atomname[i1]);
  60   a2->resname[a2->atom[i2].resnr] =
  61     put_symtab(tab,*a1->resname[a1->atom[i1].resnr]);
  62   copy_rvec(xin[i1],xout[i2]);
  63   copy_rvec(vin[i1],vout[i2]);
  64 }
  65
  66 static void rotate_x(int natom,rvec xin[],real angle,rvec xout[],
  67                      bool bZ,bool bUpsideDown,real dz)
  68 {
  69   int i;
  70   matrix mat;
  71
  72   angle *= DEG2RAD;
  73   clear_mat(mat);
  74   if (bZ) {
  75     mat[XX][XX] = cos(angle);
  76     mat[XX][YY] = sin(angle);
  77     mat[YY][XX] = -sin(angle);
  78     mat[YY][YY] = cos(angle);
  79     mat[ZZ][ZZ] = 1;
  80   }
  81   else {
  82     mat[XX][XX] = 1;
  83     mat[YY][YY] = cos(angle);
  84     mat[YY][ZZ] = sin(angle);
  85     mat[ZZ][YY] = -sin(angle);
  86     mat[ZZ][ZZ] = cos(angle);
  87   }
  88
  89   for(i=0; (i<natom); i++) {
  90     mvmul(mat,xin[i],xout[i]);
  91     if (bUpsideDown)
  92       xout[i][ZZ] *= -1;
  93     xout[i][ZZ] += dz;
  94   }
  95 }
  96
  97 static void prep_x(int natom,rvec x[],real rDist,real rAngleZ,real rAngleX)
  98 {
  99   int  i;
 100   rvec xcm;
 101   rvec *xx;
 102
 103   /* Center on Z-axis */
 104   clear_rvec(xcm);
 105   for(i=0; (i<natom); i++) {
 106     xcm[XX] += x[i][XX];
 107     xcm[YY] += x[i][YY];
 108     xcm[ZZ] += x[i][ZZ];
 109   }
 110   xcm[XX] /= natom;
 111   xcm[YY] /= natom;
 112   xcm[ZZ] /= natom;
 113   for(i=0; (i<natom); i++) {
 114     x[i][XX] -= xcm[XX];
 115     x[i][YY] -= xcm[YY];
 116     x[i][ZZ] -= xcm[ZZ];
 117   }
 118   if (rAngleZ != 0) {
 119     snew(xx,natom);
 120     rotate_x(natom,x,rAngleZ,xx,TRUE,FALSE,0);
 121     for(i=0; (i<natom); i++)
 122       copy_rvec(xx[i],x[i]);
 123     sfree(xx);
 124   }
 125   if (rAngleX != 0) {
 126     snew(xx,natom);
 127     rotate_x(natom,x,rAngleX,xx,FALSE,FALSE,0);
 128     for(i=0; (i<natom); i++)
 129       copy_rvec(xx[i],x[i]);
 130     sfree(xx);
 131   }
 132   if (rDist > 0) {
 133     for(i=0; (i<natom); i++)
 134       x[i][XX] += rDist;
 135   }
 136 }
 137
 138 int main(int argc, char *argv[])
 139 {
 140   t_symtab tab;
 141   static char *desc[] = {
 142     "hexamer takes a single input coordinate file and makes five symmetry",
 143     "related copies."
 144   };
 145 #define NPA asize(pa)
 146   t_filenm fnm[] = {
 147     { efSTX, "-f", NULL, ffREAD },
 148     { efPDB, "-o", NULL, ffWRITE }
 149   };
 150 #define NFILE asize(fnm)
 151   bool bCenter    = FALSE;
 152   bool bTrimer    = FALSE;
 153   bool bAlternate = FALSE;
 154   real rDist = 0,rAngleZ = 0,rAngleX = 0, alterz = 0;
 155   t_pargs pa[] = {
 156     { "-center",   FALSE, etBOOL,  {&bCenter},
 157       "Center molecule on Z-axis first" },
 158     { "-trimer",   FALSE, etBOOL,  {&bTrimer},
 159       "Make trimer rather than hexamer" },
 160     { "-alternate",FALSE, etBOOL,  {&bAlternate},
 161       "Turn every other molecule upside down" },
 162     { "-alterz",   FALSE, etREAL,  {&alterz},
 163       "Add this amount to Z-coordinate in every other molecule" },
 164     { "-radius",   FALSE, etREAL,  {&rDist},
 165       "Distance of protein axis from Z-axis (implies -center)" },
 166     { "-anglez",   FALSE, etREAL,  {&rAngleZ},
 167       "Initial angle of rotation around Z-axis of protein" },
 168     { "-anglex",   FALSE, etREAL,  {&rAngleX},
 169       "Initial angle of rotation around X-axis of protein" }
 170   };
 171 #define NPA asize(pa)
 172   FILE    *fp;
 173   int     i,iout,now,natom;
 174   rvec    *xin,*vin,*xout;
 175   matrix  box;
 176   t_atoms atoms,aout;
 177   char    *infile,*outfile,title[256],buf[32];
 178
 179   CopyRight(stderr,argv[0]);
 180   parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,NPA,pa,
 181                     asize(desc),desc,0,NULL);
 182   bCenter = bCenter || (rDist > 0) || bAlternate;
 183
 184   infile  = ftp2fn(efSTX,NFILE,fnm);
 185   outfile = ftp2fn(efPDB,NFILE,fnm);
 186
 187   get_stx_coordnum(infile,&natom);
 188   init_t_atoms(&atoms,natom,TRUE);
 189   snew(xin,natom);
 190   snew(xout,natom);
 191   snew(vin,natom);
 192   read_stx_conf(infile,title,&atoms,xin,vin,box);
 193   printf("Read %d atoms\n",atoms.nr);
 194
 195   if (bCenter)
 196     prep_x(atoms.nr,xin,rDist,rAngleZ,rAngleX);
 197
 198   fp = ffopen(outfile,"w");
 199   for(i=0; (i<(bTrimer ? 3 : 6)); i++) {
 200     rotate_x(atoms.nr,xin,i*(bTrimer ? 120.0 : 60.0),xout,TRUE,
 201              bAlternate && ((i % 2) != 0),alterz*(((i % 2) == 0) ? 0 : 1));
 202     sprintf(buf,"Rotated %d degrees",i*(bTrimer ? 120 : 60));
 203     write_pdbfile(fp,buf,&atoms,xout,box,'A'+i,1+i);
 204   }
 205   fclose(fp);
 206
 207   thanx(stderr);
 208
 209   return 0;
 210 }