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