src/contrib/hexamer.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2006, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, 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.
  11  *
  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.
  16  *
  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.
  21  *
  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.
  26  *
  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.
  34  *
  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.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <string.h>
  44 #include <ctype.h>
  45 #include "pdbio.h"
  46 #include "confio.h"
  47 #include "symtab.h"
  48 #include "smalloc.h"
  49 #include "symtab.h"
  50 #include "macros.h"
  51 #include "smalloc.h"
  52 #include "copyrite.h"
  53 #include "statutil.h"
  54 #include "string2.h"
  55 #include "strdb.h"
  56 #include "index.h"
  57 #include "vec.h"
  58 #include "typedefs.h"
  59 #include "gbutil.h"
  60 #include "strdb.h"
  61 #include "physics.h"
  62 #include "atomprop.h"
  63
  64 void copy_atom(t_symtab *tab,t_atoms *a1,int i1,t_atoms *a2,int i2,
  65                rvec xin[],rvec xout[],rvec vin[],rvec vout[])
  66 {
  67   a2->atom[i2]     = a1->atom[i1];
  68   a2->atomname[i2] = put_symtab(tab,*a1->atomname[i1]);
  69   a2->resname[a2->atom[i2].resnr] =
  70     put_symtab(tab,*a1->resname[a1->atom[i1].resnr]);
  71   copy_rvec(xin[i1],xout[i2]);
  72   copy_rvec(vin[i1],vout[i2]);
  73 }
  74
  75 static void rotate_x(int natom,rvec xin[],real angle,rvec xout[],
  76                      gmx_bool bZ,gmx_bool bUpsideDown,real dz)
  77 {
  78   int i;
  79   matrix mat;
  80
  81   angle *= DEG2RAD;
  82   clear_mat(mat);
  83   if (bZ) {
  84     mat[XX][XX] = cos(angle);
  85     mat[XX][YY] = sin(angle);
  86     mat[YY][XX] = -sin(angle);
  87     mat[YY][YY] = cos(angle);
  88     mat[ZZ][ZZ] = 1;
  89   }
  90   else {
  91     mat[XX][XX] = 1;
  92     mat[YY][YY] = cos(angle);
  93     mat[YY][ZZ] = sin(angle);
  94     mat[ZZ][YY] = -sin(angle);
  95     mat[ZZ][ZZ] = cos(angle);
  96   }
  97
  98   for(i=0; (i<natom); i++) {
  99     mvmul(mat,xin[i],xout[i]);
 100     if (bUpsideDown)
 101       xout[i][ZZ] *= -1;
 102     xout[i][ZZ] += dz;
 103   }
 104 }
 105
 106 static void prep_x(int natom,rvec x[],real rDist,real rAngleZ,real rAngleX)
 107 {
 108   int  i;
 109   rvec xcm;
 110   rvec *xx;
 111
 112   /* Center on Z-axis */
 113   clear_rvec(xcm);
 114   for(i=0; (i<natom); i++) {
 115     xcm[XX] += x[i][XX];
 116     xcm[YY] += x[i][YY];
 117     xcm[ZZ] += x[i][ZZ];
 118   }
 119   xcm[XX] /= natom;
 120   xcm[YY] /= natom;
 121   xcm[ZZ] /= natom;
 122   for(i=0; (i<natom); i++) {
 123     x[i][XX] -= xcm[XX];
 124     x[i][YY] -= xcm[YY];
 125     x[i][ZZ] -= xcm[ZZ];
 126   }
 127   if (rAngleZ != 0) {
 128     snew(xx,natom);
 129     rotate_x(natom,x,rAngleZ,xx,TRUE,FALSE,0);
 130     for(i=0; (i<natom); i++)
 131       copy_rvec(xx[i],x[i]);
 132     sfree(xx);
 133   }
 134   if (rAngleX != 0) {
 135     snew(xx,natom);
 136     rotate_x(natom,x,rAngleX,xx,FALSE,FALSE,0);
 137     for(i=0; (i<natom); i++)
 138       copy_rvec(xx[i],x[i]);
 139     sfree(xx);
 140   }
 141   if (rDist > 0) {
 142     for(i=0; (i<natom); i++)
 143       x[i][XX] += rDist;
 144   }
 145 }
 146
 147 int main(int argc, char *argv[])
 148 {
 149   t_symtab tab;
 150   static char *desc[] = {
 151     "[TT]hexamer[tt] takes a single input coordinate file and makes five symmetry",
 152     "related copies."
 153   };
 154 #define NPA asize(pa)
 155   t_filenm fnm[] = {
 156     { efSTX, "-f", NULL, ffREAD },
 157     { efPDB, "-o", NULL, ffWRITE }
 158   };
 159 #define NFILE asize(fnm)
 160   gmx_bool bCenter    = FALSE;
 161   gmx_bool bTrimer    = FALSE;
 162   gmx_bool bAlternate = FALSE;
 163   real rDist = 0,rAngleZ = 0,rAngleX = 0, alterz = 0;
 164   t_pargs pa[] = {
 165     { "-center",   FALSE, etBOOL,  {&bCenter},
 166       "Center molecule on Z-axis first" },
 167     { "-trimer",   FALSE, etBOOL,  {&bTrimer},
 168       "Make trimer rather than hexamer" },
 169     { "-alternate",FALSE, etBOOL,  {&bAlternate},
 170       "Turn every other molecule upside down" },
 171     { "-alterz",   FALSE, etREAL,  {&alterz},
 172       "Add this amount to Z-coordinate in every other molecule" },
 173     { "-radius",   FALSE, etREAL,  {&rDist},
 174       "Distance of protein axis from Z-axis (implies [TT]-center[tt])" },
 175     { "-anglez",   FALSE, etREAL,  {&rAngleZ},
 176       "Initial angle of rotation around Z-axis of protein" },
 177     { "-anglex",   FALSE, etREAL,  {&rAngleX},
 178       "Initial angle of rotation around X-axis of protein" }
 179   };
 180 #define NPA asize(pa)
 181   FILE    *fp;
 182   int     i,iout,now,natom;
 183   rvec    *xin,*vin,*xout;
 184   matrix  box;
 185   t_atoms atoms,aout;
 186   char    *infile,*outfile,title[256],buf[32];
 187
 188   CopyRight(stderr,argv[0]);
 189   parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa,
 190                     asize(desc),desc,0,NULL);
 191   bCenter = bCenter || (rDist > 0) || bAlternate;
 192
 193   infile  = ftp2fn(efSTX,NFILE,fnm);
 194   outfile = ftp2fn(efPDB,NFILE,fnm);
 195
 196   get_stx_coordnum(infile,&natom);
 197   init_t_atoms(&atoms,natom,TRUE);
 198   snew(xin,natom);
 199   snew(xout,natom);
 200   snew(vin,natom);
 201   read_stx_conf(infile,title,&atoms,xin,vin,box);
 202   printf("Read %d atoms\n",atoms.nr);
 203
 204   if (bCenter)
 205     prep_x(atoms.nr,xin,rDist,rAngleZ,rAngleX);
 206
 207   fp = ffopen(outfile,"w");
 208   for(i=0; (i<(bTrimer ? 3 : 6)); i++) {
 209     rotate_x(atoms.nr,xin,i*(bTrimer ? 120.0 : 60.0),xout,TRUE,
 210              bAlternate && ((i % 2) != 0),alterz*(((i % 2) == 0) ? 0 : 1));
 211     sprintf(buf,"Rotated %d degrees",i*(bTrimer ? 120 : 60));
 212     write_pdbfile(fp,buf,&atoms,xout,box,'A'+i,1+i);
 213   }
 214   ffclose(fp);
 215
 216   thanx(stderr);
 217
 218   return 0;
 219 }