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 orient(int natom,rvec *x,rvec *v, rvec angle,matrix box)
 119 {
 120   real longest,rij,rzi;
 121   int  i,j,m,max_i=0,max_j=0;
 122   rvec origin;
 123   int  temp;
 124   real alfa=0,beta=0,gamma=0;
 125   t_pbc pbc;
 126
 127   set_pbc(&pbc,-1,box);
 128
 129   /*first i am going to look for the longest atom-atom distance*/
 130   longest=dist2(&pbc,x[0],x[1]);
 131   i=0;
 132   j=1;
 133   for (i=0;(i<natom);i++) {
 134     for (j=0;(j<natom);j++) {
 135       rij=dist2(&pbc,x[i],x[j]);
 136       if (rij>longest) {
 137         max_i=i;
 138         max_j=j;
 139         longest=rij;
 140       }
 141     }
 142   }
 143   /* first check if x[max_i]<x[max_j] else swap*/
 144   if (x[max_i][2]>x[max_j][2]) {
 145     temp=max_i;
 146     max_i=max_j;
 147     max_j=temp;
 148   }
 149
 150   /*set the origin to x[i]*/
 151   for(m=0;(m<DIM);m++)
 152     origin[m]=x[max_i][m];
 153   for(i=0;(i<natom);i++)
 154     for(m=0;(m<DIM);m++)
 155       x[i][m]-=origin[m];
 156
 157   /* calculate the rotation angles alfa(x_axis) and beta(y_axis)
 158    * the rotation angles must be calculated clockwise looking
 159    * along the rotation axis to the origin*
 160    * alfa (x-axis)
 161    */
 162   alfa=atan(x[max_j][ZZ]/x[max_j][YY])-M_PI_2;
 163   beta=M_PI_2-atan(x[max_j][ZZ]/x[max_j][XX]);
 164   rotate_conf(natom,x,v,alfa,beta,gamma);
 165
 166   /* now search the longest distance for rotation along the z_axis */
 167   longest=distance_to_z(x[0]);
 168   max_i=0;
 169   for (i=1;(i<natom);i++) {
 170     rzi=distance_to_z(x[i]);
 171     if (rzi>longest) {
 172       longest = rzi;
 173       max_i=i;
 174     }
 175   }
 176   gamma=atan(x[max_i][YY]/x[max_i][XX])-M_PI_2;
 177   rotate_conf(natom,x,v,0,0,gamma);
 178   angle[0]=alfa;
 179   angle[1]=beta;
 180   angle[2]=gamma;
 181 } /*orient()*/
 182
 183
 184 void genconf(t_atoms *atoms,rvec *x,rvec *v,real *r,matrix box,ivec n_box)
 185 {
 186   int     i,ix,iy,iz,m,j,imol,offset;
 187   rvec    delta;
 188   int     nmol;
 189
 190   nmol=n_box[XX]*n_box[YY]*n_box[ZZ];
 191
 192   /*print message*/
 193   fprintf(stderr,"Generating configuration\n");
 194   imol=0;
 195   for(ix=0; (ix < n_box[XX]); ix++) {
 196     delta[XX]=ix*box[XX][XX];
 197     for(iy=0; (iy < n_box[YY]); iy++) {
 198       delta[YY]=iy*box[YY][YY];
 199       for(iz=0; (iz < n_box[ZZ]); iz++) {
 200         delta[ZZ]=iz*box[ZZ][ZZ];
 201         offset=imol*atoms->nr;
 202         for (i=0;(i < atoms->nr);i++) {
 203           for (m=0;(m < DIM);m++)
 204             x[offset+i][m]=delta[m]+x[i][m];
 205           if (v)
 206             for (m=0;(m < DIM);m++)
 207               v[offset+i][m]=v[i][m];
 208           r[offset+i]=r[i];
 209         }
 210         imol++;
 211       }
 212     }
 213   }
 214   for (i=1;(i<nmol);i++) {
 215     int offs    = i*atoms->nr;
 216     int offsres = i*atoms->nres;
 217     for (j=0;(j<atoms->nr);j++) {
 218       atoms->atomname[offs+j]  = atoms->atomname[j];
 219       atoms->atom[offs+j].resind = atoms->atom[j].resind + offsres;
 220       atoms->resinfo[atoms->atom[offs+j].resind] =
 221         atoms->resinfo[atoms->atom[j].resind];
 222       atoms->resinfo[atoms->atom[offs+j].resind].nr += offsres;
 223     }
 224   }
 225   atoms->nr*=nmol;
 226   atoms->nres*=nmol;
 227   for(i=0; i<DIM; i++)
 228     for(j=0; j<DIM; j++)
 229       box[j][i]*=n_box[j];
 230 } /*genconf()*/
 231
 232 /*gen_box() generates a box around a configuration*/
 233 void gen_box(int NTB,int natoms,rvec *x, matrix box,rvec box_space,
 234              gmx_bool bCenter)
 235 {
 236   int i,m;
 237   rvec xmin, xmax;
 238   real max_box;
 239
 240   /*calculate minimum and maximum x[0..DIM-1]*/
 241   for (m=0;(m<DIM);m++)
 242     xmin[m]=xmax[m]=x[0][m];
 243   for (i=1;(i < natoms); i++)
 244     for (m=0;m<DIM;m++) {
 245       xmin[m]=min(xmin[m],x[i][m]);
 246       xmax[m]=max(xmax[m],x[i][m]);
 247     }
 248
 249   /*calculate the new box sizes for cubic and octahedral ...*/
 250   for (m=0; (m<DIM);m++)
 251     box[m][m]=xmax[m]-xmin[m]+2*box_space[m];
 252
 253   /*calculate the box size if NTB=1 (truncated octahedron)*/
 254   if (NTB==1) {
 255     max_box=box[0][0];
 256     for(m=0;(m<DIM);m++)
 257       max_box=max(max_box,box[m][m]);
 258     for (m=0;(m<DIM);m++)
 259       box[m][m]=max_box;
 260   }
 261
 262   /*move the molecule to the center of the box*/
 263   if (bCenter)
 264     for(i=0;(i<natoms);i++)
 265       for (m=0;(m<DIM);m++) {
 266         x[i][m]+=0.5*(box[m][m]-xmin[m]-xmax[m]);
 267       }
 268
 269
 270 #ifdef DEBUG
 271   /* print data to check this */
 272   print_stat(x,natoms,box);
 273 #endif
 274 }/*gen_box()*/
 275