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, rvec x[], const int index[], matrix box)
  48 {
  49     int       m, i, j, j0, j1, jj, ai, aj;
  50     int       imin, jmin;
  51     real      fac, min_dist2;
  52     rvec      dx, xtest, box_center;
  53     int       nmol, imol_center;
  54     int*      molind;
  55     gmx_bool *bMol, *bTmp;
  56     rvec *    m_com, *m_shift;
  57     t_pbc     pbc;
  58     int*      cluster;
  59     int*      added;
  60     int       ncluster, nadded;
  61     real      tmp_r2;
  62
  63     calc_box_center(ecenter, box, box_center);
  64
  65     /* Initiate the pbc structure */
  66     std::memset(&pbc, 0, sizeof(pbc));
  67     set_pbc(&pbc, ePBC, box);
  68
  69     /* Convert atom index to molecular */
  70     nmol   = top->mols.nr;
  71     molind = top->mols.index;
  72     snew(bMol, nmol);
  73     snew(m_com, nmol);
  74     snew(m_shift, nmol);
  75     snew(cluster, nmol);
  76     snew(added, nmol);
  77     snew(bTmp, top->atoms.nr);
  78
  79     for (i = 0; (i < nrefat); i++)
  80     {
  81         /* Mark all molecules in the index */
  82         ai       = index[i];
  83         bTmp[ai] = TRUE;
  84         /* Binary search assuming the molecules are sorted */
  85         j0 = 0;
  86         j1 = nmol - 1;
  87         while (j0 < j1)
  88         {
  89             if (ai < molind[j0 + 1])
  90             {
  91                 j1 = j0;
  92             }
  93             else if (ai >= molind[j1])
  94             {
  95                 j0 = j1;
  96             }
  97             else
  98             {
  99                 jj = (j0 + j1) / 2;
 100                 if (ai < molind[jj + 1])
 101                 {
 102                     j1 = jj;
 103                 }
 104                 else
 105                 {
 106                     j0 = jj;
 107                 }
 108             }
 109         }
 110         bMol[j0] = TRUE;
 111     }
 112     /* Double check whether all atoms in all molecules that are marked are part
 113      * of the cluster. Simultaneously compute the center of geometry.
 114      */
 115     min_dist2   = 10 * gmx::square(trace(box));
 116     imol_center = -1;
 117     ncluster    = 0;
 118     for (i = 0; i < nmol; i++)
 119     {
 120         for (j = molind[i]; j < molind[i + 1]; j++)
 121         {
 122             if (bMol[i] && !bTmp[j])
 123             {
 124                 gmx_fatal(FARGS,
 125                           "Molecule %d marked for clustering but not atom %d in it - check your "
 126                           "index!",
 127                           i + 1, j + 1);
 128             }
 129             else if (!bMol[i] && bTmp[j])
 130             {
 131                 gmx_fatal(FARGS,
 132                           "Atom %d marked for clustering but not molecule %d - this is an internal "
 133                           "error...",
 134                           j + 1, i + 1);
 135             }
 136             else if (bMol[i])
 137             {
 138                 /* Make molecule whole, move 2nd and higher atom to same periodicity as 1st atom in molecule */
 139                 if (j > molind[i])
 140                 {
 141                     pbc_dx(&pbc, x[j], x[j - 1], dx);
 142                     rvec_add(x[j - 1], dx, x[j]);
 143                 }
 144                 /* Compute center of geometry of molecule - m_com[i] was zeroed when we did snew() on it! */
 145                 rvec_inc(m_com[i], x[j]);
 146             }
 147         }
 148         if (bMol[i])
 149         {
 150             /* Normalize center of geometry */
 151             fac = 1.0 / (molind[i + 1] - molind[i]);
 152             for (m = 0; (m < DIM); m++)
 153             {
 154                 m_com[i][m] *= fac;
 155             }
 156             /* Determine which molecule is closest to the center of the box */
 157             pbc_dx(&pbc, box_center, m_com[i], dx);
 158             tmp_r2 = iprod(dx, dx);
 159
 160             if (tmp_r2 < min_dist2)
 161             {
 162                 min_dist2   = tmp_r2;
 163                 imol_center = i;
 164             }
 165             cluster[ncluster++] = i;
 166         }
 167     }
 168     sfree(bTmp);
 169
 170     if (ncluster <= 0)
 171     {
 172         fprintf(stderr, "No molecules selected in the cluster\n");
 173         return;
 174     }
 175     else if (imol_center == -1)
 176     {
 177         fprintf(stderr, "No central molecules could be found\n");
 178         return;
 179     }
 180
 181     nadded            = 0;
 182     added[nadded++]   = imol_center;
 183     bMol[imol_center] = FALSE;
 184
 185     while (nadded < ncluster)
 186     {
 187         /* Find min distance between cluster molecules and those remaining to be added */
 188         min_dist2 = 10 * gmx::square(trace(box));
 189         imin      = -1;
 190         jmin      = -1;
 191         /* Loop over added mols */
 192         for (i = 0; i < nadded; i++)
 193         {
 194             ai = added[i];
 195             /* Loop over all mols */
 196             for (j = 0; j < ncluster; j++)
 197             {
 198                 aj = cluster[j];
 199                 /* check those remaining to be added */
 200                 if (bMol[aj])
 201                 {
 202                     pbc_dx(&pbc, m_com[aj], m_com[ai], dx);
 203                     tmp_r2 = iprod(dx, dx);
 204                     if (tmp_r2 < min_dist2)
 205                     {
 206                         min_dist2 = tmp_r2;
 207                         imin      = ai;
 208                         jmin      = aj;
 209                     }
 210                 }
 211             }
 212         }
 213
 214         /* Add the best molecule */
 215         added[nadded++] = jmin;
 216         bMol[jmin]      = FALSE;
 217         /* Calculate the shift from the ai molecule */
 218         pbc_dx(&pbc, m_com[jmin], m_com[imin], dx);
 219         rvec_add(m_com[imin], dx, xtest);
 220         rvec_sub(xtest, m_com[jmin], m_shift[jmin]);
 221         rvec_inc(m_com[jmin], m_shift[jmin]);
 222
 223         for (j = molind[jmin]; j < molind[jmin + 1]; j++)
 224         {
 225             rvec_inc(x[j], m_shift[jmin]);
 226         }
 227         fprintf(stdout, "\rClustering iteration %d of %d...", nadded, ncluster);
 228         fflush(stdout);
 229     }
 230
 231     sfree(added);
 232     sfree(cluster);
 233     sfree(bMol);
 234     sfree(m_com);
 235     sfree(m_shift);
 236
 237     fprintf(stdout, "\n");
 238 }
 239
 240 void put_molecule_com_in_box(int      unitcell_enum,
 241                              int      ecenter,
 242                              t_block* mols,
 243                              int      natoms,
 244                              t_atom   atom[],
 245                              int      ePBC,
 246                              matrix   box,
 247                              rvec     x[])
 248 {
 249     int    i, j;
 250     int    d;
 251     rvec   com, shift, box_center;
 252     real   m;
 253     double mtot;
 254     t_pbc  pbc;
 255
 256     calc_box_center(ecenter, box, box_center);
 257     set_pbc(&pbc, ePBC, box);
 258     if (mols->nr <= 0)
 259     {
 260         gmx_fatal(FARGS,
 261                   "There are no molecule descriptions. I need a .tpr file for this pbc option.");
 262     }
 263     for (i = 0; (i < mols->nr); i++)
 264     {
 265         /* calc COM */
 266         clear_rvec(com);
 267         mtot = 0;
 268         for (j = mols->index[i]; (j < mols->index[i + 1] && j < natoms); j++)
 269         {
 270             m = atom[j].m;
 271             for (d = 0; d < DIM; d++)
 272             {
 273                 com[d] += m * x[j][d];
 274             }
 275             mtot += m;
 276         }
 277         /* calculate final COM */
 278         svmul(1.0 / mtot, com, com);
 279
 280         /* check if COM is outside box */
 281         gmx::RVec newCom;
 282         copy_rvec(com, newCom);
 283         auto newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
 284         switch (unitcell_enum)
 285         {
 286             case euRect: put_atoms_in_box(ePBC, box, newComArrayRef); break;
 287             case euTric: put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef); break;
 288             case euCompact:
 289                 put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
 290                 break;
 291         }
 292         rvec_sub(newCom, com, shift);
 293         if (norm2(shift) > 0)
 294         {
 295             if (debug)
 296             {
 297                 fprintf(debug,
 298                         "\nShifting position of molecule %d "
 299                         "by %8.3f  %8.3f  %8.3f\n",
 300                         i + 1, shift[XX], shift[YY], shift[ZZ]);
 301             }
 302             for (j = mols->index[i]; (j < mols->index[i + 1] && j < natoms); j++)
 303             {
 304                 rvec_inc(x[j], shift);
 305             }
 306         }
 307     }
 308 }
 309
 310 void put_residue_com_in_box(int unitcell_enum, int ecenter, int natoms, t_atom atom[], int ePBC, matrix box, rvec x[])
 311 {
 312     int              i, j, res_start, res_end;
 313     int              d, presnr;
 314     real             m;
 315     double           mtot;
 316     rvec             box_center, com, shift;
 317     static const int NOTSET = -12347;
 318     calc_box_center(ecenter, box, box_center);
 319
 320     presnr    = NOTSET;
 321     res_start = 0;
 322     clear_rvec(com);
 323     mtot = 0;
 324     for (i = 0; i < natoms + 1; i++)
 325     {
 326         if (i == natoms || (presnr != atom[i].resind && presnr != NOTSET))
 327         {
 328             /* calculate final COM */
 329             res_end = i;
 330             svmul(1.0 / mtot, com, com);
 331
 332             /* check if COM is outside box */
 333             gmx::RVec newCom;
 334             copy_rvec(com, newCom);
 335             auto newComArrayRef = gmx::arrayRefFromArray(&newCom, 1);
 336             switch (unitcell_enum)
 337             {
 338                 case euRect: put_atoms_in_box(ePBC, box, newComArrayRef); break;
 339                 case euTric: put_atoms_in_triclinic_unitcell(ecenter, box, newComArrayRef); break;
 340                 case euCompact:
 341                     put_atoms_in_compact_unitcell(ePBC, ecenter, box, newComArrayRef);
 342                     break;
 343             }
 344             rvec_sub(newCom, com, shift);
 345             if (norm2(shift) != 0.0F)
 346             {
 347                 if (debug)
 348                 {
 349                     fprintf(debug,
 350                             "\nShifting position of residue %d (atoms %d-%d) "
 351                             "by %g,%g,%g\n",
 352                             atom[res_start].resind + 1, res_start + 1, res_end + 1, shift[XX],
 353                             shift[YY], shift[ZZ]);
 354                 }
 355                 for (j = res_start; j < res_end; j++)
 356                 {
 357                     rvec_inc(x[j], shift);
 358                 }
 359             }
 360             clear_rvec(com);
 361             mtot = 0;
 362
 363             /* remember start of new residue */
 364             res_start = i;
 365         }
 366         if (i < natoms)
 367         {
 368             /* calc COM */
 369             m = atom[i].m;
 370             for (d = 0; d < DIM; d++)
 371             {
 372                 com[d] += m * x[i][d];
 373             }
 374             mtot += m;
 375
 376             presnr = atom[i].resind;
 377         }
 378     }
 379 }
 380
 381 void center_x(int ecenter, rvec x[], matrix box, int n, int nc, const int ci[])
 382 {
 383     int  i, m, ai;
 384     rvec cmin, cmax, box_center, dx;
 385
 386     if (nc > 0)
 387     {
 388         copy_rvec(x[ci[0]], cmin);
 389         copy_rvec(x[ci[0]], cmax);
 390         for (i = 0; i < nc; i++)
 391         {
 392             ai = ci[i];
 393             for (m = 0; m < DIM; m++)
 394             {
 395                 if (x[ai][m] < cmin[m])
 396                 {
 397                     cmin[m] = x[ai][m];
 398                 }
 399                 else if (x[ai][m] > cmax[m])
 400                 {
 401                     cmax[m] = x[ai][m];
 402                 }
 403             }
 404         }
 405         calc_box_center(ecenter, box, box_center);
 406         for (m = 0; m < DIM; m++)
 407         {
 408             dx[m] = box_center[m] - (cmin[m] + cmax[m]) * 0.5;
 409         }
 410
 411         for (i = 0; i < n; i++)
 412         {
 413             rvec_inc(x[i], dx);
 414         }
 415     }
 416 }