src/gromacs/pbcutil/pbcmethods.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmxpre.h"
  36
  37 #include "pbcmethods.h"
  38
  39 #include <algorithm>
  40 #include <memory>
  41
  42 #include "gromacs/pbcutil/pbc.h"
  43 #include "gromacs/topology/topology.h"
  44 #include "gromacs/utility/fatalerror.h"
  45 #include "gromacs/utility/smalloc.h"
  46
  47 void calc_pbc_cluster(int ecenter, int nrefat, t_topology *top, int ePBC,
  48                       rvec x[], const int index[], matrix box)
  49 {
  50     int       m, i, j, j0, j1, jj, ai, aj;
  51     int       imin, jmin;
  52     real      fac, min_dist2;
  53     rvec      dx, xtest, box_center;
  54     int       nmol, imol_center;
  55     int      *molind;
  56     gmx_bool *bMol, *bTmp;
  57     rvec     *m_com, *m_shift;
  58     t_pbc     pbc;
  59     int      *cluster;
  60     int      *added;
  61     int       ncluster, nadded;
  62     real      tmp_r2;
  63
  64     calc_box_center(ecenter, box, box_center);
  65
  66     /* Initiate the pbc structure */
  67     std::memset(&pbc, 0, sizeof(pbc));
  68     set_pbc(&pbc, ePBC, box);
  69
  70     /* Convert atom index to molecular */
  71     nmol   = top->mols.nr;
  72     molind = top->mols.index;
  73     snew(bMol, nmol);
  74     snew(m_com, nmol);
  75     snew(m_shift, nmol);
  76     snew(cluster, nmol);
  77     snew(added, nmol);
  78     snew(bTmp, top->atoms.nr);
  79
  80     for (i = 0; (i < nrefat); i++)
  81     {
  82         /* Mark all molecules in the index */
  83         ai       = index[i];
  84         bTmp[ai] = TRUE;
  85         /* Binary search assuming the molecules are sorted */
  86         j0 = 0;
  87         j1 = nmol-1;
  88         while (j0 < j1)
  89         {
  90             if (ai < molind[j0+1])
  91             {
  92                 j1 = j0;
  93             }
  94             else if (ai >= molind[j1])
  95             {
  96                 j0 = j1;
  97             }
  98             else
  99             {
 100                 jj = (j0+j1)/2;
 101                 if (ai < molind[jj+1])
 102                 {
 103                     j1 = jj;
 104                 }
 105                 else
 106                 {
 107                     j0 = jj;
 108                 }
 109             }
 110         }
 111         bMol[j0] = TRUE;
 112     }
 113     /* Double check whether all atoms in all molecules that are marked are part
 114      * of the cluster. Simultaneously compute the center of geometry.
 115      */
 116     min_dist2   = 10*gmx::square(trace(box));
 117     imol_center = -1;
 118     ncluster    = 0;
 119     for (i = 0; i < nmol; i++)
 120     {
 121         for (j = molind[i]; j < molind[i+1]; j++)
 122         {
 123             if (bMol[i] && !bTmp[j])
 124             {
 125                 gmx_fatal(FARGS, "Molecule %d marked for clustering but not atom %d in it - check your index!", i+1, j+1);
 126             }
 127             else if (!bMol[i] && bTmp[j])
 128             {
 129                 gmx_fatal(FARGS, "Atom %d marked for clustering but not molecule %d - this is an internal error...", j+1, i+1);
 130             }
 131             else if (bMol[i])
 132             {
 133                 /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
 134                 if (j > molind[i])
 135                 {
 136                     pbc_dx(&pbc, x[j], x[j-1], dx);
 137                     rvec_add(x[j-1], dx, x[j]);
 138                 }
 139                 /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
 140                 rvec_inc(m_com[i], x[j]);
 141             }
 142         }
 143         if (bMol[i])
 144         {
 145             /* Normalize center of geometry */
 146             fac = 1.0/(molind[i+1]-molind[i]);
 147             for (m = 0; (m < DIM); m++)
 148             {
 149                 m_com[i][m] *= fac;
 150             }
 151             /* Determine which molecule is closest to the center of the box */
 152             pbc_dx(&pbc, box_center, m_com[i], dx);
 153             tmp_r2 = iprod(dx, dx);
 154
 155             if (tmp_r2 < min_dist2)
 156             {
 157                 min_dist2   = tmp_r2;
 158                 imol_center = i;
 159             }
 160             cluster[ncluster++] = i;
 161         }
 162     }
 163     sfree(bTmp);
 164
 165     if (ncluster <= 0)
 166     {
 167         fprintf(stderr, "No molecules selected in the cluster\n");
 168         return;
 169     }
 170     else if (imol_center == -1)
 171     {
 172         fprintf(stderr, "No central molecules could be found\n");
 173         return;
 174     }
 175
 176     nadded            = 0;
 177     added[nadded++]   = imol_center;
 178     bMol[imol_center] = FALSE;
 179
 180     while (nadded < ncluster)
 181     {
 182         /* Find min distance between cluster molecules and those remaining to be added */
 183         min_dist2   = 10*gmx::square(trace(box));
 184         imin        = -1;
 185         jmin        = -1;
 186         /* Loop over added mols */
 187         for (i = 0; i < nadded; i++)
 188         {
 189             ai = added[i];
 190             /* Loop over all mols */
 191             for (j = 0; j < ncluster; j++)
 192             {
 193                 aj = cluster[j];
 194                 /* check those remaining to be added */
 195                 if (bMol[aj])
 196                 {
 197                     pbc_dx(&pbc, m_com[aj], m_com[ai], dx);
 198                     tmp_r2 = iprod(dx, dx);
 199                     if (tmp_r2 < min_dist2)
 200                     {
 201                         min_dist2   = tmp_r2;
 202                         imin        = ai;
 203                         jmin        = aj;
 204                     }
 205                 }
 206             }
 207         }
 208
 209         /* Add the best molecule */
 210         added[nadded++]   = jmin;
 211         bMol[jmin]        = FALSE;
 212         /* Calculate the shift from the ai molecule */
 213         pbc_dx(&pbc, m_com[jmin], m_com[imin], dx);
 214         rvec_add(m_com[imin], dx, xtest);
 215         rvec_sub(xtest, m_com[jmin], m_shift[jmin]);
 216         rvec_inc(m_com[jmin], m_shift[jmin]);
 217
 218         for (j = molind[jmin]; j < molind[jmin+1]; j++)
 219         {
 220             rvec_inc(x[j], m_shift[jmin]);
 221         }
 222         fprintf(stdout, "\rClustering iteration %d of %d...", nadded, ncluster);
 223         fflush(stdout);
 224     }
 225
 226     sfree(added);
 227     sfree(cluster);
 228     sfree(bMol);
 229     sfree(m_com);
 230     sfree(m_shift);
 231
 232     fprintf(stdout, "\n");
 233 }
 234
 235 void put_molecule_com_in_box(int unitcell_enum, int ecenter,
 236                              t_block *mols,
 237                              int natoms, t_atom atom[],
 238                              int ePBC, matrix box, rvec x[])
 239 {
 240     int     i, j;
 241     int     d;
 242     rvec    com, shift, box_center;
 243     real    m;
 244     double  mtot;
 245     t_pbc   pbc;
 246
 247     calc_box_center(ecenter, box, box_center);
 248     set_pbc(&pbc, ePBC, box);
 249     if (mols->nr <= 0)
 250     {
 251         gmx_fatal(FARGS, "There are no molecule descriptions. I need a .tpr file for this pbc option.");
 252     }
 253     for (i = 0; (i < mols->nr); i++)
 254     {
 255         /* calc COM */
 256         clear_rvec(com);
 257         mtot = 0;
 258         for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
 259         {
 260             m = atom[j].m;
 261             for (d = 0; d < DIM; d++)
 262             {
 263                 com[d] += m*x[j][d];
 264             }
 265             mtot += m;
 266         }
 267         /* calculate final COM */
 268         svmul(1.0/mtot, com, com);
 269
 270         /* check if COM is outside box */
 271         gmx::RVec newCom;
 272         copy_rvec(com, newCom);
 273         auto      newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
 274         switch (unitcell_enum)
 275         {
 276             case euRect:
 277                 put_atoms_in_box(ePBC, box, newComArrayRef);
 278                 break;
 279             case euTric:
 280                 put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef);
 281                 break;
 282             case euCompact:
 283                 put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
 284                 break;
 285         }
 286         rvec_sub(newCom, com, shift);
 287         if (norm2(shift) > 0)
 288         {
 289             if (debug)
 290             {
 291                 fprintf(debug, "\nShifting position of molecule %d "
 292                         "by %8.3f  %8.3f  %8.3f\n", i+1,
 293                         shift[XX], shift[YY], shift[ZZ]);
 294             }
 295             for (j = mols->index[i]; (j < mols->index[i+1] && j < natoms); j++)
 296             {
 297                 rvec_inc(x[j], shift);
 298             }
 299         }
 300     }
 301 }
 302
 303 void put_residue_com_in_box(int unitcell_enum, int ecenter,
 304                             int natoms, t_atom atom[],
 305                             int ePBC, matrix box, rvec x[])
 306 {
 307     int              i, j, res_start, res_end;
 308     int              d, presnr;
 309     real             m;
 310     double           mtot;
 311     rvec             box_center, com, shift;
 312     static const int NOTSET = -12347;
 313     calc_box_center(ecenter, box, box_center);
 314
 315     presnr    = NOTSET;
 316     res_start = 0;
 317     clear_rvec(com);
 318     mtot = 0;
 319     for (i = 0; i < natoms+1; i++)
 320     {
 321         if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET))
 322         {
 323             /* calculate final COM */
 324             res_end = i;
 325             svmul(1.0/mtot, com, com);
 326
 327             /* check if COM is outside box */
 328             gmx::RVec newCom;
 329             copy_rvec(com, newCom);
 330             auto      newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
 331             switch (unitcell_enum)
 332             {
 333                 case euRect:
 334                     put_atoms_in_box(ePBC, box, newComArrayRef);
 335                     break;
 336                 case euTric:
 337                     put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef);
 338                     break;
 339                 case euCompact:
 340                     put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
 341                     break;
 342             }
 343             rvec_sub(newCom, com, shift);
 344             if (norm2(shift) != 0.0F)
 345             {
 346                 if (debug)
 347                 {
 348                     fprintf(debug, "\nShifting position of residue %d (atoms %d-%d) "
 349                             "by %g,%g,%g\n", atom[res_start].resind+1,
 350                             res_start+1, res_end+1, shift[XX], shift[YY], shift[ZZ]);
 351                 }
 352                 for (j = res_start; j < res_end; j++)
 353                 {
 354                     rvec_inc(x[j], shift);
 355                 }
 356             }
 357             clear_rvec(com);
 358             mtot = 0;
 359
 360             /* remember start of new residue */
 361             res_start = i;
 362         }
 363         if (i < natoms)
 364         {
 365             /* calc COM */
 366             m = atom[i].m;
 367             for (d = 0; d < DIM; d++)
 368             {
 369                 com[d] += m*x[i][d];
 370             }
 371             mtot += m;
 372
 373             presnr = atom[i].resind;
 374         }
 375     }
 376 }
 377
 378 void center_x(int ecenter, rvec x[], matrix box, int n, int nc, const int ci[])
 379 {
 380     int  i, m, ai;
 381     rvec cmin, cmax, box_center, dx;
 382
 383     if (nc > 0)
 384     {
 385         copy_rvec(x[ci[0]], cmin);
 386         copy_rvec(x[ci[0]], cmax);
 387         for (i = 0; i < nc; i++)
 388         {
 389             ai = ci[i];
 390             for (m = 0; m < DIM; m++)
 391             {
 392                 if (x[ai][m] < cmin[m])
 393                 {
 394                     cmin[m] = x[ai][m];
 395                 }
 396                 else if (x[ai][m] > cmax[m])
 397                 {
 398                     cmax[m] = x[ai][m];
 399                 }
 400             }
 401         }
 402         calc_box_center(ecenter, box, box_center);
 403         for (m = 0; m < DIM; m++)
 404         {
 405             dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
 406         }
 407
 408         for (i = 0; i < n; i++)
 409         {
 410             rvec_inc(x[i], dx);
 411         }
 412     }
 413 }