src/gromacs/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     {
  69         copy_rvec(x[i], x_old);
  70         /*calculate new x[i] by rotation alfa around the x-axis*/
  71         x[i][XX] =   x_old[XX];
  72         x[i][YY] =             cos(alfa)*x_old[YY] - sin(alfa)*x_old[ZZ];
  73         x[i][ZZ] =             sin(alfa)*x_old[YY] + cos(alfa)*x_old[ZZ];
  74         copy_rvec(x[i], x_old);
  75         /*calculate new x[i] by rotation beta around the y-axis*/
  76         x[i][XX] =   cos(beta)*x_old[XX]           + sin(beta)*x_old[ZZ];
  77         x[i][YY] =                       x_old[YY];
  78         x[i][ZZ] = -sin(beta)*x_old[XX]           + cos(beta)*x_old[ZZ];
  79         copy_rvec(x[i], x_old);
  80         /*calculate new x[i] by rotation gamma around the z-axis*/
  81         x[i][XX] = x_old[XX]*cos(gamma) - x_old[YY]*sin(gamma);
  82         x[i][YY] = x_old[XX]*sin(gamma) + x_old[YY]*cos(gamma);
  83         x[i][ZZ] =                                             x_old[ZZ];
  84     }
  85 }
  86
  87 static void low_rotate_conf_indexed(int nindex, atom_id *index, rvec *x, real alfa, real beta, real gamma)
  88 {
  89     int  i;
  90     rvec x_old;
  91
  92     for (i = 0; i < nindex; i++)
  93     {
  94         copy_rvec(x[index[i]], x_old);
  95         /*calculate new x[index[i]] by rotation alfa around the x-axis*/
  96         x[index[i]][XX] =   x_old[XX];
  97         x[index[i]][YY] =             cos(alfa)*x_old[YY] - sin(alfa)*x_old[ZZ];
  98         x[index[i]][ZZ] =             sin(alfa)*x_old[YY] + cos(alfa)*x_old[ZZ];
  99         copy_rvec(x[index[i]], x_old);
 100         /*calculate new x[index[i]] by rotation beta around the y-axis*/
 101         x[index[i]][XX] =   cos(beta)*x_old[XX]           + sin(beta)*x_old[ZZ];
 102         x[index[i]][YY] =                       x_old[YY];
 103         x[index[i]][ZZ] = -sin(beta)*x_old[XX]           + cos(beta)*x_old[ZZ];
 104         copy_rvec(x[index[i]], x_old);
 105         /*calculate new x[index[i]] by rotation gamma around the z-axis*/
 106         x[index[i]][XX] = x_old[XX]*cos(gamma) - x_old[YY]*sin(gamma);
 107         x[index[i]][YY] = x_old[XX]*sin(gamma) + x_old[YY]*cos(gamma);
 108         x[index[i]][ZZ] =                                             x_old[ZZ];
 109     }
 110 }
 111
 112 void rotate_conf(int natom, rvec *x, rvec *v, real alfa, real beta, real gamma)
 113 {
 114     if (x)
 115     {
 116         low_rotate_conf(natom, x, alfa, beta, gamma);
 117     }
 118     if (v)
 119     {
 120         low_rotate_conf(natom, v, alfa, beta, gamma);
 121     }
 122 }
 123
 124 void orient(int natom, rvec *x, rvec *v, rvec angle, matrix box)
 125 {
 126     real  longest, rij, rzi;
 127     int   i, j, m, max_i = 0, max_j = 0;
 128     rvec  origin;
 129     int   temp;
 130     real  alfa = 0, beta = 0, gamma = 0;
 131     t_pbc pbc;
 132
 133     set_pbc(&pbc, -1, box);
 134
 135     /*first i am going to look for the longest atom-atom distance*/
 136     longest = dist2(&pbc, x[0], x[1]);
 137     i       = 0;
 138     j       = 1;
 139     for (i = 0; (i < natom); i++)
 140     {
 141         for (j = 0; (j < natom); j++)
 142         {
 143             rij = dist2(&pbc, x[i], x[j]);
 144             if (rij > longest)
 145             {
 146                 max_i   = i;
 147                 max_j   = j;
 148                 longest = rij;
 149             }
 150         }
 151     }
 152     /* first check if x[max_i]<x[max_j] else swap*/
 153     if (x[max_i][2] > x[max_j][2])
 154     {
 155         temp  = max_i;
 156         max_i = max_j;
 157         max_j = temp;
 158     }
 159
 160     /*set the origin to x[i]*/
 161     for (m = 0; (m < DIM); m++)
 162     {
 163         origin[m] = x[max_i][m];
 164     }
 165     for (i = 0; (i < natom); i++)
 166     {
 167         for (m = 0; (m < DIM); m++)
 168         {
 169             x[i][m] -= origin[m];
 170         }
 171     }
 172
 173     /* calculate the rotation angles alfa(x_axis) and beta(y_axis)
 174      * the rotation angles must be calculated clockwise looking
 175      * along the rotation axis to the origin*
 176      * alfa (x-axis)
 177      */
 178     alfa = atan(x[max_j][ZZ]/x[max_j][YY])-M_PI_2;
 179     beta = M_PI_2-atan(x[max_j][ZZ]/x[max_j][XX]);
 180     rotate_conf(natom, x, v, alfa, beta, gamma);
 181
 182     /* now search the longest distance for rotation along the z_axis */
 183     longest = distance_to_z(x[0]);
 184     max_i   = 0;
 185     for (i = 1; (i < natom); i++)
 186     {
 187         rzi = distance_to_z(x[i]);
 188         if (rzi > longest)
 189         {
 190             longest = rzi;
 191             max_i   = i;
 192         }
 193     }
 194     gamma = atan(x[max_i][YY]/x[max_i][XX])-M_PI_2;
 195     rotate_conf(natom, x, v, 0, 0, gamma);
 196     angle[0] = alfa;
 197     angle[1] = beta;
 198     angle[2] = gamma;
 199 } /*orient()*/
 200
 201
 202 void genconf(t_atoms *atoms, rvec *x, rvec *v, real *r, matrix box, ivec n_box)
 203 {
 204     int     i, ix, iy, iz, m, j, imol, offset;
 205     rvec    delta;
 206     int     nmol;
 207
 208     nmol = n_box[XX]*n_box[YY]*n_box[ZZ];
 209
 210     /*print message*/
 211     fprintf(stderr, "Generating configuration\n");
 212     imol = 0;
 213     for (ix = 0; (ix < n_box[XX]); ix++)
 214     {
 215         delta[XX] = ix*box[XX][XX];
 216         for (iy = 0; (iy < n_box[YY]); iy++)
 217         {
 218             delta[YY] = iy*box[YY][YY];
 219             for (iz = 0; (iz < n_box[ZZ]); iz++)
 220             {
 221                 delta[ZZ] = iz*box[ZZ][ZZ];
 222                 offset    = imol*atoms->nr;
 223                 for (i = 0; (i < atoms->nr); i++)
 224                 {
 225                     for (m = 0; (m < DIM); m++)
 226                     {
 227                         x[offset+i][m] = delta[m]+x[i][m];
 228                     }
 229                     if (v)
 230                     {
 231                         for (m = 0; (m < DIM); m++)
 232                         {
 233                             v[offset+i][m] = v[i][m];
 234                         }
 235                     }
 236                     r[offset+i] = r[i];
 237                 }
 238                 imol++;
 239             }
 240         }
 241     }
 242     for (i = 1; (i < nmol); i++)
 243     {
 244         int offs    = i*atoms->nr;
 245         int offsres = i*atoms->nres;
 246         for (j = 0; (j < atoms->nr); j++)
 247         {
 248             atoms->atomname[offs+j]                    = atoms->atomname[j];
 249             atoms->atom[offs+j].resind                 = atoms->atom[j].resind + offsres;
 250             atoms->resinfo[atoms->atom[offs+j].resind] =
 251                 atoms->resinfo[atoms->atom[j].resind];
 252             atoms->resinfo[atoms->atom[offs+j].resind].nr += offsres;
 253         }
 254     }
 255     atoms->nr   *= nmol;
 256     atoms->nres *= nmol;
 257     for (i = 0; i < DIM; i++)
 258     {
 259         for (j = 0; j < DIM; j++)
 260         {
 261             box[j][i] *= n_box[j];
 262         }
 263     }
 264 } /*genconf()*/
 265
 266 /*gen_box() generates a box around a configuration*/
 267 void gen_box(int NTB, int natoms, rvec *x, matrix box, rvec box_space,
 268              gmx_bool bCenter)
 269 {
 270     int  i, m;
 271     rvec xmin, xmax;
 272     real max_box;
 273
 274     /*calculate minimum and maximum x[0..DIM-1]*/
 275     for (m = 0; (m < DIM); m++)
 276     {
 277         xmin[m] = xmax[m] = x[0][m];
 278     }
 279     for (i = 1; (i < natoms); i++)
 280     {
 281         for (m = 0; m < DIM; m++)
 282         {
 283             xmin[m] = min(xmin[m], x[i][m]);
 284             xmax[m] = max(xmax[m], x[i][m]);
 285         }
 286     }
 287
 288     /*calculate the new box sizes for cubic and octahedral ...*/
 289     for (m = 0; (m < DIM); m++)
 290     {
 291         box[m][m] = xmax[m]-xmin[m]+2*box_space[m];
 292     }
 293
 294     /*calculate the box size if NTB=1 (truncated octahedron)*/
 295     if (NTB == 1)
 296     {
 297         max_box = box[0][0];
 298         for (m = 0; (m < DIM); m++)
 299         {
 300             max_box = max(max_box, box[m][m]);
 301         }
 302         for (m = 0; (m < DIM); m++)
 303         {
 304             box[m][m] = max_box;
 305         }
 306     }
 307
 308     /*move the molecule to the center of the box*/
 309     if (bCenter)
 310     {
 311         for (i = 0; (i < natoms); i++)
 312         {
 313             for (m = 0; (m < DIM); m++)
 314             {
 315                 x[i][m] += 0.5*(box[m][m]-xmin[m]-xmax[m]);
 316             }
 317         }
 318     }
 319
 320
 321 #ifdef DEBUG
 322     /* print data to check this */
 323     print_stat(x, natoms, box);
 324 #endif
 325 } /*gen_box()*/