src/gmxlib/gbutil.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-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, 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 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <math.h>
  44 #include "macros.h"
  45 #include "vec.h"
  46 #include "gmx_fatal.h"
  47 #include "gstat.h"
  48 #include "pbc.h"
  49 #include "gbutil.h"
  50
  51 static real dist2(t_pbc *pbc,rvec x,rvec y)
  52 {
  53   rvec dx;
  54
  55   pbc_dx(pbc,x,y,dx);
  56
  57   return norm2(dx);
  58 }
  59
  60 real distance_to_z(rvec x)
  61 {
  62   return (sqr(x[XX])+sqr(x[YY]));
  63 } /*distance_to_z()*/
  64
  65 static void low_rotate_conf(int natom,rvec *x,real alfa, real beta,real gamma)
  66 {
  67   int  i;
  68   rvec x_old;
  69
  70   for (i=0; i<natom; i++) {
  71     copy_rvec(x[i],x_old);
  72     /*calculate new x[i] by rotation alfa around the x-axis*/
  73     x[i][XX]=   x_old[XX];
  74     x[i][YY]=             cos(alfa)*x_old[YY] - sin(alfa)*x_old[ZZ];
  75     x[i][ZZ]=             sin(alfa)*x_old[YY] + cos(alfa)*x_old[ZZ];
  76     copy_rvec(x[i],x_old);
  77     /*calculate new x[i] by rotation beta around the y-axis*/
  78     x[i][XX]=   cos(beta)*x_old[XX]           + sin(beta)*x_old[ZZ];
  79     x[i][YY]=                       x_old[YY];
  80     x[i][ZZ]= - sin(beta)*x_old[XX]           + cos(beta)*x_old[ZZ];
  81     copy_rvec(x[i],x_old);
  82     /*calculate new x[i] by rotation gamma around the z-axis*/
  83     x[i][XX]= x_old[XX]*cos(gamma) - x_old[YY]*sin(gamma);
  84     x[i][YY]= x_old[XX]*sin(gamma) + x_old[YY]*cos(gamma);
  85     x[i][ZZ]=                                             x_old[ZZ];
  86   }
  87 }
  88
  89 static void low_rotate_conf_indexed(int nindex,atom_id *index,rvec *x,real alfa, real beta,real gamma)
  90 {
  91   int  i;
  92   rvec x_old;
  93
  94   for (i=0; i<nindex; i++) {
  95     copy_rvec(x[index[i]],x_old);
  96     /*calculate new x[index[i]] by rotation alfa around the x-axis*/
  97     x[index[i]][XX]=   x_old[XX];
  98     x[index[i]][YY]=             cos(alfa)*x_old[YY] - sin(alfa)*x_old[ZZ];
  99     x[index[i]][ZZ]=             sin(alfa)*x_old[YY] + cos(alfa)*x_old[ZZ];
 100     copy_rvec(x[index[i]],x_old);
 101     /*calculate new x[index[i]] by rotation beta around the y-axis*/
 102     x[index[i]][XX]=   cos(beta)*x_old[XX]           + sin(beta)*x_old[ZZ];
 103     x[index[i]][YY]=                       x_old[YY];
 104     x[index[i]][ZZ]= - sin(beta)*x_old[XX]           + cos(beta)*x_old[ZZ];
 105     copy_rvec(x[index[i]],x_old);
 106     /*calculate new x[index[i]] by rotation gamma around the z-axis*/
 107     x[index[i]][XX]= x_old[XX]*cos(gamma) - x_old[YY]*sin(gamma);
 108     x[index[i]][YY]= x_old[XX]*sin(gamma) + x_old[YY]*cos(gamma);
 109     x[index[i]][ZZ]=                                             x_old[ZZ];
 110   }
 111 }
 112
 113 void rotate_conf(int natom,rvec *x,rvec *v,real alfa, real beta,real gamma)
 114 {
 115   if (x)
 116     low_rotate_conf(natom,x,alfa,beta,gamma);
 117   if (v)
 118     low_rotate_conf(natom,v,alfa,beta,gamma);
 119 }
 120
 121 void orient(int natom,rvec *x,rvec *v, rvec angle,matrix box)
 122 {
 123   real longest,rij,rzi;
 124   int  i,j,m,max_i=0,max_j=0;
 125   rvec origin;
 126   int  temp;
 127   real alfa=0,beta=0,gamma=0;
 128   t_pbc pbc;
 129
 130   set_pbc(&pbc,-1,box);
 131
 132   /*first i am going to look for the longest atom-atom distance*/
 133   longest=dist2(&pbc,x[0],x[1]);
 134   i=0;
 135   j=1;
 136   for (i=0;(i<natom);i++) {
 137     for (j=0;(j<natom);j++) {
 138       rij=dist2(&pbc,x[i],x[j]);
 139       if (rij>longest) {
 140         max_i=i;
 141         max_j=j;
 142         longest=rij;
 143       }
 144     }
 145   }
 146   /* first check if x[max_i]<x[max_j] else swap*/
 147   if (x[max_i][2]>x[max_j][2]) {
 148     temp=max_i;
 149     max_i=max_j;
 150     max_j=temp;
 151   }
 152
 153   /*set the origin to x[i]*/
 154   for(m=0;(m<DIM);m++)
 155     origin[m]=x[max_i][m];
 156   for(i=0;(i<natom);i++)
 157     for(m=0;(m<DIM);m++)
 158       x[i][m]-=origin[m];
 159
 160   /* calculate the rotation angles alfa(x_axis) and beta(y_axis)
 161    * the rotation angles must be calculated clockwise looking
 162    * along the rotation axis to the origin*
 163    * alfa (x-axis)
 164    */
 165   alfa=atan(x[max_j][ZZ]/x[max_j][YY])-M_PI_2;
 166   beta=M_PI_2-atan(x[max_j][ZZ]/x[max_j][XX]);
 167   rotate_conf(natom,x,v,alfa,beta,gamma);
 168
 169   /* now search the longest distance for rotation along the z_axis */
 170   longest=distance_to_z(x[0]);
 171   max_i=0;
 172   for (i=1;(i<natom);i++) {
 173     rzi=distance_to_z(x[i]);
 174     if (rzi>longest) {
 175       longest = rzi;
 176       max_i=i;
 177     }
 178   }
 179   gamma=atan(x[max_i][YY]/x[max_i][XX])-M_PI_2;
 180   rotate_conf(natom,x,v,0,0,gamma);
 181   angle[0]=alfa;
 182   angle[1]=beta;
 183   angle[2]=gamma;
 184 } /*orient()*/
 185
 186
 187 void genconf(t_atoms *atoms,rvec *x,rvec *v,real *r,matrix box,ivec n_box)
 188 {
 189   int     i,ix,iy,iz,m,j,imol,offset;
 190   rvec    delta;
 191   int     nmol;
 192
 193   nmol=n_box[XX]*n_box[YY]*n_box[ZZ];
 194
 195   /*print message*/
 196   fprintf(stderr,"Generating configuration\n");
 197   imol=0;
 198   for(ix=0; (ix < n_box[XX]); ix++) {
 199     delta[XX]=ix*box[XX][XX];
 200     for(iy=0; (iy < n_box[YY]); iy++) {
 201       delta[YY]=iy*box[YY][YY];
 202       for(iz=0; (iz < n_box[ZZ]); iz++) {
 203         delta[ZZ]=iz*box[ZZ][ZZ];
 204         offset=imol*atoms->nr;
 205         for (i=0;(i < atoms->nr);i++) {
 206           for (m=0;(m < DIM);m++)
 207             x[offset+i][m]=delta[m]+x[i][m];
 208           if (v)
 209             for (m=0;(m < DIM);m++)
 210               v[offset+i][m]=v[i][m];
 211           r[offset+i]=r[i];
 212         }
 213         imol++;
 214       }
 215     }
 216   }
 217   for (i=1;(i<nmol);i++) {
 218     int offs    = i*atoms->nr;
 219     int offsres = i*atoms->nres;
 220     for (j=0;(j<atoms->nr);j++) {
 221       atoms->atomname[offs+j]  = atoms->atomname[j];
 222       atoms->atom[offs+j].resind = atoms->atom[j].resind + offsres;
 223       atoms->resinfo[atoms->atom[offs+j].resind] =
 224         atoms->resinfo[atoms->atom[j].resind];
 225       atoms->resinfo[atoms->atom[offs+j].resind].nr += offsres;
 226     }
 227   }
 228   atoms->nr*=nmol;
 229   atoms->nres*=nmol;
 230   for(i=0; i<DIM; i++)
 231     for(j=0; j<DIM; j++)
 232       box[j][i]*=n_box[j];
 233 } /*genconf()*/
 234
 235 /*gen_box() generates a box around a configuration*/
 236 void gen_box(int NTB,int natoms,rvec *x, matrix box,rvec box_space,
 237              gmx_bool bCenter)
 238 {
 239   int i,m;
 240   rvec xmin, xmax;
 241   real max_box;
 242
 243   /*calculate minimum and maximum x[0..DIM-1]*/
 244   for (m=0;(m<DIM);m++)
 245     xmin[m]=xmax[m]=x[0][m];
 246   for (i=1;(i < natoms); i++)
 247     for (m=0;m<DIM;m++) {
 248       xmin[m]=min(xmin[m],x[i][m]);
 249       xmax[m]=max(xmax[m],x[i][m]);
 250     }
 251
 252   /*calculate the new box sizes for cubic and octahedral ...*/
 253   for (m=0; (m<DIM);m++)
 254     box[m][m]=xmax[m]-xmin[m]+2*box_space[m];
 255
 256   /*calculate the box size if NTB=1 (truncated octahedron)*/
 257   if (NTB==1) {
 258     max_box=box[0][0];
 259     for(m=0;(m<DIM);m++)
 260       max_box=max(max_box,box[m][m]);
 261     for (m=0;(m<DIM);m++)
 262       box[m][m]=max_box;
 263   }
 264
 265   /*move the molecule to the center of the box*/
 266   if (bCenter)
 267     for(i=0;(i<natoms);i++)
 268       for (m=0;(m<DIM);m++) {
 269         x[i][m]+=0.5*(box[m][m]-xmin[m]-xmax[m]);
 270       }
 271
 272
 273 #ifdef DEBUG
 274   /* print data to check this */
 275   print_stat(x,natoms,box);
 276 #endif
 277 }/*gen_box()*/
 278