src/gromacs/pulling/pull_rotation.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2008, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "pull_rotation.h"
  41
  42 #include "config.h"
  43
  44 #include <cstdio>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include <algorithm>
  49 #include <memory>
  50
  51 #include "gromacs/commandline/filenm.h"
  52 #include "gromacs/domdec/dlbtiming.h"
  53 #include "gromacs/domdec/domdec_struct.h"
  54 #include "gromacs/domdec/ga2la.h"
  55 #include "gromacs/domdec/localatomset.h"
  56 #include "gromacs/domdec/localatomsetmanager.h"
  57 #include "gromacs/fileio/gmxfio.h"
  58 #include "gromacs/fileio/xvgr.h"
  59 #include "gromacs/gmxlib/network.h"
  60 #include "gromacs/linearalgebra/nrjac.h"
  61 #include "gromacs/math/functions.h"
  62 #include "gromacs/math/units.h"
  63 #include "gromacs/math/utilities.h"
  64 #include "gromacs/math/vec.h"
  65 #include "gromacs/mdlib/groupcoord.h"
  66 #include "gromacs/mdlib/stat.h"
  67 #include "gromacs/mdrunutility/handlerestart.h"
  68 #include "gromacs/mdtypes/commrec.h"
  69 #include "gromacs/mdtypes/inputrec.h"
  70 #include "gromacs/mdtypes/md_enums.h"
  71 #include "gromacs/mdtypes/mdrunoptions.h"
  72 #include "gromacs/mdtypes/state.h"
  73 #include "gromacs/pbcutil/pbc.h"
  74 #include "gromacs/timing/cyclecounter.h"
  75 #include "gromacs/timing/wallcycle.h"
  76 #include "gromacs/topology/mtop_lookup.h"
  77 #include "gromacs/topology/mtop_util.h"
  78 #include "gromacs/utility/basedefinitions.h"
  79 #include "gromacs/utility/fatalerror.h"
  80 #include "gromacs/utility/pleasecite.h"
  81 #include "gromacs/utility/smalloc.h"
  82
  83 static const std::string RotStr = { "Enforced rotation:" };
  84
  85 /* Set the minimum weight for the determination of the slab centers */
  86 #define WEIGHT_MIN (10 * GMX_FLOAT_MIN)
  87
  88 //! Helper structure for sorting positions along rotation vector
  89 struct sort_along_vec_t
  90 {
  91     //! Projection of xc on the rotation vector
  92     real xcproj;
  93     //! Index of xc
  94     int ind;
  95     //! Mass
  96     real m;
  97     //! Position
  98     rvec x;
  99     //! Reference position
 100     rvec x_ref;
 101 };
 102
 103
 104 //! Enforced rotation / flexible: determine the angle of each slab
 105 struct gmx_slabdata
 106 {
 107     //! Number of atoms belonging to this slab
 108     int nat;
 109     /*! \brief The positions belonging to this slab.
 110      *
 111      * In general, this should be all positions of the whole
 112      * rotation group, but we leave those away that have a small
 113      * enough weight. */
 114     rvec* x;
 115     //! Same for reference
 116     rvec* ref;
 117     //! The weight for each atom
 118     real* weight;
 119 };
 120
 121
 122 //! Helper structure for potential fitting
 123 struct gmx_potfit
 124 {
 125     /*! \brief Set of angles for which the potential is calculated.
 126      *
 127      * The optimum fit is determined as the angle for with the
 128      * potential is minimal. */
 129     real* degangle;
 130     //! Potential for the different angles
 131     real* V;
 132     //! Rotation matrix corresponding to the angles
 133     matrix* rotmat;
 134 };
 135
 136
 137 //! Enforced rotation data for a single rotation group
 138 struct gmx_enfrotgrp
 139 {
 140     //! Input parameters for this group
 141     const t_rotgrp* rotg = nullptr;
 142     //! Index of this group within the set of groups
 143     int groupIndex;
 144     //! Rotation angle in degrees
 145     real degangle;
 146     //! Rotation matrix
 147     matrix rotmat;
 148     //! The atoms subject to enforced rotation
 149     std::unique_ptr<gmx::LocalAtomSet> atomSet;
 150
 151     //! The normalized rotation vector
 152     rvec vec;
 153     //! Rotation potential for this rotation group
 154     real V;
 155     //! Array to store the forces on the local atoms resulting from enforced rotation potential
 156     rvec* f_rot_loc;
 157
 158     /* Collective coordinates for the whole rotation group */
 159     //! Length of each x_rotref vector after x_rotref has been put into origin
 160     real* xc_ref_length;
 161     //! Center of the rotation group positions, may be mass weighted
 162     rvec xc_center;
 163     //! Center of the rotation group reference positions
 164     rvec xc_ref_center;
 165     //! Current (collective) positions
 166     rvec* xc;
 167     //! Current (collective) shifts
 168     ivec* xc_shifts;
 169     //! Extra shifts since last DD step
 170     ivec* xc_eshifts;
 171     //! Old (collective) positions
 172     rvec* xc_old;
 173     //! Normalized form of the current positions
 174     rvec* xc_norm;
 175     //! Reference positions (sorted in the same order as xc when sorted)
 176     rvec* xc_ref_sorted;
 177     //! Where is a position found after sorting?
 178     int* xc_sortind;
 179     //! Collective masses
 180     real* mc;
 181     //! Collective masses sorted
 182     real* mc_sorted;
 183     //! one over the total mass of the rotation group
 184     real invmass;
 185
 186     //! Torque in the direction of rotation vector
 187     real torque_v;
 188     //! Actual angle of the whole rotation group
 189     real angle_v;
 190     /* Fixed rotation only */
 191     //! Weights for angle determination
 192     real weight_v;
 193     //! Local reference coords, correctly rotated
 194     rvec* xr_loc;
 195     //! Local current coords, correct PBC image
 196     rvec* x_loc_pbc;
 197     //! Masses of the current local atoms
 198     real* m_loc;
 199
 200     /* Flexible rotation only */
 201     //! For this many slabs memory is allocated
 202     int nslabs_alloc;
 203     //! Lowermost slab for that the calculation needs to be performed at a given time step
 204     int slab_first;
 205     //! Uppermost slab ...
 206     int slab_last;
 207     //! First slab for which ref. center is stored
 208     int slab_first_ref;
 209     //! Last ...
 210     int slab_last_ref;
 211     //! Slab buffer region around reference slabs
 212     int slab_buffer;
 213     //! First relevant atom for a slab
 214     int* firstatom;
 215     //! Last relevant atom for a slab
 216     int* lastatom;
 217     //! Gaussian-weighted slab center
 218     rvec* slab_center;
 219     //! Gaussian-weighted slab center for the reference positions
 220     rvec* slab_center_ref;
 221     //! Sum of gaussian weights in a slab
 222     real* slab_weights;
 223     //! Torque T = r x f for each slab. torque_v = m.v = angular momentum in the direction of v
 224     real* slab_torque_v;
 225     //! min_gaussian from t_rotgrp is the minimum value the gaussian must have so that the force is actually evaluated. max_beta is just another way to put it
 226     real max_beta;
 227     //! Precalculated gaussians for a single atom
 228     real* gn_atom;
 229     //! Tells to which slab each precalculated gaussian belongs
 230     int* gn_slabind;
 231     //! Inner sum of the flexible2 potential per slab; this is precalculated for optimization reasons
 232     rvec* slab_innersumvec;
 233     //! Holds atom positions and gaussian weights of atoms belonging to a slab
 234     gmx_slabdata* slab_data;
 235
 236     /* For potential fits with varying angle: */
 237     //! Used for fit type 'potential'
 238     gmx_potfit* PotAngleFit;
 239 };
 240
 241
 242 //! Enforced rotation data for all groups
 243 struct gmx_enfrot
 244 {
 245     //! Input parameters.
 246     const t_rot* rot = nullptr;
 247     //! Output period for main rotation outfile
 248     int nstrout;
 249     //! Output period for per-slab data
 250     int nstsout;
 251     //! Output file for rotation data
 252     FILE* out_rot = nullptr;
 253     //! Output file for torque data
 254     FILE* out_torque = nullptr;
 255     //! Output file for slab angles for flexible type
 256     FILE* out_angles = nullptr;
 257     //! Output file for slab centers
 258     FILE* out_slabs = nullptr;
 259     //! Allocation size of buf
 260     int bufsize = 0;
 261     //! Coordinate buffer variable for sorting
 262     rvec* xbuf = nullptr;
 263     //! Masses buffer variable for sorting
 264     real* mbuf = nullptr;
 265     //! Buffer variable needed for position sorting
 266     sort_along_vec_t* data = nullptr;
 267     //! MPI buffer
 268     real* mpi_inbuf = nullptr;
 269     //! MPI buffer
 270     real* mpi_outbuf = nullptr;
 271     //! Allocation size of in & outbuf
 272     int mpi_bufsize = 0;
 273     //! If true, append output files
 274     gmx_bool restartWithAppending = false;
 275     //! Used to skip first output when appending to avoid duplicate entries in rotation outfiles
 276     gmx_bool bOut = false;
 277     //! Stores working data per group
 278     std::vector<gmx_enfrotgrp> enfrotgrp;
 279     ~gmx_enfrot();
 280 };
 281
 282 gmx_enfrot::~gmx_enfrot()
 283 {
 284     if (out_rot)
 285     {
 286         gmx_fio_fclose(out_rot);
 287     }
 288     if (out_slabs)
 289     {
 290         gmx_fio_fclose(out_slabs);
 291     }
 292     if (out_angles)
 293     {
 294         gmx_fio_fclose(out_angles);
 295     }
 296     if (out_torque)
 297     {
 298         gmx_fio_fclose(out_torque);
 299     }
 300 }
 301
 302 namespace gmx
 303 {
 304
 305 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
 306
 307 class EnforcedRotation::Impl
 308 {
 309 public:
 310     gmx_enfrot enforcedRotation_;
 311 };
 312
 313 EnforcedRotation::EnforcedRotation() : impl_(new Impl) {}
 314
 315 EnforcedRotation::~EnforcedRotation() = default;
 316
 317 gmx_enfrot* EnforcedRotation::getLegacyEnfrot()
 318 {
 319     return &impl_->enforcedRotation_;
 320 }
 321
 322 } // namespace gmx
 323
 324 /* Activate output of forces for correctness checks */
 325 /* #define PRINT_FORCES */
 326 #ifdef PRINT_FORCES
 327 #    define PRINT_FORCE_J                       \
 328         fprintf(stderr,                         \
 329                 "f%d = %15.8f %15.8f %15.8f\n", \
 330                 erg->xc_ref_ind[j],             \
 331                 erg->f_rot_loc[j][XX],          \
 332                 erg->f_rot_loc[j][YY],          \
 333                 erg->f_rot_loc[j][ZZ]);
 334 #    define PRINT_POT_TAU                   \
 335         if (MASTER(cr))                     \
 336         {                                   \
 337             fprintf(stderr,                 \
 338                     "potential = %15.8f\n"  \
 339                     "torque    = %15.8f\n", \
 340                     erg->V,                 \
 341                     erg->torque_v);         \
 342         }
 343 #else
 344 #    define PRINT_FORCE_J
 345 #    define PRINT_POT_TAU
 346 #endif
 347
 348 /* Shortcuts for often used queries */
 349 #define ISFLEX(rg)                                                                                         \
 350     (((rg)->eType == EnforcedRotationGroupType::Flex) || ((rg)->eType == EnforcedRotationGroupType::Flext) \
 351      || ((rg)->eType == EnforcedRotationGroupType::Flex2)                                                  \
 352      || ((rg)->eType == EnforcedRotationGroupType::Flex2t))
 353 #define ISCOLL(rg)                                                                                         \
 354     (((rg)->eType == EnforcedRotationGroupType::Flex) || ((rg)->eType == EnforcedRotationGroupType::Flext) \
 355      || ((rg)->eType == EnforcedRotationGroupType::Flex2)                                                  \
 356      || ((rg)->eType == EnforcedRotationGroupType::Flex2t)                                                 \
 357      || ((rg)->eType == EnforcedRotationGroupType::Rmpf)                                                   \
 358      || ((rg)->eType == EnforcedRotationGroupType::Rm2pf))
 359
 360
 361 /* Does any of the rotation groups use slab decomposition? */
 362 static gmx_bool HaveFlexibleGroups(const t_rot* rot)
 363 {
 364     for (int g = 0; g < rot->ngrp; g++)
 365     {
 366         if (ISFLEX(&rot->grp[g]))
 367         {
 368             return TRUE;
 369         }
 370     }
 371
 372     return FALSE;
 373 }
 374
 375
 376 /* Is for any group the fit angle determined by finding the minimum of the
 377  * rotation potential? */
 378 static gmx_bool HavePotFitGroups(const t_rot* rot)
 379 {
 380     for (int g = 0; g < rot->ngrp; g++)
 381     {
 382         if (RotationGroupFitting::Pot == rot->grp[g].eFittype)
 383         {
 384             return TRUE;
 385         }
 386     }
 387
 388     return FALSE;
 389 }
 390
 391
 392 static double** allocate_square_matrix(int dim)
 393 {
 394     int      i;
 395     double** mat = nullptr;
 396
 397
 398     snew(mat, dim);
 399     for (i = 0; i < dim; i++)
 400     {
 401         snew(mat[i], dim);
 402     }
 403
 404     return mat;
 405 }
 406
 407
 408 static void free_square_matrix(double** mat, int dim)
 409 {
 410     int i;
 411
 412
 413     for (i = 0; i < dim; i++)
 414     {
 415         sfree(mat[i]);
 416     }
 417     sfree(mat);
 418 }
 419
 420
 421 /* Return the angle for which the potential is minimal */
 422 static real get_fitangle(const gmx_enfrotgrp* erg)
 423 {
 424     int         i;
 425     real        fitangle = -999.9;
 426     real        pot_min  = GMX_FLOAT_MAX;
 427     gmx_potfit* fit;
 428
 429
 430     fit = erg->PotAngleFit;
 431
 432     for (i = 0; i < erg->rotg->PotAngle_nstep; i++)
 433     {
 434         if (fit->V[i] < pot_min)
 435         {
 436             pot_min  = fit->V[i];
 437             fitangle = fit->degangle[i];
 438         }
 439     }
 440
 441     return fitangle;
 442 }
 443
 444
 445 /* Reduce potential angle fit data for this group at this time step? */
 446 static inline gmx_bool bPotAngle(const gmx_enfrot* er, const t_rotgrp* rotg, int64_t step)
 447 {
 448     return ((RotationGroupFitting::Pot == rotg->eFittype)
 449             && (do_per_step(step, er->nstsout) || do_per_step(step, er->nstrout)));
 450 }
 451
 452 /* Reduce slab torqe data for this group at this time step? */
 453 static inline gmx_bool bSlabTau(const gmx_enfrot* er, const t_rotgrp* rotg, int64_t step)
 454 {
 455     return ((ISFLEX(rotg)) && do_per_step(step, er->nstsout));
 456 }
 457
 458 /* Output rotation energy, torques, etc. for each rotation group */
 459 static void reduce_output(const t_commrec* cr, gmx_enfrot* er, real t, int64_t step)
 460 {
 461     int      i, islab, nslabs = 0;
 462     int      count; /* MPI element counter                               */
 463     real     fitangle;
 464     gmx_bool bFlex;
 465
 466     /* Fill the MPI buffer with stuff to reduce. If items are added for reduction
 467      * here, the MPI buffer size has to be enlarged also in calc_mpi_bufsize() */
 468     if (PAR(cr))
 469     {
 470         count = 0;
 471         for (auto& ergRef : er->enfrotgrp)
 472         {
 473             gmx_enfrotgrp*  erg    = &ergRef;
 474             const t_rotgrp* rotg   = erg->rotg;
 475             nslabs                 = erg->slab_last - erg->slab_first + 1;
 476             er->mpi_inbuf[count++] = erg->V;
 477             er->mpi_inbuf[count++] = erg->torque_v;
 478             er->mpi_inbuf[count++] = erg->angle_v;
 479             er->mpi_inbuf[count++] =
 480                     erg->weight_v; /* weights are not needed for flex types, but this is just a single value */
 481
 482             if (bPotAngle(er, rotg, step))
 483             {
 484                 for (i = 0; i < rotg->PotAngle_nstep; i++)
 485                 {
 486                     er->mpi_inbuf[count++] = erg->PotAngleFit->V[i];
 487                 }
 488             }
 489             if (bSlabTau(er, rotg, step))
 490             {
 491                 for (i = 0; i < nslabs; i++)
 492                 {
 493                     er->mpi_inbuf[count++] = erg->slab_torque_v[i];
 494                 }
 495             }
 496         }
 497         if (count > er->mpi_bufsize)
 498         {
 499             gmx_fatal(FARGS, "%s MPI buffer overflow, please report this error.", RotStr.c_str());
 500         }
 501
 502 #if GMX_MPI
 503         MPI_Reduce(er->mpi_inbuf, er->mpi_outbuf, count, GMX_MPI_REAL, MPI_SUM, MASTERRANK(cr), cr->mpi_comm_mygroup);
 504 #endif
 505
 506         /* Copy back the reduced data from the buffer on the master */
 507         if (MASTER(cr))
 508         {
 509             count = 0;
 510             for (auto& ergRef : er->enfrotgrp)
 511             {
 512                 gmx_enfrotgrp*  erg  = &ergRef;
 513                 const t_rotgrp* rotg = erg->rotg;
 514                 nslabs               = erg->slab_last - erg->slab_first + 1;
 515                 erg->V               = er->mpi_outbuf[count++];
 516                 erg->torque_v        = er->mpi_outbuf[count++];
 517                 erg->angle_v         = er->mpi_outbuf[count++];
 518                 erg->weight_v        = er->mpi_outbuf[count++];
 519
 520                 if (bPotAngle(er, rotg, step))
 521                 {
 522                     for (int i = 0; i < rotg->PotAngle_nstep; i++)
 523                     {
 524                         erg->PotAngleFit->V[i] = er->mpi_outbuf[count++];
 525                     }
 526                 }
 527                 if (bSlabTau(er, rotg, step))
 528                 {
 529                     for (int i = 0; i < nslabs; i++)
 530                     {
 531                         erg->slab_torque_v[i] = er->mpi_outbuf[count++];
 532                     }
 533                 }
 534             }
 535         }
 536     }
 537
 538     /* Output */
 539     if (MASTER(cr))
 540     {
 541         /* Angle and torque for each rotation group */
 542         for (auto& ergRef : er->enfrotgrp)
 543         {
 544             gmx_enfrotgrp*  erg  = &ergRef;
 545             const t_rotgrp* rotg = erg->rotg;
 546             bFlex                = ISFLEX(rotg);
 547
 548             /* Output to main rotation output file: */
 549             if (do_per_step(step, er->nstrout))
 550             {
 551                 if (RotationGroupFitting::Pot == rotg->eFittype)
 552                 {
 553                     fitangle = get_fitangle(erg);
 554                 }
 555                 else
 556                 {
 557                     if (bFlex)
 558                     {
 559                         fitangle = erg->angle_v; /* RMSD fit angle */
 560                     }
 561                     else
 562                     {
 563                         fitangle = (erg->angle_v / erg->weight_v) * 180.0 * M_1_PI;
 564                     }
 565                 }
 566                 fprintf(er->out_rot, "%12.4f", fitangle);
 567                 fprintf(er->out_rot, "%12.3e", erg->torque_v);
 568                 fprintf(er->out_rot, "%12.3e", erg->V);
 569             }
 570
 571             if (do_per_step(step, er->nstsout))
 572             {
 573                 /* Output to torque log file: */
 574                 if (bFlex)
 575                 {
 576                     fprintf(er->out_torque, "%12.3e%6d", t, erg->groupIndex);
 577                     for (int i = erg->slab_first; i <= erg->slab_last; i++)
 578                     {
 579                         islab = i - erg->slab_first; /* slab index */
 580                         /* Only output if enough weight is in slab */
 581                         if (erg->slab_weights[islab] > rotg->min_gaussian)
 582                         {
 583                             fprintf(er->out_torque, "%6d%12.3e", i, erg->slab_torque_v[islab]);
 584                         }
 585                     }
 586                     fprintf(er->out_torque, "\n");
 587                 }
 588
 589                 /* Output to angles log file: */
 590                 if (RotationGroupFitting::Pot == rotg->eFittype)
 591                 {
 592                     fprintf(er->out_angles, "%12.3e%6d%12.4f", t, erg->groupIndex, erg->degangle);
 593                     /* Output energies at a set of angles around the reference angle */
 594                     for (int i = 0; i < rotg->PotAngle_nstep; i++)
 595                     {
 596                         fprintf(er->out_angles, "%12.3e", erg->PotAngleFit->V[i]);
 597                     }
 598                     fprintf(er->out_angles, "\n");
 599                 }
 600             }
 601         }
 602         if (do_per_step(step, er->nstrout))
 603         {
 604             fprintf(er->out_rot, "\n");
 605         }
 606     }
 607 }
 608
 609
 610 /* Add the forces from enforced rotation potential to the local forces.
 611  * Should be called after the SR forces have been evaluated */
 612 real add_rot_forces(gmx_enfrot* er, gmx::ArrayRef<gmx::RVec> force, const t_commrec* cr, int64_t step, real t)
 613 {
 614     real Vrot = 0.0; /* If more than one rotation group is present, Vrot
 615                         assembles the local parts from all groups         */
 616
 617     /* Loop over enforced rotation groups (usually 1, though)
 618      * Apply the forces from rotation potentials */
 619     for (auto& ergRef : er->enfrotgrp)
 620     {
 621         gmx_enfrotgrp* erg = &ergRef;
 622         Vrot += erg->V; /* add the local parts from the nodes */
 623         const auto& localRotationGroupIndex = erg->atomSet->localIndex();
 624         for (gmx::index l = 0; l < localRotationGroupIndex.ssize(); l++)
 625         {
 626             /* Get the right index of the local force */
 627             int ii = localRotationGroupIndex[l];
 628             /* Add */
 629             rvec_inc(force[ii], erg->f_rot_loc[l]);
 630         }
 631     }
 632
 633     /* Reduce energy,torque, angles etc. to get the sum values (per rotation group)
 634      * on the master and output these values to file. */
 635     if ((do_per_step(step, er->nstrout) || do_per_step(step, er->nstsout)) && er->bOut)
 636     {
 637         reduce_output(cr, er, t, step);
 638     }
 639
 640     /* When appending, er->bOut is FALSE the first time to avoid duplicate entries */
 641     er->bOut = TRUE;
 642
 643     PRINT_POT_TAU
 644
 645     return Vrot;
 646 }
 647
 648
 649 /* The Gaussian norm is chosen such that the sum of the gaussian functions
 650  * over the slabs is approximately 1.0 everywhere */
 651 #define GAUSS_NORM 0.569917543430618
 652
 653
 654 /* Calculate the maximum beta that leads to a gaussian larger min_gaussian,
 655  * also does some checks
 656  */
 657 static double calc_beta_max(real min_gaussian, real slab_dist)
 658 {
 659     double sigma;
 660     double arg;
 661
 662
 663     /* Actually the next two checks are already made in grompp */
 664     if (slab_dist <= 0)
 665     {
 666         gmx_fatal(FARGS, "Slab distance of flexible rotation groups must be >=0 !");
 667     }
 668     if (min_gaussian <= 0)
 669     {
 670         gmx_fatal(FARGS, "Cutoff value for Gaussian must be > 0. (You requested %f)", min_gaussian);
 671     }
 672
 673     /* Define the sigma value */
 674     sigma = 0.7 * slab_dist;
 675
 676     /* Calculate the argument for the logarithm and check that the log() result is negative or 0 */
 677     arg = min_gaussian / GAUSS_NORM;
 678     if (arg > 1.0)
 679     {
 680         gmx_fatal(FARGS, "min_gaussian of flexible rotation groups must be <%g", GAUSS_NORM);
 681     }
 682
 683     return std::sqrt(-2.0 * sigma * sigma * log(min_gaussian / GAUSS_NORM));
 684 }
 685
 686
 687 static inline real calc_beta(rvec curr_x, const gmx_enfrotgrp* erg, int n)
 688 {
 689     return iprod(curr_x, erg->vec) - erg->rotg->slab_dist * n;
 690 }
 691
 692
 693 static inline real gaussian_weight(rvec curr_x, const gmx_enfrotgrp* erg, int n)
 694 {
 695     const real norm = GAUSS_NORM;
 696     real       sigma;
 697
 698
 699     /* Define the sigma value */
 700     sigma = 0.7 * erg->rotg->slab_dist;
 701     /* Calculate the Gaussian value of slab n for position curr_x */
 702     return norm * exp(-0.5 * gmx::square(calc_beta(curr_x, erg, n) / sigma));
 703 }
 704
 705
 706 /* Returns the weight in a single slab, also calculates the Gaussian- and mass-
 707  * weighted sum of positions for that slab */
 708 static real get_slab_weight(int j, const gmx_enfrotgrp* erg, rvec xc[], const real mc[], rvec* x_weighted_sum)
 709 {
 710     rvec curr_x;           /* The position of an atom                      */
 711     rvec curr_x_weighted;  /* The gaussian-weighted position               */
 712     real gaussian;         /* A single gaussian weight                     */
 713     real wgauss;           /* gaussian times current mass                  */
 714     real slabweight = 0.0; /* The sum of weights in the slab               */
 715
 716     clear_rvec(*x_weighted_sum);
 717
 718     /* Loop over all atoms in the rotation group */
 719     for (int i = 0; i < erg->rotg->nat; i++)
 720     {
 721         copy_rvec(xc[i], curr_x);
 722         gaussian = gaussian_weight(curr_x, erg, j);
 723         wgauss   = gaussian * mc[i];
 724         svmul(wgauss, curr_x, curr_x_weighted);
 725         rvec_add(*x_weighted_sum, curr_x_weighted, *x_weighted_sum);
 726         slabweight += wgauss;
 727     } /* END of loop over rotation group atoms */
 728
 729     return slabweight;
 730 }
 731
 732
 733 static void get_slab_centers(gmx_enfrotgrp* erg,  /* Enforced rotation group working data */
 734                              rvec*          xc,   /* The rotation group positions; will
 735                                                      typically be enfrotgrp->xc, but at first call
 736                                                      it is enfrotgrp->xc_ref                      */
 737                              real*    mc,         /* The masses of the rotation group atoms       */
 738                              real     time,       /* Used for output only                         */
 739                              FILE*    out_slabs,  /* For outputting center per slab information   */
 740                              gmx_bool bOutStep,   /* Is this an output step?                      */
 741                              gmx_bool bReference) /* If this routine is called from
 742                                                      init_rot_group we need to store
 743                                                      the reference slab centers                   */
 744 {
 745     /* Loop over slabs */
 746     for (int j = erg->slab_first; j <= erg->slab_last; j++)
 747     {
 748         int slabIndex                = j - erg->slab_first;
 749         erg->slab_weights[slabIndex] = get_slab_weight(j, erg, xc, mc, &erg->slab_center[slabIndex]);
 750
 751         /* We can do the calculations ONLY if there is weight in the slab! */
 752         if (erg->slab_weights[slabIndex] > WEIGHT_MIN)
 753         {
 754             svmul(1.0 / erg->slab_weights[slabIndex], erg->slab_center[slabIndex], erg->slab_center[slabIndex]);
 755         }
 756         else
 757         {
 758             /* We need to check this here, since we divide through slab_weights
 759              * in the flexible low-level routines! */
 760             gmx_fatal(FARGS, "Not enough weight in slab %d. Slab center cannot be determined!", j);
 761         }
 762
 763         /* At first time step: save the centers of the reference structure */
 764         if (bReference)
 765         {
 766             copy_rvec(erg->slab_center[slabIndex], erg->slab_center_ref[slabIndex]);
 767         }
 768     } /* END of loop over slabs */
 769
 770     /* Output on the master */
 771     if ((nullptr != out_slabs) && bOutStep)
 772     {
 773         fprintf(out_slabs, "%12.3e%6d", time, erg->groupIndex);
 774         for (int j = erg->slab_first; j <= erg->slab_last; j++)
 775         {
 776             int slabIndex = j - erg->slab_first;
 777             fprintf(out_slabs,
 778                     "%6d%12.3e%12.3e%12.3e",
 779                     j,
 780                     erg->slab_center[slabIndex][XX],
 781                     erg->slab_center[slabIndex][YY],
 782                     erg->slab_center[slabIndex][ZZ]);
 783         }
 784         fprintf(out_slabs, "\n");
 785     }
 786 }
 787
 788
 789 static void calc_rotmat(const rvec vec,
 790                         real   degangle, /* Angle alpha of rotation at time t in degrees       */
 791                         matrix rotmat)   /* Rotation matrix                                    */
 792 {
 793     real radangle;            /* Rotation angle in radians */
 794     real cosa;                /* cosine alpha              */
 795     real sina;                /* sine alpha                */
 796     real OMcosa;              /* 1 - cos(alpha)            */
 797     real dumxy, dumxz, dumyz; /* save computations         */
 798     rvec rot_vec;             /* Rotate around rot_vec ... */
 799
 800
 801     radangle = degangle * M_PI / 180.0;
 802     copy_rvec(vec, rot_vec);
 803
 804     /* Precompute some variables: */
 805     cosa   = cos(radangle);
 806     sina   = sin(radangle);
 807     OMcosa = 1.0 - cosa;
 808     dumxy  = rot_vec[XX] * rot_vec[YY] * OMcosa;
 809     dumxz  = rot_vec[XX] * rot_vec[ZZ] * OMcosa;
 810     dumyz  = rot_vec[YY] * rot_vec[ZZ] * OMcosa;
 811
 812     /* Construct the rotation matrix for this rotation group: */
 813     /* 1st column: */
 814     rotmat[XX][XX] = cosa + rot_vec[XX] * rot_vec[XX] * OMcosa;
 815     rotmat[YY][XX] = dumxy + rot_vec[ZZ] * sina;
 816     rotmat[ZZ][XX] = dumxz - rot_vec[YY] * sina;
 817     /* 2nd column: */
 818     rotmat[XX][YY] = dumxy - rot_vec[ZZ] * sina;
 819     rotmat[YY][YY] = cosa + rot_vec[YY] * rot_vec[YY] * OMcosa;
 820     rotmat[ZZ][YY] = dumyz + rot_vec[XX] * sina;
 821     /* 3rd column: */
 822     rotmat[XX][ZZ] = dumxz + rot_vec[YY] * sina;
 823     rotmat[YY][ZZ] = dumyz - rot_vec[XX] * sina;
 824     rotmat[ZZ][ZZ] = cosa + rot_vec[ZZ] * rot_vec[ZZ] * OMcosa;
 825
 826 #ifdef PRINTMATRIX
 827     int iii, jjj;
 828
 829     for (iii = 0; iii < 3; iii++)
 830     {
 831         for (jjj = 0; jjj < 3; jjj++)
 832         {
 833             fprintf(stderr, " %10.8f ", rotmat[iii][jjj]);
 834         }
 835         fprintf(stderr, "\n");
 836     }
 837 #endif
 838 }
 839
 840
 841 /* Calculates torque on the rotation axis tau = position x force */
 842 static inline real torque(const rvec rotvec, /* rotation vector; MUST be normalized!          */
 843                           rvec force, /* force                                                */
 844                           rvec x,     /* position of atom on which the force acts             */
 845                           rvec pivot) /* pivot point of rotation axis                         */
 846 {
 847     rvec vectmp, tau;
 848
 849
 850     /* Subtract offset */
 851     rvec_sub(x, pivot, vectmp);
 852
 853     /* position x force */
 854     cprod(vectmp, force, tau);
 855
 856     /* Return the part of the torque which is parallel to the rotation vector */
 857     return iprod(tau, rotvec);
 858 }
 859
 860
 861 /* Right-aligned output of value with standard width */
 862 static void print_aligned(FILE* fp, char const* str)
 863 {
 864     fprintf(fp, "%12s", str);
 865 }
 866
 867
 868 /* Right-aligned output of value with standard short width */
 869 static void print_aligned_short(FILE* fp, char const* str)
 870 {
 871     fprintf(fp, "%6s", str);
 872 }
 873
 874
 875 static FILE* open_output_file(const char* fn, int steps, const char what[])
 876 {
 877     FILE* fp;
 878
 879
 880     fp = gmx_ffopen(fn, "w");
 881
 882     fprintf(fp, "# Output of %s is written in intervals of %d time step%s.\n#\n", what, steps, steps > 1 ? "s" : "");
 883
 884     return fp;
 885 }
 886
 887
 888 /* Open output file for slab center data. Call on master only */
 889 static FILE* open_slab_out(const char* fn, gmx_enfrot* er)
 890 {
 891     FILE* fp;
 892
 893     if (er->restartWithAppending)
 894     {
 895         fp = gmx_fio_fopen(fn, "a");
 896     }
 897     else
 898     {
 899         fp = open_output_file(fn, er->nstsout, "gaussian weighted slab centers");
 900
 901         for (auto& ergRef : er->enfrotgrp)
 902         {
 903             gmx_enfrotgrp* erg = &ergRef;
 904             if (ISFLEX(erg->rotg))
 905             {
 906                 fprintf(fp,
 907                         "# Rotation group %d (%s), slab distance %f nm, %s.\n",
 908                         erg->groupIndex,
 909                         enumValueToString(erg->rotg->eType),
 910                         erg->rotg->slab_dist,
 911                         erg->rotg->bMassW ? "centers of mass" : "geometrical centers");
 912             }
 913         }
 914
 915         fprintf(fp, "# Reference centers are listed first (t=-1).\n");
 916         fprintf(fp, "# The following columns have the syntax:\n");
 917         fprintf(fp, "#     ");
 918         print_aligned_short(fp, "t");
 919         print_aligned_short(fp, "grp");
 920         /* Print legend for the first two entries only ... */
 921         for (int i = 0; i < 2; i++)
 922         {
 923             print_aligned_short(fp, "slab");
 924             print_aligned(fp, "X center");
 925             print_aligned(fp, "Y center");
 926             print_aligned(fp, "Z center");
 927         }
 928         fprintf(fp, " ...\n");
 929         fflush(fp);
 930     }
 931
 932     return fp;
 933 }
 934
 935
 936 /* Adds 'buf' to 'str' */
 937 static void add_to_string(char** str, char* buf)
 938 {
 939     int len;
 940
 941
 942     len = strlen(*str) + strlen(buf) + 1;
 943     srenew(*str, len);
 944     strcat(*str, buf);
 945 }
 946
 947
 948 static void add_to_string_aligned(char** str, char* buf)
 949 {
 950     char buf_aligned[STRLEN];
 951
 952     sprintf(buf_aligned, "%12s", buf);
 953     add_to_string(str, buf_aligned);
 954 }
 955
 956
 957 /* Open output file and print some general information about the rotation groups.
 958  * Call on master only */
 959 static FILE* open_rot_out(const char* fn, const gmx_output_env_t* oenv, gmx_enfrot* er)
 960 {
 961     FILE*        fp;
 962     int          nsets;
 963     const char** setname;
 964     char         buf[50], buf2[75];
 965     gmx_bool     bFlex;
 966     char*        LegendStr = nullptr;
 967     const t_rot* rot       = er->rot;
 968
 969     if (er->restartWithAppending)
 970     {
 971         fp = gmx_fio_fopen(fn, "a");
 972     }
 973     else
 974     {
 975         fp = xvgropen(fn,
 976                       "Rotation angles and energy",
 977                       "Time (ps)",
 978                       "angles (degrees) and energies (kJ/mol)",
 979                       oenv);
 980         fprintf(fp,
 981                 "# Output of enforced rotation data is written in intervals of %d time "
 982                 "step%s.\n#\n",
 983                 er->nstrout,
 984                 er->nstrout > 1 ? "s" : "");
 985         fprintf(fp,
 986                 "# The scalar tau is the torque (kJ/mol) in the direction of the rotation vector "
 987                 "v.\n");
 988         fprintf(fp, "# To obtain the vectorial torque, multiply tau with the group's rot-vec.\n");
 989         fprintf(fp,
 990                 "# For flexible groups, tau(t,n) from all slabs n have been summed in a single "
 991                 "value tau(t) here.\n");
 992         fprintf(fp, "# The torques tau(t,n) are found in the rottorque.log (-rt) output file\n");
 993
 994         for (int g = 0; g < rot->ngrp; g++)
 995         {
 996             const t_rotgrp*      rotg = &rot->grp[g];
 997             const gmx_enfrotgrp* erg  = &er->enfrotgrp[g];
 998             bFlex                     = ISFLEX(rotg);
 999
1000             fprintf(fp, "#\n");
1001             fprintf(fp, "# ROTATION GROUP %d, potential type '%s':\n", g, enumValueToString(rotg->eType));
1002             fprintf(fp, "# rot-massw%d          %s\n", g, booleanValueToString(rotg->bMassW));
1003             fprintf(fp,
1004                     "# rot-vec%d            %12.5e %12.5e %12.5e\n",
1005                     g,
1006                     erg->vec[XX],
1007                     erg->vec[YY],
1008                     erg->vec[ZZ]);
1009             fprintf(fp, "# rot-rate%d           %12.5e degrees/ps\n", g, rotg->rate);
1010             fprintf(fp, "# rot-k%d              %12.5e kJ/(mol*nm^2)\n", g, rotg->k);
1011             if (rotg->eType == EnforcedRotationGroupType::Iso || rotg->eType == EnforcedRotationGroupType::Pm
1012                 || rotg->eType == EnforcedRotationGroupType::Rm
1013                 || rotg->eType == EnforcedRotationGroupType::Rm2)
1014             {
1015                 fprintf(fp,
1016                         "# rot-pivot%d          %12.5e %12.5e %12.5e  nm\n",
1017                         g,
1018                         rotg->pivot[XX],
1019                         rotg->pivot[YY],
1020                         rotg->pivot[ZZ]);
1021             }
1022
1023             if (bFlex)
1024             {
1025                 fprintf(fp, "# rot-slab-distance%d   %f nm\n", g, rotg->slab_dist);
1026                 fprintf(fp, "# rot-min-gaussian%d   %12.5e\n", g, rotg->min_gaussian);
1027             }
1028
1029             /* Output the centers of the rotation groups for the pivot-free potentials */
1030             if ((rotg->eType == EnforcedRotationGroupType::Isopf)
1031                 || (rotg->eType == EnforcedRotationGroupType::Pmpf)
1032                 || (rotg->eType == EnforcedRotationGroupType::Rmpf)
1033                 || (rotg->eType == EnforcedRotationGroupType::Rm2pf
1034                     || (rotg->eType == EnforcedRotationGroupType::Flext)
1035                     || (rotg->eType == EnforcedRotationGroupType::Flex2t)))
1036             {
1037                 fprintf(fp,
1038                         "# ref. grp. %d center  %12.5e %12.5e %12.5e\n",
1039                         g,
1040                         erg->xc_ref_center[XX],
1041                         erg->xc_ref_center[YY],
1042                         erg->xc_ref_center[ZZ]);
1043
1044                 fprintf(fp,
1045                         "# grp. %d init.center  %12.5e %12.5e %12.5e\n",
1046                         g,
1047                         erg->xc_center[XX],
1048                         erg->xc_center[YY],
1049                         erg->xc_center[ZZ]);
1050             }
1051
1052             if ((rotg->eType == EnforcedRotationGroupType::Rm2)
1053                 || (rotg->eType == EnforcedRotationGroupType::Flex2)
1054                 || (rotg->eType == EnforcedRotationGroupType::Flex2t))
1055             {
1056                 fprintf(fp, "# rot-eps%d            %12.5e nm^2\n", g, rotg->eps);
1057             }
1058             if (RotationGroupFitting::Pot == rotg->eFittype)
1059             {
1060                 fprintf(fp, "#\n");
1061                 fprintf(fp,
1062                         "# theta_fit%d is determined by first evaluating the potential for %d "
1063                         "angles around theta_ref%d.\n",
1064                         g,
1065                         rotg->PotAngle_nstep,
1066                         g);
1067                 fprintf(fp,
1068                         "# The fit angle is the one with the smallest potential. It is given as "
1069                         "the deviation\n");
1070                 fprintf(fp,
1071                         "# from the reference angle, i.e. if theta_ref=X and theta_fit=Y, then the "
1072                         "angle with\n");
1073                 fprintf(fp,
1074                         "# minimal value of the potential is X+Y. Angular resolution is %g "
1075                         "degrees.\n",
1076                         rotg->PotAngle_step);
1077             }
1078         }
1079
1080         /* Print a nice legend */
1081         snew(LegendStr, 1);
1082         LegendStr[0] = '\0';
1083         sprintf(buf, "#     %6s", "time");
1084         add_to_string_aligned(&LegendStr, buf);
1085
1086         nsets = 0;
1087         snew(setname, 4 * rot->ngrp);
1088
1089         for (int g = 0; g < rot->ngrp; g++)
1090         {
1091             sprintf(buf, "theta_ref%d", g);
1092             add_to_string_aligned(&LegendStr, buf);
1093
1094             sprintf(buf2, "%s (degrees)", buf);
1095             setname[nsets] = gmx_strdup(buf2);
1096             nsets++;
1097         }
1098         for (int g = 0; g < rot->ngrp; g++)
1099         {
1100             const t_rotgrp* rotg = &rot->grp[g];
1101             bFlex                = ISFLEX(rotg);
1102
1103             /* For flexible axis rotation we use RMSD fitting to determine the
1104              * actual angle of the rotation group */
1105             if (bFlex || RotationGroupFitting::Pot == rotg->eFittype)
1106             {
1107                 sprintf(buf, "theta_fit%d", g);
1108             }
1109             else
1110             {
1111                 sprintf(buf, "theta_av%d", g);
1112             }
1113             add_to_string_aligned(&LegendStr, buf);
1114             sprintf(buf2, "%s (degrees)", buf);
1115             setname[nsets] = gmx_strdup(buf2);
1116             nsets++;
1117
1118             sprintf(buf, "tau%d", g);
1119             add_to_string_aligned(&LegendStr, buf);
1120             sprintf(buf2, "%s (kJ/mol)", buf);
1121             setname[nsets] = gmx_strdup(buf2);
1122             nsets++;
1123
1124             sprintf(buf, "energy%d", g);
1125             add_to_string_aligned(&LegendStr, buf);
1126             sprintf(buf2, "%s (kJ/mol)", buf);
1127             setname[nsets] = gmx_strdup(buf2);
1128             nsets++;
1129         }
1130         fprintf(fp, "#\n");
1131
1132         if (nsets > 1)
1133         {
1134             xvgr_legend(fp, nsets, setname, oenv);
1135         }
1136         sfree(setname);
1137
1138         fprintf(fp, "#\n# Legend for the following data columns:\n");
1139         fprintf(fp, "%s\n", LegendStr);
1140         sfree(LegendStr);
1141
1142         fflush(fp);
1143     }
1144
1145     return fp;
1146 }
1147
1148
1149 /* Call on master only */
1150 static FILE* open_angles_out(const char* fn, gmx_enfrot* er)
1151 {
1152     FILE*        fp;
1153     char         buf[100];
1154     const t_rot* rot = er->rot;
1155
1156     if (er->restartWithAppending)
1157     {
1158         fp = gmx_fio_fopen(fn, "a");
1159     }
1160     else
1161     {
1162         /* Open output file and write some information about it's structure: */
1163         fp = open_output_file(fn, er->nstsout, "rotation group angles");
1164         fprintf(fp, "# All angles given in degrees, time in ps.\n");
1165         for (int g = 0; g < rot->ngrp; g++)
1166         {
1167             const t_rotgrp*      rotg = &rot->grp[g];
1168             const gmx_enfrotgrp* erg  = &er->enfrotgrp[g];
1169
1170             /* Output for this group happens only if potential type is flexible or
1171              * if fit type is potential! */
1172             if (ISFLEX(rotg) || (RotationGroupFitting::Pot == rotg->eFittype))
1173             {
1174                 if (ISFLEX(rotg))
1175                 {
1176                     sprintf(buf, " slab distance %f nm, ", rotg->slab_dist);
1177                 }
1178                 else
1179                 {
1180                     buf[0] = '\0';
1181                 }
1182
1183                 fprintf(fp,
1184                         "#\n# ROTATION GROUP %d '%s',%s fit type '%s'.\n",
1185                         g,
1186                         enumValueToString(rotg->eType),
1187                         buf,
1188                         enumValueToString(rotg->eFittype));
1189
1190                 /* Special type of fitting using the potential minimum. This is
1191                  * done for the whole group only, not for the individual slabs. */
1192                 if (RotationGroupFitting::Pot == rotg->eFittype)
1193                 {
1194                     fprintf(fp,
1195                             "#    To obtain theta_fit%d, the potential is evaluated for %d angles "
1196                             "around theta_ref%d\n",
1197                             g,
1198                             rotg->PotAngle_nstep,
1199                             g);
1200                     fprintf(fp,
1201                             "#    The fit angle in the rotation standard outfile is the one with "
1202                             "minimal energy E(theta_fit) [kJ/mol].\n");
1203                     fprintf(fp, "#\n");
1204                 }
1205
1206                 fprintf(fp, "# Legend for the group %d data columns:\n", g);
1207                 fprintf(fp, "#     ");
1208                 print_aligned_short(fp, "time");
1209                 print_aligned_short(fp, "grp");
1210                 print_aligned(fp, "theta_ref");
1211
1212                 if (RotationGroupFitting::Pot == rotg->eFittype)
1213                 {
1214                     /* Output the set of angles around the reference angle */
1215                     for (int i = 0; i < rotg->PotAngle_nstep; i++)
1216                     {
1217                         sprintf(buf, "E(%g)", erg->PotAngleFit->degangle[i]);
1218                         print_aligned(fp, buf);
1219                     }
1220                 }
1221                 else
1222                 {
1223                     /* Output fit angle for each slab */
1224                     print_aligned_short(fp, "slab");
1225                     print_aligned_short(fp, "atoms");
1226                     print_aligned(fp, "theta_fit");
1227                     print_aligned_short(fp, "slab");
1228                     print_aligned_short(fp, "atoms");
1229                     print_aligned(fp, "theta_fit");
1230                     fprintf(fp, " ...");
1231                 }
1232                 fprintf(fp, "\n");
1233             }
1234         }
1235         fflush(fp);
1236     }
1237
1238     return fp;
1239 }
1240
1241
1242 /* Open torque output file and write some information about it's structure.
1243  * Call on master only */
1244 static FILE* open_torque_out(const char* fn, gmx_enfrot* er)
1245 {
1246     FILE*        fp;
1247     const t_rot* rot = er->rot;
1248
1249     if (er->restartWithAppending)
1250     {
1251         fp = gmx_fio_fopen(fn, "a");
1252     }
1253     else
1254     {
1255         fp = open_output_file(fn, er->nstsout, "torques");
1256
1257         for (int g = 0; g < rot->ngrp; g++)
1258         {
1259             const t_rotgrp*      rotg = &rot->grp[g];
1260             const gmx_enfrotgrp* erg  = &er->enfrotgrp[g];
1261             if (ISFLEX(rotg))
1262             {
1263                 fprintf(fp,
1264                         "# Rotation group %d (%s), slab distance %f nm.\n",
1265                         g,
1266                         enumValueToString(rotg->eType),
1267                         rotg->slab_dist);
1268                 fprintf(fp,
1269                         "# The scalar tau is the torque (kJ/mol) in the direction of the rotation "
1270                         "vector.\n");
1271                 fprintf(fp, "# To obtain the vectorial torque, multiply tau with\n");
1272                 fprintf(fp,
1273                         "# rot-vec%d            %10.3e %10.3e %10.3e\n",
1274                         g,
1275                         erg->vec[XX],
1276                         erg->vec[YY],
1277                         erg->vec[ZZ]);
1278                 fprintf(fp, "#\n");
1279             }
1280         }
1281         fprintf(fp, "# Legend for the following data columns: (tau=torque for that slab):\n");
1282         fprintf(fp, "#     ");
1283         print_aligned_short(fp, "t");
1284         print_aligned_short(fp, "grp");
1285         print_aligned_short(fp, "slab");
1286         print_aligned(fp, "tau");
1287         print_aligned_short(fp, "slab");
1288         print_aligned(fp, "tau");
1289         fprintf(fp, " ...\n");
1290         fflush(fp);
1291     }
1292
1293     return fp;
1294 }
1295
1296
1297 static void swap_val(double* vec, int i, int j)
1298 {
1299     double tmp = vec[j];
1300
1301
1302     vec[j] = vec[i];
1303     vec[i] = tmp;
1304 }
1305
1306
1307 static void swap_col(double** mat, int i, int j)
1308 {
1309     double tmp[3] = { mat[0][j], mat[1][j], mat[2][j] };
1310
1311
1312     mat[0][j] = mat[0][i];
1313     mat[1][j] = mat[1][i];
1314     mat[2][j] = mat[2][i];
1315
1316     mat[0][i] = tmp[0];
1317     mat[1][i] = tmp[1];
1318     mat[2][i] = tmp[2];
1319 }
1320
1321
1322 /* Eigenvectors are stored in columns of eigen_vec */
1323 static void diagonalize_symmetric(double** matrix, double** eigen_vec, double eigenval[3])
1324 {
1325     int n_rot;
1326
1327
1328     jacobi(matrix, 3, eigenval, eigen_vec, &n_rot);
1329
1330     /* sort in ascending order */
1331     if (eigenval[0] > eigenval[1])
1332     {
1333         swap_val(eigenval, 0, 1);
1334         swap_col(eigen_vec, 0, 1);
1335     }
1336     if (eigenval[1] > eigenval[2])
1337     {
1338         swap_val(eigenval, 1, 2);
1339         swap_col(eigen_vec, 1, 2);
1340     }
1341     if (eigenval[0] > eigenval[1])
1342     {
1343         swap_val(eigenval, 0, 1);
1344         swap_col(eigen_vec, 0, 1);
1345     }
1346 }
1347
1348
1349 static void align_with_z(rvec* s, /* Structure to align */
1350                          int   natoms,
1351                          rvec  axis)
1352 {
1353     int    i, j, k;
1354     rvec   zet         = { 0.0, 0.0, 1.0 };
1355     rvec   rot_axis    = { 0.0, 0.0, 0.0 };
1356     rvec*  rotated_str = nullptr;
1357     real   ooanorm;
1358     real   angle;
1359     matrix rotmat;
1360
1361
1362     snew(rotated_str, natoms);
1363
1364     /* Normalize the axis */
1365     ooanorm = 1.0 / norm(axis);
1366     svmul(ooanorm, axis, axis);
1367
1368     /* Calculate the angle for the fitting procedure */
1369     cprod(axis, zet, rot_axis);
1370     angle = acos(axis[2]);
1371     if (angle < 0.0)
1372     {
1373         angle += M_PI;
1374     }
1375
1376     /* Calculate the rotation matrix */
1377     calc_rotmat(rot_axis, angle * 180.0 / M_PI, rotmat);
1378
1379     /* Apply the rotation matrix to s */
1380     for (i = 0; i < natoms; i++)
1381     {
1382         for (j = 0; j < 3; j++)
1383         {
1384             for (k = 0; k < 3; k++)
1385             {
1386                 rotated_str[i][j] += rotmat[j][k] * s[i][k];
1387             }
1388         }
1389     }
1390
1391     /* Rewrite the rotated structure to s */
1392     for (i = 0; i < natoms; i++)
1393     {
1394         for (j = 0; j < 3; j++)
1395         {
1396             s[i][j] = rotated_str[i][j];
1397         }
1398     }
1399
1400     sfree(rotated_str);
1401 }
1402
1403
1404 static void calc_correl_matrix(rvec* Xstr, rvec* Ystr, double** Rmat, int natoms)
1405 {
1406     int i, j, k;
1407
1408
1409     for (i = 0; i < 3; i++)
1410     {
1411         for (j = 0; j < 3; j++)
1412         {
1413             Rmat[i][j] = 0.0;
1414         }
1415     }
1416
1417     for (i = 0; i < 3; i++)
1418     {
1419         for (j = 0; j < 3; j++)
1420         {
1421             for (k = 0; k < natoms; k++)
1422             {
1423                 Rmat[i][j] += Ystr[k][i] * Xstr[k][j];
1424             }
1425         }
1426     }
1427 }
1428
1429
1430 static void weigh_coords(rvec* str, real* weight, int natoms)
1431 {
1432     int i, j;
1433
1434
1435     for (i = 0; i < natoms; i++)
1436     {
1437         for (j = 0; j < 3; j++)
1438         {
1439             str[i][j] *= std::sqrt(weight[i]);
1440         }
1441     }
1442 }
1443
1444
1445 static real opt_angle_analytic(rvec*      ref_s,
1446                                rvec*      act_s,
1447                                real*      weight,
1448                                int        natoms,
1449                                const rvec ref_com,
1450                                const rvec act_com,
1451                                rvec       axis)
1452 {
1453     int      i, j, k;
1454     rvec*    ref_s_1 = nullptr;
1455     rvec*    act_s_1 = nullptr;
1456     rvec     shift;
1457     double **Rmat, **RtR, **eigvec;
1458     double   eigval[3];
1459     double   V[3][3], WS[3][3];
1460     double   rot_matrix[3][3];
1461     double   opt_angle;
1462
1463
1464     /* Do not change the original coordinates */
1465     snew(ref_s_1, natoms);
1466     snew(act_s_1, natoms);
1467     for (i = 0; i < natoms; i++)
1468     {
1469         copy_rvec(ref_s[i], ref_s_1[i]);
1470         copy_rvec(act_s[i], act_s_1[i]);
1471     }
1472
1473     /* Translate the structures to the origin */
1474     shift[XX] = -ref_com[XX];
1475     shift[YY] = -ref_com[YY];
1476     shift[ZZ] = -ref_com[ZZ];
1477     translate_x(ref_s_1, natoms, shift);
1478
1479     shift[XX] = -act_com[XX];
1480     shift[YY] = -act_com[YY];
1481     shift[ZZ] = -act_com[ZZ];
1482     translate_x(act_s_1, natoms, shift);
1483
1484     /* Align rotation axis with z */
1485     align_with_z(ref_s_1, natoms, axis);
1486     align_with_z(act_s_1, natoms, axis);
1487
1488     /* Correlation matrix */
1489     Rmat = allocate_square_matrix(3);
1490
1491     for (i = 0; i < natoms; i++)
1492     {
1493         ref_s_1[i][2] = 0.0;
1494         act_s_1[i][2] = 0.0;
1495     }
1496
1497     /* Weight positions with sqrt(weight) */
1498     if (nullptr != weight)
1499     {
1500         weigh_coords(ref_s_1, weight, natoms);
1501         weigh_coords(act_s_1, weight, natoms);
1502     }
1503
1504     /* Calculate correlation matrices R=YXt (X=ref_s; Y=act_s) */
1505     calc_correl_matrix(ref_s_1, act_s_1, Rmat, natoms);
1506
1507     /* Calculate RtR */
1508     RtR = allocate_square_matrix(3);
1509     for (i = 0; i < 3; i++)
1510     {
1511         for (j = 0; j < 3; j++)
1512         {
1513             for (k = 0; k < 3; k++)
1514             {
1515                 RtR[i][j] += Rmat[k][i] * Rmat[k][j];
1516             }
1517         }
1518     }
1519     /* Diagonalize RtR */
1520     snew(eigvec, 3);
1521     for (i = 0; i < 3; i++)
1522     {
1523         snew(eigvec[i], 3);
1524     }
1525
1526     diagonalize_symmetric(RtR, eigvec, eigval);
1527     swap_col(eigvec, 0, 1);
1528     swap_col(eigvec, 1, 2);
1529     swap_val(eigval, 0, 1);
1530     swap_val(eigval, 1, 2);
1531
1532     /* Calculate V */
1533     for (i = 0; i < 3; i++)
1534     {
1535         for (j = 0; j < 3; j++)
1536         {
1537             V[i][j]  = 0.0;
1538             WS[i][j] = 0.0;
1539         }
1540     }
1541
1542     for (i = 0; i < 2; i++)
1543     {
1544         for (j = 0; j < 2; j++)
1545         {
1546             WS[i][j] = eigvec[i][j] / std::sqrt(eigval[j]);
1547         }
1548     }
1549
1550     for (i = 0; i < 3; i++)
1551     {
1552         for (j = 0; j < 3; j++)
1553         {
1554             for (k = 0; k < 3; k++)
1555             {
1556                 V[i][j] += Rmat[i][k] * WS[k][j];
1557             }
1558         }
1559     }
1560     free_square_matrix(Rmat, 3);
1561
1562     /* Calculate optimal rotation matrix */
1563     for (i = 0; i < 3; i++)
1564     {
1565         for (j = 0; j < 3; j++)
1566         {
1567             rot_matrix[i][j] = 0.0;
1568         }
1569     }
1570
1571     for (i = 0; i < 3; i++)
1572     {
1573         for (j = 0; j < 3; j++)
1574         {
1575             for (k = 0; k < 3; k++)
1576             {
1577                 rot_matrix[i][j] += eigvec[i][k] * V[j][k];
1578             }
1579         }
1580     }
1581     rot_matrix[2][2] = 1.0;
1582
1583     /* In some cases abs(rot_matrix[0][0]) can be slighly larger
1584      * than unity due to numerical inacurracies. To be able to calculate
1585      * the acos function, we put these values back in range. */
1586     if (rot_matrix[0][0] > 1.0)
1587     {
1588         rot_matrix[0][0] = 1.0;
1589     }
1590     else if (rot_matrix[0][0] < -1.0)
1591     {
1592         rot_matrix[0][0] = -1.0;
1593     }
1594
1595     /* Determine the optimal rotation angle: */
1596     opt_angle = (-1.0) * acos(rot_matrix[0][0]) * 180.0 / M_PI;
1597     if (rot_matrix[0][1] < 0.0)
1598     {
1599         opt_angle = (-1.0) * opt_angle;
1600     }
1601
1602     /* Give back some memory */
1603     free_square_matrix(RtR, 3);
1604     sfree(ref_s_1);
1605     sfree(act_s_1);
1606     for (i = 0; i < 3; i++)
1607     {
1608         sfree(eigvec[i]);
1609     }
1610     sfree(eigvec);
1611
1612     return static_cast<real>(opt_angle);
1613 }
1614
1615
1616 /* Determine angle of the group by RMSD fit to the reference */
1617 /* Not parallelized, call this routine only on the master */
1618 static real flex_fit_angle(gmx_enfrotgrp* erg)
1619 {
1620     rvec* fitcoords = nullptr;
1621     rvec  center;   /* Center of positions passed to the fit routine */
1622     real  fitangle; /* Angle of the rotation group derived by fitting */
1623     rvec  coord;
1624     real  scal;
1625
1626     /* Get the center of the rotation group.
1627      * Note, again, erg->xc has been sorted in do_flexible */
1628     get_center(erg->xc, erg->mc_sorted, erg->rotg->nat, center);
1629
1630     /* === Determine the optimal fit angle for the rotation group === */
1631     if (erg->rotg->eFittype == RotationGroupFitting::Norm)
1632     {
1633         /* Normalize every position to it's reference length */
1634         for (int i = 0; i < erg->rotg->nat; i++)
1635         {
1636             /* Put the center of the positions into the origin */
1637             rvec_sub(erg->xc[i], center, coord);
1638             /* Determine the scaling factor for the length: */
1639             scal = erg->xc_ref_length[erg->xc_sortind[i]] / norm(coord);
1640             /* Get position, multiply with the scaling factor and save  */
1641             svmul(scal, coord, erg->xc_norm[i]);
1642         }
1643         fitcoords = erg->xc_norm;
1644     }
1645     else
1646     {
1647         fitcoords = erg->xc;
1648     }
1649     /* From the point of view of the current positions, the reference has rotated
1650      * backwards. Since we output the angle relative to the fixed reference,
1651      * we need the minus sign. */
1652     fitangle = -opt_angle_analytic(
1653             erg->xc_ref_sorted, fitcoords, erg->mc_sorted, erg->rotg->nat, erg->xc_ref_center, center, erg->vec);
1654
1655     return fitangle;
1656 }
1657
1658
1659 /* Determine actual angle of each slab by RMSD fit to the reference */
1660 /* Not parallelized, call this routine only on the master */
1661 static void flex_fit_angle_perslab(gmx_enfrotgrp* erg, double t, real degangle, FILE* fp)
1662 {
1663     rvec curr_x, ref_x;
1664     rvec act_center; /* Center of actual positions that are passed to the fit routine */
1665     rvec ref_center; /* Same for the reference positions */
1666     real fitangle;   /* Angle of a slab derived from an RMSD fit to
1667                       * the reference structure at t=0  */
1668     gmx_slabdata* sd;
1669     real          OOm_av; /* 1/average_mass of a rotation group atom */
1670     real          m_rel;  /* Relative mass of a rotation group atom  */
1671
1672
1673     /* Average mass of a rotation group atom: */
1674     OOm_av = erg->invmass * erg->rotg->nat;
1675
1676     /**********************************/
1677     /* First collect the data we need */
1678     /**********************************/
1679
1680     /* Collect the data for the individual slabs */
1681     for (int n = erg->slab_first; n <= erg->slab_last; n++)
1682     {
1683         int slabIndex = n - erg->slab_first; /* slab index */
1684         sd            = &(erg->slab_data[slabIndex]);
1685         sd->nat       = erg->lastatom[slabIndex] - erg->firstatom[slabIndex] + 1;
1686         int ind       = 0;
1687
1688         /* Loop over the relevant atoms in the slab */
1689         for (int l = erg->firstatom[slabIndex]; l <= erg->lastatom[slabIndex]; l++)
1690         {
1691             /* Current position of this atom: x[ii][XX/YY/ZZ] */
1692             copy_rvec(erg->xc[l], curr_x);
1693
1694             /* The (unrotated) reference position of this atom is copied to ref_x.
1695              * Beware, the xc coords have been sorted in do_flexible */
1696             copy_rvec(erg->xc_ref_sorted[l], ref_x);
1697
1698             /* Save data for doing angular RMSD fit later */
1699             /* Save the current atom position */
1700             copy_rvec(curr_x, sd->x[ind]);
1701             /* Save the corresponding reference position */
1702             copy_rvec(ref_x, sd->ref[ind]);
1703
1704             /* Maybe also mass-weighting was requested. If yes, additionally
1705              * multiply the weights with the relative mass of the atom. If not,
1706              * multiply with unity. */
1707             m_rel = erg->mc_sorted[l] * OOm_av;
1708
1709             /* Save the weight for this atom in this slab */
1710             sd->weight[ind] = gaussian_weight(curr_x, erg, n) * m_rel;
1711
1712             /* Next atom in this slab */
1713             ind++;
1714         }
1715     }
1716
1717     /******************************/
1718     /* Now do the fit calculation */
1719     /******************************/
1720
1721     fprintf(fp, "%12.3e%6d%12.3f", t, erg->groupIndex, degangle);
1722
1723     /* === Now do RMSD fitting for each slab === */
1724     /* We require at least SLAB_MIN_ATOMS in a slab, such that the fit makes sense. */
1725 #define SLAB_MIN_ATOMS 4
1726
1727     for (int n = erg->slab_first; n <= erg->slab_last; n++)
1728     {
1729         int slabIndex = n - erg->slab_first; /* slab index */
1730         sd            = &(erg->slab_data[slabIndex]);
1731         if (sd->nat >= SLAB_MIN_ATOMS)
1732         {
1733             /* Get the center of the slabs reference and current positions */
1734             get_center(sd->ref, sd->weight, sd->nat, ref_center);
1735             get_center(sd->x, sd->weight, sd->nat, act_center);
1736             if (erg->rotg->eFittype == RotationGroupFitting::Norm)
1737             {
1738                 /* Normalize every position to it's reference length
1739                  * prior to performing the fit */
1740                 for (int i = 0; i < sd->nat; i++) /* Center */
1741                 {
1742                     rvec_dec(sd->ref[i], ref_center);
1743                     rvec_dec(sd->x[i], act_center);
1744                     /* Normalize x_i such that it gets the same length as ref_i */
1745                     svmul(norm(sd->ref[i]) / norm(sd->x[i]), sd->x[i], sd->x[i]);
1746                 }
1747                 /* We already subtracted the centers */
1748                 clear_rvec(ref_center);
1749                 clear_rvec(act_center);
1750             }
1751             fitangle = -opt_angle_analytic(
1752                     sd->ref, sd->x, sd->weight, sd->nat, ref_center, act_center, erg->vec);
1753             fprintf(fp, "%6d%6d%12.3f", n, sd->nat, fitangle);
1754         }
1755     }
1756     fprintf(fp, "\n");
1757
1758 #undef SLAB_MIN_ATOMS
1759 }
1760
1761
1762 /* Shift x with is */
1763 static inline void shift_single_coord(const matrix box, rvec x, const ivec is)
1764 {
1765     int tx, ty, tz;
1766
1767
1768     tx = is[XX];
1769     ty = is[YY];
1770     tz = is[ZZ];
1771
1772     if (TRICLINIC(box))
1773     {
1774         x[XX] += tx * box[XX][XX] + ty * box[YY][XX] + tz * box[ZZ][XX];
1775         x[YY] += ty * box[YY][YY] + tz * box[ZZ][YY];
1776         x[ZZ] += tz * box[ZZ][ZZ];
1777     }
1778     else
1779     {
1780         x[XX] += tx * box[XX][XX];
1781         x[YY] += ty * box[YY][YY];
1782         x[ZZ] += tz * box[ZZ][ZZ];
1783     }
1784 }
1785
1786
1787 /* Determine the 'home' slab of this atom which is the
1788  * slab with the highest Gaussian weight of all */
1789 static inline int get_homeslab(rvec curr_x, /* The position for which the home slab shall be determined */
1790                                const rvec rotvec, /* The rotation vector */
1791                                real       slabdist)     /* The slab distance */
1792 {
1793     real dist;
1794
1795
1796     /* The distance of the atom to the coordinate center (where the
1797      * slab with index 0) is */
1798     dist = iprod(rotvec, curr_x);
1799
1800     return gmx::roundToInt(dist / slabdist);
1801 }
1802
1803
1804 /* For a local atom determine the relevant slabs, i.e. slabs in
1805  * which the gaussian is larger than min_gaussian
1806  */
1807 static int get_single_atom_gaussians(rvec curr_x, gmx_enfrotgrp* erg)
1808 {
1809
1810     /* Determine the 'home' slab of this atom: */
1811     int homeslab = get_homeslab(curr_x, erg->vec, erg->rotg->slab_dist);
1812
1813     /* First determine the weight in the atoms home slab: */
1814     real g                 = gaussian_weight(curr_x, erg, homeslab);
1815     int  count             = 0;
1816     erg->gn_atom[count]    = g;
1817     erg->gn_slabind[count] = homeslab;
1818     count++;
1819
1820
1821     /* Determine the max slab */
1822     int slab = homeslab;
1823     while (g > erg->rotg->min_gaussian)
1824     {
1825         slab++;
1826         g                      = gaussian_weight(curr_x, erg, slab);
1827         erg->gn_slabind[count] = slab;
1828         erg->gn_atom[count]    = g;
1829         count++;
1830     }
1831     count--;
1832
1833     /* Determine the min slab */
1834     slab = homeslab;
1835     do
1836     {
1837         slab--;
1838         g                      = gaussian_weight(curr_x, erg, slab);
1839         erg->gn_slabind[count] = slab;
1840         erg->gn_atom[count]    = g;
1841         count++;
1842     } while (g > erg->rotg->min_gaussian);
1843     count--;
1844
1845     return count;
1846 }
1847
1848
1849 static void flex2_precalc_inner_sum(const gmx_enfrotgrp* erg)
1850 {
1851     rvec xi;       /* positions in the i-sum                        */
1852     rvec xcn, ycn; /* the current and the reference slab centers    */
1853     real gaussian_xi;
1854     rvec yi0;
1855     rvec rin; /* Helper variables                              */
1856     real fac, fac2;
1857     rvec innersumvec;
1858     real OOpsii, OOpsiistar;
1859     real sin_rin; /* s_ii.r_ii */
1860     rvec s_in, tmpvec, tmpvec2;
1861     real mi, wi; /* Mass-weighting of the positions                 */
1862     real N_M;    /* N/M                                             */
1863
1864
1865     N_M = erg->rotg->nat * erg->invmass;
1866
1867     /* Loop over all slabs that contain something */
1868     for (int n = erg->slab_first; n <= erg->slab_last; n++)
1869     {
1870         int slabIndex = n - erg->slab_first; /* slab index */
1871
1872         /* The current center of this slab is saved in xcn: */
1873         copy_rvec(erg->slab_center[slabIndex], xcn);
1874         /* ... and the reference center in ycn: */
1875         copy_rvec(erg->slab_center_ref[slabIndex + erg->slab_buffer], ycn);
1876
1877         /*** D. Calculate the whole inner sum used for second and third sum */
1878         /* For slab n, we need to loop over all atoms i again. Since we sorted
1879          * the atoms with respect to the rotation vector, we know that it is sufficient
1880          * to calculate from firstatom to lastatom only. All other contributions will
1881          * be very small. */
1882         clear_rvec(innersumvec);
1883         for (int i = erg->firstatom[slabIndex]; i <= erg->lastatom[slabIndex]; i++)
1884         {
1885             /* Coordinate xi of this atom */
1886             copy_rvec(erg->xc[i], xi);
1887
1888             /* The i-weights */
1889             gaussian_xi = gaussian_weight(xi, erg, n);
1890             mi          = erg->mc_sorted[i]; /* need the sorted mass here */
1891             wi          = N_M * mi;
1892
1893             /* Calculate rin */
1894             copy_rvec(erg->xc_ref_sorted[i], yi0); /* Reference position yi0   */
1895             rvec_sub(yi0, ycn, tmpvec2);           /* tmpvec2 = yi0 - ycn      */
1896             mvmul(erg->rotmat, tmpvec2, rin);      /* rin = Omega.(yi0 - ycn)  */
1897
1898             /* Calculate psi_i* and sin */
1899             rvec_sub(xi, xcn, tmpvec2); /* tmpvec2 = xi - xcn       */
1900
1901             /* In rare cases, when an atom position coincides with a slab center
1902              * (tmpvec2 == 0) we cannot compute the vector product for s_in.
1903              * However, since the atom is located directly on the pivot, this
1904              * slab's contribution to the force on that atom will be zero
1905              * anyway. Therefore, we continue with the next atom. */
1906             if (gmx_numzero(norm(tmpvec2))) /* 0 == norm(xi - xcn) */
1907             {
1908                 continue;
1909             }
1910
1911             cprod(erg->vec, tmpvec2, tmpvec);            /* tmpvec = v x (xi - xcn)  */
1912             OOpsiistar = norm2(tmpvec) + erg->rotg->eps; /* OOpsii* = 1/psii* = |v x (xi-xcn)|^2 + eps */
1913             OOpsii = norm(tmpvec);                       /* OOpsii = 1 / psii = |v x (xi - xcn)| */
1914
1915             /*                           *         v x (xi - xcn)          */
1916             unitv(tmpvec, s_in); /*  sin = ----------------         */
1917                                  /*        |v x (xi - xcn)|         */
1918
1919             sin_rin = iprod(s_in, rin); /* sin_rin = sin . rin             */
1920
1921             /* Now the whole sum */
1922             fac = OOpsii / OOpsiistar;
1923             svmul(fac, rin, tmpvec);
1924             fac2 = fac * fac * OOpsii;
1925             svmul(fac2 * sin_rin, s_in, tmpvec2);
1926             rvec_dec(tmpvec, tmpvec2);
1927
1928             svmul(wi * gaussian_xi * sin_rin, tmpvec, tmpvec2);
1929
1930             rvec_inc(innersumvec, tmpvec2);
1931         } /* now we have the inner sum, used both for sum2 and sum3 */
1932
1933         /* Save it to be used in do_flex2_lowlevel */
1934         copy_rvec(innersumvec, erg->slab_innersumvec[slabIndex]);
1935     } /* END of loop over slabs */
1936 }
1937
1938
1939 static void flex_precalc_inner_sum(const gmx_enfrotgrp* erg)
1940 {
1941     rvec xi;       /* position                                      */
1942     rvec xcn, ycn; /* the current and the reference slab centers    */
1943     rvec qin, rin; /* q_i^n and r_i^n                               */
1944     real bin;
1945     rvec tmpvec;
1946     rvec innersumvec; /* Inner part of sum_n2                          */
1947     real gaussian_xi; /* Gaussian weight gn(xi)                        */
1948     real mi, wi;      /* Mass-weighting of the positions               */
1949     real N_M;         /* N/M                                           */
1950
1951     N_M = erg->rotg->nat * erg->invmass;
1952
1953     /* Loop over all slabs that contain something */
1954     for (int n = erg->slab_first; n <= erg->slab_last; n++)
1955     {
1956         int slabIndex = n - erg->slab_first; /* slab index */
1957
1958         /* The current center of this slab is saved in xcn: */
1959         copy_rvec(erg->slab_center[slabIndex], xcn);
1960         /* ... and the reference center in ycn: */
1961         copy_rvec(erg->slab_center_ref[slabIndex + erg->slab_buffer], ycn);
1962
1963         /* For slab n, we need to loop over all atoms i again. Since we sorted
1964          * the atoms with respect to the rotation vector, we know that it is sufficient
1965          * to calculate from firstatom to lastatom only. All other contributions will
1966          * be very small. */
1967         clear_rvec(innersumvec);
1968         for (int i = erg->firstatom[slabIndex]; i <= erg->lastatom[slabIndex]; i++)
1969         {
1970             /* Coordinate xi of this atom */
1971             copy_rvec(erg->xc[i], xi);
1972
1973             /* The i-weights */
1974             gaussian_xi = gaussian_weight(xi, erg, n);
1975             mi          = erg->mc_sorted[i]; /* need the sorted mass here */
1976             wi          = N_M * mi;
1977
1978             /* Calculate rin and qin */
1979             rvec_sub(erg->xc_ref_sorted[i], ycn, tmpvec); /* tmpvec = yi0-ycn */
1980
1981             /* In rare cases, when an atom position coincides with a slab center
1982              * (tmpvec == 0) we cannot compute the vector product for qin.
1983              * However, since the atom is located directly on the pivot, this
1984              * slab's contribution to the force on that atom will be zero
1985              * anyway. Therefore, we continue with the next atom. */
1986             if (gmx_numzero(norm(tmpvec))) /* 0 == norm(yi0 - ycn) */
1987             {
1988                 continue;
1989             }
1990
1991             mvmul(erg->rotmat, tmpvec, rin); /* rin = Omega.(yi0 - ycn)  */
1992             cprod(erg->vec, rin, tmpvec);    /* tmpvec = v x Omega*(yi0-ycn) */
1993
1994             /*                                *        v x Omega*(yi0-ycn)    */
1995             unitv(tmpvec, qin); /* qin = ---------------------   */
1996                                 /*       |v x Omega*(yi0-ycn)|   */
1997
1998             /* Calculate bin */
1999             rvec_sub(xi, xcn, tmpvec); /* tmpvec = xi-xcn          */
2000             bin = iprod(qin, tmpvec);  /* bin  = qin*(xi-xcn)      */
2001
2002             svmul(wi * gaussian_xi * bin, qin, tmpvec);
2003
2004             /* Add this contribution to the inner sum: */
2005             rvec_add(innersumvec, tmpvec, innersumvec);
2006         } /* now we have the inner sum vector S^n for this slab */
2007           /* Save it to be used in do_flex_lowlevel */
2008         copy_rvec(innersumvec, erg->slab_innersumvec[slabIndex]);
2009     }
2010 }
2011
2012
2013 static real do_flex2_lowlevel(gmx_enfrotgrp*                 erg,
2014                               real                           sigma, /* The Gaussian width sigma */
2015                               gmx::ArrayRef<const gmx::RVec> coords,
2016                               gmx_bool                       bOutstepRot,
2017                               gmx_bool                       bOutstepSlab,
2018                               const matrix                   box)
2019 {
2020     int  count, ii, iigrp;
2021     rvec xj;          /* position in the i-sum                         */
2022     rvec yj0;         /* the reference position in the j-sum           */
2023     rvec xcn, ycn;    /* the current and the reference slab centers    */
2024     real V;           /* This node's part of the rotation pot. energy  */
2025     real gaussian_xj; /* Gaussian weight                               */
2026     real beta;
2027
2028     real numerator, fit_numerator;
2029     rvec rjn, fit_rjn; /* Helper variables                              */
2030     real fac, fac2;
2031
2032     real     OOpsij, OOpsijstar;
2033     real     OOsigma2; /* 1/(sigma^2)                                   */
2034     real     sjn_rjn;
2035     real     betasigpsi;
2036     rvec     sjn, tmpvec, tmpvec2, yj0_ycn;
2037     rvec     sum1vec_part, sum1vec, sum2vec_part, sum2vec, sum3vec, sum4vec, innersumvec;
2038     real     sum3, sum4;
2039     real     mj, wj; /* Mass-weighting of the positions               */
2040     real     N_M;    /* N/M                                           */
2041     real     Wjn;    /* g_n(x_j) m_j / Mjn                            */
2042     gmx_bool bCalcPotFit;
2043
2044     /* To calculate the torque per slab */
2045     rvec slab_force; /* Single force from slab n on one atom          */
2046     rvec slab_sum1vec_part;
2047     real slab_sum3part, slab_sum4part;
2048     rvec slab_sum1vec, slab_sum2vec, slab_sum3vec, slab_sum4vec;
2049
2050     /* Pre-calculate the inner sums, so that we do not have to calculate
2051      * them again for every atom */
2052     flex2_precalc_inner_sum(erg);
2053
2054     bCalcPotFit = (bOutstepRot || bOutstepSlab) && (RotationGroupFitting::Pot == erg->rotg->eFittype);
2055
2056     /********************************************************/
2057     /* Main loop over all local atoms of the rotation group */
2058     /********************************************************/
2059     N_M                                      = erg->rotg->nat * erg->invmass;
2060     V                                        = 0.0;
2061     OOsigma2                                 = 1.0 / (sigma * sigma);
2062     const auto& localRotationGroupIndex      = erg->atomSet->localIndex();
2063     const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
2064
2065     for (gmx::index j = 0; j < localRotationGroupIndex.ssize(); j++)
2066     {
2067         /* Local index of a rotation group atom  */
2068         ii = localRotationGroupIndex[j];
2069         /* Position of this atom in the collective array */
2070         iigrp = collectiveRotationGroupIndex[j];
2071         /* Mass-weighting */
2072         mj = erg->mc[iigrp]; /* need the unsorted mass here */
2073         wj = N_M * mj;
2074
2075         /* Current position of this atom: x[ii][XX/YY/ZZ]
2076          * Note that erg->xc_center contains the center of mass in case the flex2-t
2077          * potential was chosen. For the flex2 potential erg->xc_center must be
2078          * zero. */
2079         rvec_sub(coords[ii], erg->xc_center, xj);
2080
2081         /* Shift this atom such that it is near its reference */
2082         shift_single_coord(box, xj, erg->xc_shifts[iigrp]);
2083
2084         /* Determine the slabs to loop over, i.e. the ones with contributions
2085          * larger than min_gaussian */
2086         count = get_single_atom_gaussians(xj, erg);
2087
2088         clear_rvec(sum1vec_part);
2089         clear_rvec(sum2vec_part);
2090         sum3 = 0.0;
2091         sum4 = 0.0;
2092         /* Loop over the relevant slabs for this atom */
2093         for (int ic = 0; ic < count; ic++)
2094         {
2095             int n = erg->gn_slabind[ic];
2096
2097             /* Get the precomputed Gaussian value of curr_slab for curr_x */
2098             gaussian_xj = erg->gn_atom[ic];
2099
2100             int slabIndex = n - erg->slab_first; /* slab index */
2101
2102             /* The (unrotated) reference position of this atom is copied to yj0: */
2103             copy_rvec(erg->rotg->x_ref[iigrp], yj0);
2104
2105             beta = calc_beta(xj, erg, n);
2106
2107             /* The current center of this slab is saved in xcn: */
2108             copy_rvec(erg->slab_center[slabIndex], xcn);
2109             /* ... and the reference center in ycn: */
2110             copy_rvec(erg->slab_center_ref[slabIndex + erg->slab_buffer], ycn);
2111
2112             rvec_sub(yj0, ycn, yj0_ycn); /* yj0_ycn = yj0 - ycn      */
2113
2114             /* Rotate: */
2115             mvmul(erg->rotmat, yj0_ycn, rjn); /* rjn = Omega.(yj0 - ycn)  */
2116
2117             /* Subtract the slab center from xj */
2118             rvec_sub(xj, xcn, tmpvec2); /* tmpvec2 = xj - xcn       */
2119
2120             /* In rare cases, when an atom position coincides with a slab center
2121              * (tmpvec2 == 0) we cannot compute the vector product for sjn.
2122              * However, since the atom is located directly on the pivot, this
2123              * slab's contribution to the force on that atom will be zero
2124              * anyway. Therefore, we directly move on to the next slab.       */
2125             if (gmx_numzero(norm(tmpvec2))) /* 0 == norm(xj - xcn) */
2126             {
2127                 continue;
2128             }
2129
2130             /* Calculate sjn */
2131             cprod(erg->vec, tmpvec2, tmpvec); /* tmpvec = v x (xj - xcn)  */
2132
2133             OOpsijstar = norm2(tmpvec) + erg->rotg->eps; /* OOpsij* = 1/psij* = |v x (xj-xcn)|^2 + eps */
2134
2135             numerator = gmx::square(iprod(tmpvec, rjn));
2136
2137             /*********************************/
2138             /* Add to the rotation potential */
2139             /*********************************/
2140             V += 0.5 * erg->rotg->k * wj * gaussian_xj * numerator / OOpsijstar;
2141
2142             /* If requested, also calculate the potential for a set of angles
2143              * near the current reference angle */
2144             if (bCalcPotFit)
2145             {
2146                 for (int ifit = 0; ifit < erg->rotg->PotAngle_nstep; ifit++)
2147                 {
2148                     mvmul(erg->PotAngleFit->rotmat[ifit], yj0_ycn, fit_rjn);
2149                     fit_numerator = gmx::square(iprod(tmpvec, fit_rjn));
2150                     erg->PotAngleFit->V[ifit] +=
2151                             0.5 * erg->rotg->k * wj * gaussian_xj * fit_numerator / OOpsijstar;
2152                 }
2153             }
2154
2155             /*************************************/
2156             /* Now calculate the force on atom j */
2157             /*************************************/
2158
2159             OOpsij = norm(tmpvec); /* OOpsij = 1 / psij = |v x (xj - xcn)| */
2160
2161             /*                              *         v x (xj - xcn)          */
2162             unitv(tmpvec, sjn); /*  sjn = ----------------         */
2163                                 /*        |v x (xj - xcn)|         */
2164
2165             sjn_rjn = iprod(sjn, rjn); /* sjn_rjn = sjn . rjn             */
2166
2167
2168             /*** A. Calculate the first of the four sum terms: ****************/
2169             fac = OOpsij / OOpsijstar;
2170             svmul(fac, rjn, tmpvec);
2171             fac2 = fac * fac * OOpsij;
2172             svmul(fac2 * sjn_rjn, sjn, tmpvec2);
2173             rvec_dec(tmpvec, tmpvec2);
2174             fac2 = wj * gaussian_xj; /* also needed for sum4 */
2175             svmul(fac2 * sjn_rjn, tmpvec, slab_sum1vec_part);
2176             /********************/
2177             /*** Add to sum1: ***/
2178             /********************/
2179             rvec_inc(sum1vec_part, slab_sum1vec_part); /* sum1 still needs to vector multiplied with v */
2180
2181             /*** B. Calculate the forth of the four sum terms: ****************/
2182             betasigpsi = beta * OOsigma2 * OOpsij; /* this is also needed for sum3 */
2183             /********************/
2184             /*** Add to sum4: ***/
2185             /********************/
2186             slab_sum4part = fac2 * betasigpsi * fac * sjn_rjn
2187                             * sjn_rjn; /* Note that fac is still valid from above */
2188             sum4 += slab_sum4part;
2189
2190             /*** C. Calculate Wjn for second and third sum */
2191             /* Note that we can safely divide by slab_weights since we check in
2192              * get_slab_centers that it is non-zero. */
2193             Wjn = gaussian_xj * mj / erg->slab_weights[slabIndex];
2194
2195             /* We already have precalculated the inner sum for slab n */
2196             copy_rvec(erg->slab_innersumvec[slabIndex], innersumvec);
2197
2198             /* Weigh the inner sum vector with Wjn */
2199             svmul(Wjn, innersumvec, innersumvec);
2200
2201             /*** E. Calculate the second of the four sum terms: */
2202             /********************/
2203             /*** Add to sum2: ***/
2204             /********************/
2205             rvec_inc(sum2vec_part, innersumvec); /* sum2 still needs to be vector crossproduct'ed with v */
2206
2207             /*** F. Calculate the third of the four sum terms: */
2208             slab_sum3part = betasigpsi * iprod(sjn, innersumvec);
2209             sum3 += slab_sum3part; /* still needs to be multiplied with v */
2210
2211             /*** G. Calculate the torque on the local slab's axis: */
2212             if (bOutstepRot)
2213             {
2214                 /* Sum1 */
2215                 cprod(slab_sum1vec_part, erg->vec, slab_sum1vec);
2216                 /* Sum2 */
2217                 cprod(innersumvec, erg->vec, slab_sum2vec);
2218                 /* Sum3 */
2219                 svmul(slab_sum3part, erg->vec, slab_sum3vec);
2220                 /* Sum4 */
2221                 svmul(slab_sum4part, erg->vec, slab_sum4vec);
2222
2223                 /* The force on atom ii from slab n only: */
2224                 for (int m = 0; m < DIM; m++)
2225                 {
2226                     slab_force[m] = erg->rotg->k
2227                                     * (-slab_sum1vec[m] + slab_sum2vec[m] - slab_sum3vec[m]
2228                                        + 0.5 * slab_sum4vec[m]);
2229                 }
2230
2231                 erg->slab_torque_v[slabIndex] += torque(erg->vec, slab_force, xj, xcn);
2232             }
2233         } /* END of loop over slabs */
2234
2235         /* Construct the four individual parts of the vector sum: */
2236         cprod(sum1vec_part, erg->vec, sum1vec); /* sum1vec =   { } x v  */
2237         cprod(sum2vec_part, erg->vec, sum2vec); /* sum2vec =   { } x v  */
2238         svmul(sum3, erg->vec, sum3vec);         /* sum3vec =   { } . v  */
2239         svmul(sum4, erg->vec, sum4vec);         /* sum4vec =   { } . v  */
2240
2241         /* Store the additional force so that it can be added to the force
2242          * array after the normal forces have been evaluated */
2243         for (int m = 0; m < DIM; m++)
2244         {
2245             erg->f_rot_loc[j][m] =
2246                     erg->rotg->k * (-sum1vec[m] + sum2vec[m] - sum3vec[m] + 0.5 * sum4vec[m]);
2247         }
2248
2249 #ifdef SUM_PARTS
2250         fprintf(stderr,
2251                 "sum1: %15.8f %15.8f %15.8f\n",
2252                 -erg->rotg->k * sum1vec[XX],
2253                 -erg->rotg->k * sum1vec[YY],
2254                 -erg->rotg->k * sum1vec[ZZ]);
2255         fprintf(stderr,
2256                 "sum2: %15.8f %15.8f %15.8f\n",
2257                 erg->rotg->k * sum2vec[XX],
2258                 erg->rotg->k * sum2vec[YY],
2259                 erg->rotg->k * sum2vec[ZZ]);
2260         fprintf(stderr,
2261                 "sum3: %15.8f %15.8f %15.8f\n",
2262                 -erg->rotg->k * sum3vec[XX],
2263                 -erg->rotg->k * sum3vec[YY],
2264                 -erg->rotg->k * sum3vec[ZZ]);
2265         fprintf(stderr,
2266                 "sum4: %15.8f %15.8f %15.8f\n",
2267                 0.5 * erg->rotg->k * sum4vec[XX],
2268                 0.5 * erg->rotg->k * sum4vec[YY],
2269                 0.5 * erg->rotg->k * sum4vec[ZZ]);
2270 #endif
2271
2272         PRINT_FORCE_J
2273
2274     } /* END of loop over local atoms */
2275
2276     return V;
2277 }
2278
2279
2280 static real do_flex_lowlevel(gmx_enfrotgrp* erg,
2281                              real sigma, /* The Gaussian width sigma                      */
2282                              gmx::ArrayRef<const gmx::RVec> coords,
2283                              gmx_bool                       bOutstepRot,
2284                              gmx_bool                       bOutstepSlab,
2285                              const matrix                   box)
2286 {
2287     int      count, iigrp;
2288     rvec     xj, yj0;        /* current and reference position                */
2289     rvec     xcn, ycn;       /* the current and the reference slab centers    */
2290     rvec     yj0_ycn;        /* yj0 - ycn                                     */
2291     rvec     xj_xcn;         /* xj - xcn                                      */
2292     rvec     qjn, fit_qjn;   /* q_i^n                                         */
2293     rvec     sum_n1, sum_n2; /* Two contributions to the rotation force       */
2294     rvec     innersumvec;    /* Inner part of sum_n2                          */
2295     rvec     s_n;
2296     rvec     force_n;                /* Single force from slab n on one atom          */
2297     rvec     force_n1, force_n2;     /* First and second part of force_n              */
2298     rvec     tmpvec, tmpvec2, tmp_f; /* Helper variables                              */
2299     real     V;                      /* The rotation potential energy                 */
2300     real     OOsigma2;               /* 1/(sigma^2)                                   */
2301     real     beta;                   /* beta_n(xj)                                    */
2302     real     bjn, fit_bjn;           /* b_j^n                                         */
2303     real     gaussian_xj;            /* Gaussian weight gn(xj)                        */
2304     real     betan_xj_sigma2;
2305     real     mj, wj; /* Mass-weighting of the positions               */
2306     real     N_M;    /* N/M                                           */
2307     gmx_bool bCalcPotFit;
2308
2309     /* Pre-calculate the inner sums, so that we do not have to calculate
2310      * them again for every atom */
2311     flex_precalc_inner_sum(erg);
2312
2313     bCalcPotFit = (bOutstepRot || bOutstepSlab) && (RotationGroupFitting::Pot == erg->rotg->eFittype);
2314
2315     /********************************************************/
2316     /* Main loop over all local atoms of the rotation group */
2317     /********************************************************/
2318     OOsigma2                                 = 1.0 / (sigma * sigma);
2319     N_M                                      = erg->rotg->nat * erg->invmass;
2320     V                                        = 0.0;
2321     const auto& localRotationGroupIndex      = erg->atomSet->localIndex();
2322     const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
2323
2324     for (gmx::index j = 0; j < localRotationGroupIndex.ssize(); j++)
2325     {
2326         /* Local index of a rotation group atom  */
2327         int ii = localRotationGroupIndex[j];
2328         /* Position of this atom in the collective array */
2329         iigrp = collectiveRotationGroupIndex[j];
2330         /* Mass-weighting */
2331         mj = erg->mc[iigrp]; /* need the unsorted mass here */
2332         wj = N_M * mj;
2333
2334         /* Current position of this atom: x[ii][XX/YY/ZZ]
2335          * Note that erg->xc_center contains the center of mass in case the flex-t
2336          * potential was chosen. For the flex potential erg->xc_center must be
2337          * zero. */
2338         rvec_sub(coords[ii], erg->xc_center, xj);
2339
2340         /* Shift this atom such that it is near its reference */
2341         shift_single_coord(box, xj, erg->xc_shifts[iigrp]);
2342
2343         /* Determine the slabs to loop over, i.e. the ones with contributions
2344          * larger than min_gaussian */
2345         count = get_single_atom_gaussians(xj, erg);
2346
2347         clear_rvec(sum_n1);
2348         clear_rvec(sum_n2);
2349
2350         /* Loop over the relevant slabs for this atom */
2351         for (int ic = 0; ic < count; ic++)
2352         {
2353             int n = erg->gn_slabind[ic];
2354
2355             /* Get the precomputed Gaussian for xj in slab n */
2356             gaussian_xj = erg->gn_atom[ic];
2357
2358             int slabIndex = n - erg->slab_first; /* slab index */
2359
2360             /* The (unrotated) reference position of this atom is saved in yj0: */
2361             copy_rvec(erg->rotg->x_ref[iigrp], yj0);
2362
2363             beta = calc_beta(xj, erg, n);
2364
2365             /* The current center of this slab is saved in xcn: */
2366             copy_rvec(erg->slab_center[slabIndex], xcn);
2367             /* ... and the reference center in ycn: */
2368             copy_rvec(erg->slab_center_ref[slabIndex + erg->slab_buffer], ycn);
2369
2370             rvec_sub(yj0, ycn, yj0_ycn); /* yj0_ycn = yj0 - ycn */
2371
2372             /* In rare cases, when an atom position coincides with a reference slab
2373              * center (yj0_ycn == 0) we cannot compute the normal vector qjn.
2374              * However, since the atom is located directly on the pivot, this
2375              * slab's contribution to the force on that atom will be zero
2376              * anyway. Therefore, we directly move on to the next slab.       */
2377             if (gmx_numzero(norm(yj0_ycn))) /* 0 == norm(yj0 - ycn) */
2378             {
2379                 continue;
2380             }
2381
2382             /* Rotate: */
2383             mvmul(erg->rotmat, yj0_ycn, tmpvec2); /* tmpvec2= Omega.(yj0-ycn) */
2384
2385             /* Subtract the slab center from xj */
2386             rvec_sub(xj, xcn, xj_xcn); /* xj_xcn = xj - xcn         */
2387
2388             /* Calculate qjn */
2389             cprod(erg->vec, tmpvec2, tmpvec); /* tmpvec= v x Omega.(yj0-ycn) */
2390
2391             /*                         *         v x Omega.(yj0-ycn)    */
2392             unitv(tmpvec, qjn); /*  qjn = ---------------------   */
2393                                 /*        |v x Omega.(yj0-ycn)|   */
2394
2395             bjn = iprod(qjn, xj_xcn); /* bjn = qjn * (xj - xcn) */
2396
2397             /*********************************/
2398             /* Add to the rotation potential */
2399             /*********************************/
2400             V += 0.5 * erg->rotg->k * wj * gaussian_xj * gmx::square(bjn);
2401
2402             /* If requested, also calculate the potential for a set of angles
2403              * near the current reference angle */
2404             if (bCalcPotFit)
2405             {
2406                 for (int ifit = 0; ifit < erg->rotg->PotAngle_nstep; ifit++)
2407                 {
2408                     /* As above calculate Omega.(yj0-ycn), now for the other angles */
2409                     mvmul(erg->PotAngleFit->rotmat[ifit], yj0_ycn, tmpvec2); /* tmpvec2= Omega.(yj0-ycn) */
2410                     /* As above calculate qjn */
2411                     cprod(erg->vec, tmpvec2, tmpvec); /* tmpvec= v x Omega.(yj0-ycn) */
2412                     /*                                                        *             v x Omega.(yj0-ycn) */
2413                     unitv(tmpvec, fit_qjn);           /*  fit_qjn = ---------------------   */
2414                                                       /*            |v x Omega.(yj0-ycn)|   */
2415                     fit_bjn = iprod(fit_qjn, xj_xcn); /* fit_bjn = fit_qjn * (xj - xcn) */
2416                     /* Add to the rotation potential for this angle */
2417                     erg->PotAngleFit->V[ifit] +=
2418                             0.5 * erg->rotg->k * wj * gaussian_xj * gmx::square(fit_bjn);
2419                 }
2420             }
2421
2422             /****************************************************************/
2423             /* sum_n1 will typically be the main contribution to the force: */
2424             /****************************************************************/
2425             betan_xj_sigma2 = beta * OOsigma2; /*  beta_n(xj)/sigma^2  */
2426
2427             /* The next lines calculate
2428              *  qjn - (bjn*beta(xj)/(2sigma^2))v  */
2429             svmul(bjn * 0.5 * betan_xj_sigma2, erg->vec, tmpvec2);
2430             rvec_sub(qjn, tmpvec2, tmpvec);
2431
2432             /* Multiply with gn(xj)*bjn: */
2433             svmul(gaussian_xj * bjn, tmpvec, tmpvec2);
2434
2435             /* Sum over n: */
2436             rvec_inc(sum_n1, tmpvec2);
2437
2438             /* We already have precalculated the Sn term for slab n */
2439             copy_rvec(erg->slab_innersumvec[slabIndex], s_n);
2440             /*                                                             *          beta_n(xj) */
2441             svmul(betan_xj_sigma2 * iprod(s_n, xj_xcn), erg->vec, tmpvec); /* tmpvec = ---------- s_n (xj-xcn) */
2442             /*            sigma^2               */
2443
2444             rvec_sub(s_n, tmpvec, innersumvec);
2445
2446             /* We can safely divide by slab_weights since we check in get_slab_centers
2447              * that it is non-zero. */
2448             svmul(gaussian_xj / erg->slab_weights[slabIndex], innersumvec, innersumvec);
2449
2450             rvec_add(sum_n2, innersumvec, sum_n2);
2451
2452             /* Calculate the torque: */
2453             if (bOutstepRot)
2454             {
2455                 /* The force on atom ii from slab n only: */
2456                 svmul(-erg->rotg->k * wj, tmpvec2, force_n1);    /* part 1 */
2457                 svmul(erg->rotg->k * mj, innersumvec, force_n2); /* part 2 */
2458                 rvec_add(force_n1, force_n2, force_n);
2459                 erg->slab_torque_v[slabIndex] += torque(erg->vec, force_n, xj, xcn);
2460             }
2461         } /* END of loop over slabs */
2462
2463         /* Put both contributions together: */
2464         svmul(wj, sum_n1, sum_n1);
2465         svmul(mj, sum_n2, sum_n2);
2466         rvec_sub(sum_n2, sum_n1, tmp_f); /* F = -grad V */
2467
2468         /* Store the additional force so that it can be added to the force
2469          * array after the normal forces have been evaluated */
2470         for (int m = 0; m < DIM; m++)
2471         {
2472             erg->f_rot_loc[j][m] = erg->rotg->k * tmp_f[m];
2473         }
2474
2475         PRINT_FORCE_J
2476
2477     } /* END of loop over local atoms */
2478
2479     return V;
2480 }
2481
2482 static void sort_collective_coordinates(gmx_enfrotgrp* erg,
2483                                         sort_along_vec_t* data) /* Buffer for sorting the positions */
2484 {
2485     /* The projection of the position vector on the rotation vector is
2486      * the relevant value for sorting. Fill the 'data' structure */
2487     for (int i = 0; i < erg->rotg->nat; i++)
2488     {
2489         data[i].xcproj = iprod(erg->xc[i], erg->vec); /* sort criterium */
2490         data[i].m      = erg->mc[i];
2491         data[i].ind    = i;
2492         copy_rvec(erg->xc[i], data[i].x);
2493         copy_rvec(erg->rotg->x_ref[i], data[i].x_ref);
2494     }
2495     /* Sort the 'data' structure */
2496     std::sort(data, data + erg->rotg->nat, [](const sort_along_vec_t& a, const sort_along_vec_t& b) {
2497         return a.xcproj < b.xcproj;
2498     });
2499
2500     /* Copy back the sorted values */
2501     for (int i = 0; i < erg->rotg->nat; i++)
2502     {
2503         copy_rvec(data[i].x, erg->xc[i]);
2504         copy_rvec(data[i].x_ref, erg->xc_ref_sorted[i]);
2505         erg->mc_sorted[i]  = data[i].m;
2506         erg->xc_sortind[i] = data[i].ind;
2507     }
2508 }
2509
2510
2511 /* For each slab, get the first and the last index of the sorted atom
2512  * indices */
2513 static void get_firstlast_atom_per_slab(const gmx_enfrotgrp* erg)
2514 {
2515     real beta;
2516
2517     /* Find the first atom that needs to enter the calculation for each slab */
2518     int n = erg->slab_first; /* slab */
2519     int i = 0;               /* start with the first atom */
2520     do
2521     {
2522         /* Find the first atom that significantly contributes to this slab */
2523         do /* move forward in position until a large enough beta is found */
2524         {
2525             beta = calc_beta(erg->xc[i], erg, n);
2526             i++;
2527         } while ((beta < -erg->max_beta) && (i < erg->rotg->nat));
2528         i--;
2529         int slabIndex             = n - erg->slab_first; /* slab index */
2530         erg->firstatom[slabIndex] = i;
2531         /* Proceed to the next slab */
2532         n++;
2533     } while (n <= erg->slab_last);
2534
2535     /* Find the last atom for each slab */
2536     n = erg->slab_last;     /* start with last slab */
2537     i = erg->rotg->nat - 1; /* start with the last atom */
2538     do
2539     {
2540         do /* move backward in position until a large enough beta is found */
2541         {
2542             beta = calc_beta(erg->xc[i], erg, n);
2543             i--;
2544         } while ((beta > erg->max_beta) && (i > -1));
2545         i++;
2546         int slabIndex            = n - erg->slab_first; /* slab index */
2547         erg->lastatom[slabIndex] = i;
2548         /* Proceed to the next slab */
2549         n--;
2550     } while (n >= erg->slab_first);
2551 }
2552
2553
2554 /* Determine the very first and very last slab that needs to be considered
2555  * For the first slab that needs to be considered, we have to find the smallest
2556  * n that obeys:
2557  *
2558  * x_first * v - n*Delta_x <= beta_max
2559  *
2560  * slab index n, slab distance Delta_x, rotation vector v. For the last slab we
2561  * have to find the largest n that obeys
2562  *
2563  * x_last * v - n*Delta_x >= -beta_max
2564  *
2565  */
2566 static inline int get_first_slab(const gmx_enfrotgrp* erg,
2567                                  rvec firstatom) /* First atom after sorting along the rotation vector v */
2568 {
2569     /* Find the first slab for the first atom */
2570     return static_cast<int>(ceil(
2571             static_cast<double>((iprod(firstatom, erg->vec) - erg->max_beta) / erg->rotg->slab_dist)));
2572 }
2573
2574
2575 static inline int get_last_slab(const gmx_enfrotgrp* erg, rvec lastatom) /* Last atom along v */
2576 {
2577     /* Find the last slab for the last atom */
2578     return static_cast<int>(floor(
2579             static_cast<double>((iprod(lastatom, erg->vec) + erg->max_beta) / erg->rotg->slab_dist)));
2580 }
2581
2582
2583 static void get_firstlast_slab_check(gmx_enfrotgrp* erg, /* The rotation group (data only accessible in this file) */
2584                                      rvec firstatom, /* First atom after sorting along the rotation vector v */
2585                                      rvec lastatom) /* Last atom along v */
2586 {
2587     erg->slab_first = get_first_slab(erg, firstatom);
2588     erg->slab_last  = get_last_slab(erg, lastatom);
2589
2590     /* Calculate the slab buffer size, which changes when slab_first changes */
2591     erg->slab_buffer = erg->slab_first - erg->slab_first_ref;
2592
2593     /* Check whether we have reference data to compare against */
2594     if (erg->slab_first < erg->slab_first_ref)
2595     {
2596         gmx_fatal(FARGS,
2597                   "%s No reference data for first slab (n=%d), unable to proceed.",
2598                   RotStr.c_str(),
2599                   erg->slab_first);
2600     }
2601
2602     /* Check whether we have reference data to compare against */
2603     if (erg->slab_last > erg->slab_last_ref)
2604     {
2605         gmx_fatal(FARGS,
2606                   "%s No reference data for last slab (n=%d), unable to proceed.",
2607                   RotStr.c_str(),
2608                   erg->slab_last);
2609     }
2610 }
2611
2612
2613 /* Enforced rotation with a flexible axis */
2614 static void do_flexible(gmx_bool       bMaster,
2615                         gmx_enfrot*    enfrot, /* Other rotation data                        */
2616                         gmx_enfrotgrp* erg,
2617                         gmx::ArrayRef<const gmx::RVec> coords, /* The local positions */
2618                         const matrix                   box,
2619                         double   t,            /* Time in picoseconds                        */
2620                         gmx_bool bOutstepRot,  /* Output to main rotation output file        */
2621                         gmx_bool bOutstepSlab) /* Output per-slab data                       */
2622 {
2623     int  nslabs;
2624     real sigma; /* The Gaussian width sigma */
2625
2626     /* Define the sigma value */
2627     sigma = 0.7 * erg->rotg->slab_dist;
2628
2629     /* Sort the collective coordinates erg->xc along the rotation vector. This is
2630      * an optimization for the inner loop. */
2631     sort_collective_coordinates(erg, enfrot->data);
2632
2633     /* Determine the first relevant slab for the first atom and the last
2634      * relevant slab for the last atom */
2635     get_firstlast_slab_check(erg, erg->xc[0], erg->xc[erg->rotg->nat - 1]);
2636
2637     /* Determine for each slab depending on the min_gaussian cutoff criterium,
2638      * a first and a last atom index inbetween stuff needs to be calculated */
2639     get_firstlast_atom_per_slab(erg);
2640
2641     /* Determine the gaussian-weighted center of positions for all slabs */
2642     get_slab_centers(erg, erg->xc, erg->mc_sorted, t, enfrot->out_slabs, bOutstepSlab, FALSE);
2643
2644     /* Clear the torque per slab from last time step: */
2645     nslabs = erg->slab_last - erg->slab_first + 1;
2646     for (int l = 0; l < nslabs; l++)
2647     {
2648         erg->slab_torque_v[l] = 0.0;
2649     }
2650
2651     /* Call the rotational forces kernel */
2652     if (erg->rotg->eType == EnforcedRotationGroupType::Flex
2653         || erg->rotg->eType == EnforcedRotationGroupType::Flext)
2654     {
2655         erg->V = do_flex_lowlevel(erg, sigma, coords, bOutstepRot, bOutstepSlab, box);
2656     }
2657     else if (erg->rotg->eType == EnforcedRotationGroupType::Flex2
2658              || erg->rotg->eType == EnforcedRotationGroupType::Flex2t)
2659     {
2660         erg->V = do_flex2_lowlevel(erg, sigma, coords, bOutstepRot, bOutstepSlab, box);
2661     }
2662     else
2663     {
2664         gmx_fatal(FARGS, "Unknown flexible rotation type");
2665     }
2666
2667     /* Determine angle by RMSD fit to the reference - Let's hope this */
2668     /* only happens once in a while, since this is not parallelized! */
2669     if (bMaster && (RotationGroupFitting::Pot != erg->rotg->eFittype))
2670     {
2671         if (bOutstepRot)
2672         {
2673             /* Fit angle of the whole rotation group */
2674             erg->angle_v = flex_fit_angle(erg);
2675         }
2676         if (bOutstepSlab)
2677         {
2678             /* Fit angle of each slab */
2679             flex_fit_angle_perslab(erg, t, erg->degangle, enfrot->out_angles);
2680         }
2681     }
2682
2683     /* Lump together the torques from all slabs: */
2684     erg->torque_v = 0.0;
2685     for (int l = 0; l < nslabs; l++)
2686     {
2687         erg->torque_v += erg->slab_torque_v[l];
2688     }
2689 }
2690
2691
2692 /* Calculate the angle between reference and actual rotation group atom,
2693  * both projected into a plane perpendicular to the rotation vector: */
2694 static void angle(const gmx_enfrotgrp* erg,
2695                   rvec                 x_act,
2696                   rvec                 x_ref,
2697                   real*                alpha,
2698                   real* weight) /* atoms near the rotation axis should count less than atoms far away */
2699 {
2700     rvec xp, xrp; /* current and reference positions projected on a plane perpendicular to pg->vec */
2701     rvec dum;
2702
2703
2704     /* Project x_ref and x into a plane through the origin perpendicular to rot_vec: */
2705     /* Project x_ref: xrp = x_ref - (vec * x_ref) * vec */
2706     svmul(iprod(erg->vec, x_ref), erg->vec, dum);
2707     rvec_sub(x_ref, dum, xrp);
2708     /* Project x_act: */
2709     svmul(iprod(erg->vec, x_act), erg->vec, dum);
2710     rvec_sub(x_act, dum, xp);
2711
2712     /* Retrieve information about which vector precedes. gmx_angle always
2713      * returns a positive angle. */
2714     cprod(xp, xrp, dum); /* if reference precedes, this is pointing into the same direction as vec */
2715
2716     if (iprod(erg->vec, dum) >= 0)
2717     {
2718         *alpha = -gmx_angle(xrp, xp);
2719     }
2720     else
2721     {
2722         *alpha = +gmx_angle(xrp, xp);
2723     }
2724
2725     /* Also return the weight */
2726     *weight = norm(xp);
2727 }
2728
2729
2730 /* Project first vector onto a plane perpendicular to the second vector
2731  * dr = dr - (dr.v)v
2732  * Note that v must be of unit length.
2733  */
2734 static inline void project_onto_plane(rvec dr, const rvec v)
2735 {
2736     rvec tmp;
2737
2738
2739     svmul(iprod(dr, v), v, tmp); /* tmp = (dr.v)v */
2740     rvec_dec(dr, tmp);           /* dr = dr - (dr.v)v */
2741 }
2742
2743
2744 /* Fixed rotation: The rotation reference group rotates around the v axis. */
2745 /* The atoms of the actual rotation group are attached with imaginary  */
2746 /* springs to the reference atoms.                                     */
2747 static void do_fixed(gmx_enfrotgrp* erg,
2748                      gmx_bool       bOutstepRot, /* Output to main rotation output file        */
2749                      gmx_bool       bOutstepSlab)      /* Output per-slab data                       */
2750 {
2751     rvec     dr;
2752     rvec     tmp_f;  /* Force */
2753     real     alpha;  /* a single angle between an actual and a reference position */
2754     real     weight; /* single weight for a single angle */
2755     rvec     xi_xc;  /* xi - xc */
2756     gmx_bool bCalcPotFit;
2757     rvec     fit_xr_loc;
2758
2759     /* for mass weighting: */
2760     real wi;   /* Mass-weighting of the positions */
2761     real N_M;  /* N/M */
2762     real k_wi; /* k times wi */
2763
2764     gmx_bool bProject;
2765
2766     bProject = (erg->rotg->eType == EnforcedRotationGroupType::Pm)
2767                || (erg->rotg->eType == EnforcedRotationGroupType::Pmpf);
2768     bCalcPotFit = (bOutstepRot || bOutstepSlab) && (RotationGroupFitting::Pot == erg->rotg->eFittype);
2769
2770     N_M                                      = erg->rotg->nat * erg->invmass;
2771     const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
2772     /* Each process calculates the forces on its local atoms */
2773     for (size_t j = 0; j < erg->atomSet->numAtomsLocal(); j++)
2774     {
2775         /* Calculate (x_i-x_c) resp. (x_i-u) */
2776         rvec_sub(erg->x_loc_pbc[j], erg->xc_center, xi_xc);
2777
2778         /* Calculate Omega*(y_i-y_c)-(x_i-x_c) */
2779         rvec_sub(erg->xr_loc[j], xi_xc, dr);
2780
2781         if (bProject)
2782         {
2783             project_onto_plane(dr, erg->vec);
2784         }
2785
2786         /* Mass-weighting */
2787         wi = N_M * erg->m_loc[j];
2788
2789         /* Store the additional force so that it can be added to the force
2790          * array after the normal forces have been evaluated */
2791         k_wi = erg->rotg->k * wi;
2792         for (int m = 0; m < DIM; m++)
2793         {
2794             tmp_f[m]             = k_wi * dr[m];
2795             erg->f_rot_loc[j][m] = tmp_f[m];
2796             erg->V += 0.5 * k_wi * gmx::square(dr[m]);
2797         }
2798
2799         /* If requested, also calculate the potential for a set of angles
2800          * near the current reference angle */
2801         if (bCalcPotFit)
2802         {
2803             for (int ifit = 0; ifit < erg->rotg->PotAngle_nstep; ifit++)
2804             {
2805                 /* Index of this rotation group atom with respect to the whole rotation group */
2806                 int jj = collectiveRotationGroupIndex[j];
2807
2808                 /* Rotate with the alternative angle. Like rotate_local_reference(),
2809                  * just for a single local atom */
2810                 mvmul(erg->PotAngleFit->rotmat[ifit], erg->rotg->x_ref[jj], fit_xr_loc); /* fit_xr_loc = Omega*(y_i-y_c) */
2811
2812                 /* Calculate Omega*(y_i-y_c)-(x_i-x_c) */
2813                 rvec_sub(fit_xr_loc, xi_xc, dr);
2814
2815                 if (bProject)
2816                 {
2817                     project_onto_plane(dr, erg->vec);
2818                 }
2819
2820                 /* Add to the rotation potential for this angle: */
2821                 erg->PotAngleFit->V[ifit] += 0.5 * k_wi * norm2(dr);
2822             }
2823         }
2824
2825         if (bOutstepRot)
2826         {
2827             /* Add to the torque of this rotation group */
2828             erg->torque_v += torque(erg->vec, tmp_f, erg->x_loc_pbc[j], erg->xc_center);
2829
2830             /* Calculate the angle between reference and actual rotation group atom. */
2831             angle(erg, xi_xc, erg->xr_loc[j], &alpha, &weight); /* angle in rad, weighted */
2832             erg->angle_v += alpha * weight;
2833             erg->weight_v += weight;
2834         }
2835         /* If you want enforced rotation to contribute to the virial,
2836          * activate the following lines:
2837             if (MASTER(cr))
2838             {
2839                Add the rotation contribution to the virial
2840               for(j=0; j<DIM; j++)
2841                 for(m=0;m<DIM;m++)
2842                   vir[j][m] += 0.5*f[ii][j]*dr[m];
2843             }
2844          */
2845
2846         PRINT_FORCE_J
2847
2848     } /* end of loop over local rotation group atoms */
2849 }
2850
2851
2852 /* Calculate the radial motion potential and forces */
2853 static void do_radial_motion(gmx_enfrotgrp* erg,
2854                              gmx_bool bOutstepRot,  /* Output to main rotation output file        */
2855                              gmx_bool bOutstepSlab) /* Output per-slab data                       */
2856 {
2857     rvec     tmp_f;  /* Force */
2858     real     alpha;  /* a single angle between an actual and a reference position */
2859     real     weight; /* single weight for a single angle */
2860     rvec     xj_u;   /* xj - u */
2861     rvec     tmpvec, fit_tmpvec;
2862     real     fac, fac2, sum = 0.0;
2863     rvec     pj;
2864     gmx_bool bCalcPotFit;
2865
2866     /* For mass weighting: */
2867     real wj;  /* Mass-weighting of the positions */
2868     real N_M; /* N/M */
2869
2870     bCalcPotFit = (bOutstepRot || bOutstepSlab) && (RotationGroupFitting::Pot == erg->rotg->eFittype);
2871
2872     N_M                                      = erg->rotg->nat * erg->invmass;
2873     const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
2874     /* Each process calculates the forces on its local atoms */
2875     for (size_t j = 0; j < erg->atomSet->numAtomsLocal(); j++)
2876     {
2877         /* Calculate (xj-u) */
2878         rvec_sub(erg->x_loc_pbc[j], erg->xc_center, xj_u); /* xj_u = xj-u */
2879
2880         /* Calculate Omega.(yj0-u) */
2881         cprod(erg->vec, erg->xr_loc[j], tmpvec); /* tmpvec = v x Omega.(yj0-u) */
2882
2883         /*                       *         v x Omega.(yj0-u)     */
2884         unitv(tmpvec, pj); /*  pj = ---------------------   */
2885                            /*       | v x Omega.(yj0-u) |   */
2886
2887         fac  = iprod(pj, xj_u); /* fac = pj.(xj-u) */
2888         fac2 = fac * fac;
2889
2890         /* Mass-weighting */
2891         wj = N_M * erg->m_loc[j];
2892
2893         /* Store the additional force so that it can be added to the force
2894          * array after the normal forces have been evaluated */
2895         svmul(-erg->rotg->k * wj * fac, pj, tmp_f);
2896         copy_rvec(tmp_f, erg->f_rot_loc[j]);
2897         sum += wj * fac2;
2898
2899         /* If requested, also calculate the potential for a set of angles
2900          * near the current reference angle */
2901         if (bCalcPotFit)
2902         {
2903             for (int ifit = 0; ifit < erg->rotg->PotAngle_nstep; ifit++)
2904             {
2905                 /* Index of this rotation group atom with respect to the whole rotation group */
2906                 int jj = collectiveRotationGroupIndex[j];
2907
2908                 /* Rotate with the alternative angle. Like rotate_local_reference(),
2909                  * just for a single local atom */
2910                 mvmul(erg->PotAngleFit->rotmat[ifit], erg->rotg->x_ref[jj], fit_tmpvec); /* fit_tmpvec = Omega*(yj0-u) */
2911
2912                 /* Calculate Omega.(yj0-u) */
2913                 cprod(erg->vec, fit_tmpvec, tmpvec); /* tmpvec = v x Omega.(yj0-u) */
2914                 /*                                     *         v x Omega.(yj0-u)     */
2915                 unitv(tmpvec, pj); /*  pj = ---------------------   */
2916                                    /*       | v x Omega.(yj0-u) |   */
2917
2918                 fac  = iprod(pj, xj_u); /* fac = pj.(xj-u) */
2919                 fac2 = fac * fac;
2920
2921                 /* Add to the rotation potential for this angle: */
2922                 erg->PotAngleFit->V[ifit] += 0.5 * erg->rotg->k * wj * fac2;
2923             }
2924         }
2925
2926         if (bOutstepRot)
2927         {
2928             /* Add to the torque of this rotation group */
2929             erg->torque_v += torque(erg->vec, tmp_f, erg->x_loc_pbc[j], erg->xc_center);
2930
2931             /* Calculate the angle between reference and actual rotation group atom. */
2932             angle(erg, xj_u, erg->xr_loc[j], &alpha, &weight); /* angle in rad, weighted */
2933             erg->angle_v += alpha * weight;
2934             erg->weight_v += weight;
2935         }
2936
2937         PRINT_FORCE_J
2938
2939     } /* end of loop over local rotation group atoms */
2940     erg->V = 0.5 * erg->rotg->k * sum;
2941 }
2942
2943
2944 /* Calculate the radial motion pivot-free potential and forces */
2945 static void do_radial_motion_pf(gmx_enfrotgrp*                 erg,
2946                                 gmx::ArrayRef<const gmx::RVec> coords, /* The positions */
2947                                 const matrix box, /* The simulation box                         */
2948                                 gmx_bool     bOutstepRot, /* Output to main rotation output file  */
2949                                 gmx_bool     bOutstepSlab)    /* Output per-slab data */
2950 {
2951     rvec     xj;      /* Current position */
2952     rvec     xj_xc;   /* xj  - xc  */
2953     rvec     yj0_yc0; /* yj0 - yc0 */
2954     rvec     tmp_f;   /* Force */
2955     real     alpha;   /* a single angle between an actual and a reference position */
2956     real     weight;  /* single weight for a single angle */
2957     rvec     tmpvec, tmpvec2;
2958     rvec     innersumvec; /* Precalculation of the inner sum */
2959     rvec     innersumveckM;
2960     real     fac, fac2, V = 0.0;
2961     rvec     qi, qj;
2962     gmx_bool bCalcPotFit;
2963
2964     /* For mass weighting: */
2965     real mj, wi, wj; /* Mass-weighting of the positions */
2966     real N_M;        /* N/M */
2967
2968     bCalcPotFit = (bOutstepRot || bOutstepSlab) && (RotationGroupFitting::Pot == erg->rotg->eFittype);
2969
2970     N_M = erg->rotg->nat * erg->invmass;
2971
2972     /* Get the current center of the rotation group: */
2973     get_center(erg->xc, erg->mc, erg->rotg->nat, erg->xc_center);
2974
2975     /* Precalculate Sum_i [ wi qi.(xi-xc) qi ] which is needed for every single j */
2976     clear_rvec(innersumvec);
2977     for (int i = 0; i < erg->rotg->nat; i++)
2978     {
2979         /* Mass-weighting */
2980         wi = N_M * erg->mc[i];
2981
2982         /* Calculate qi. Note that xc_ref_center has already been subtracted from
2983          * x_ref in init_rot_group.*/
2984         mvmul(erg->rotmat, erg->rotg->x_ref[i], tmpvec); /* tmpvec  = Omega.(yi0-yc0) */
2985
2986         cprod(erg->vec, tmpvec, tmpvec2); /* tmpvec2 = v x Omega.(yi0-yc0) */
2987
2988         /*                                             *         v x Omega.(yi0-yc0)     */
2989         unitv(tmpvec2, qi); /*  qi = -----------------------   */
2990                             /*       | v x Omega.(yi0-yc0) |   */
2991
2992         rvec_sub(erg->xc[i], erg->xc_center, tmpvec); /* tmpvec = xi-xc */
2993
2994         svmul(wi * iprod(qi, tmpvec), qi, tmpvec2);
2995
2996         rvec_inc(innersumvec, tmpvec2);
2997     }
2998     svmul(erg->rotg->k * erg->invmass, innersumvec, innersumveckM);
2999
3000     /* Each process calculates the forces on its local atoms */
3001     const auto& localRotationGroupIndex      = erg->atomSet->localIndex();
3002     const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
3003     for (gmx::index j = 0; j < localRotationGroupIndex.ssize(); j++)
3004     {
3005         /* Local index of a rotation group atom  */
3006         int ii = localRotationGroupIndex[j];
3007         /* Position of this atom in the collective array */
3008         int iigrp = collectiveRotationGroupIndex[j];
3009         /* Mass-weighting */
3010         mj = erg->mc[iigrp]; /* need the unsorted mass here */
3011         wj = N_M * mj;
3012
3013         /* Current position of this atom: x[ii][XX/YY/ZZ] */
3014         copy_rvec(coords[ii], xj);
3015
3016         /* Shift this atom such that it is near its reference */
3017         shift_single_coord(box, xj, erg->xc_shifts[iigrp]);
3018
3019         /* The (unrotated) reference position is yj0. yc0 has already
3020          * been subtracted in init_rot_group */
3021         copy_rvec(erg->rotg->x_ref[iigrp], yj0_yc0); /* yj0_yc0 = yj0 - yc0      */
3022
3023         /* Calculate Omega.(yj0-yc0) */
3024         mvmul(erg->rotmat, yj0_yc0, tmpvec2); /* tmpvec2 = Omega.(yj0 - yc0)  */
3025
3026         cprod(erg->vec, tmpvec2, tmpvec); /* tmpvec = v x Omega.(yj0-yc0) */
3027
3028         /*                     *         v x Omega.(yj0-yc0)     */
3029         unitv(tmpvec, qj); /*  qj = -----------------------   */
3030                            /*       | v x Omega.(yj0-yc0) |   */
3031
3032         /* Calculate (xj-xc) */
3033         rvec_sub(xj, erg->xc_center, xj_xc); /* xj_xc = xj-xc */
3034
3035         fac  = iprod(qj, xj_xc); /* fac = qj.(xj-xc) */
3036         fac2 = fac * fac;
3037
3038         /* Store the additional force so that it can be added to the force
3039          * array after the normal forces have been evaluated */
3040         svmul(-erg->rotg->k * wj * fac, qj, tmp_f); /* part 1 of force */
3041         svmul(mj, innersumveckM, tmpvec);           /* part 2 of force */
3042         rvec_inc(tmp_f, tmpvec);
3043         copy_rvec(tmp_f, erg->f_rot_loc[j]);
3044         V += wj * fac2;
3045
3046         /* If requested, also calculate the potential for a set of angles
3047          * near the current reference angle */
3048         if (bCalcPotFit)
3049         {
3050             for (int ifit = 0; ifit < erg->rotg->PotAngle_nstep; ifit++)
3051             {
3052                 /* Rotate with the alternative angle. Like rotate_local_reference(),
3053                  * just for a single local atom */
3054                 mvmul(erg->PotAngleFit->rotmat[ifit], yj0_yc0, tmpvec2); /* tmpvec2 = Omega*(yj0-yc0) */
3055
3056                 /* Calculate Omega.(yj0-u) */
3057                 cprod(erg->vec, tmpvec2, tmpvec); /* tmpvec = v x Omega.(yj0-yc0) */
3058                 /*                                  *         v x Omega.(yj0-yc0)     */
3059                 unitv(tmpvec, qj); /*  qj = -----------------------   */
3060                                    /*       | v x Omega.(yj0-yc0) |   */
3061
3062                 fac  = iprod(qj, xj_xc); /* fac = qj.(xj-xc) */
3063                 fac2 = fac * fac;
3064
3065                 /* Add to the rotation potential for this angle: */
3066                 erg->PotAngleFit->V[ifit] += 0.5 * erg->rotg->k * wj * fac2;
3067             }
3068         }
3069
3070         if (bOutstepRot)
3071         {
3072             /* Add to the torque of this rotation group */
3073             erg->torque_v += torque(erg->vec, tmp_f, xj, erg->xc_center);
3074
3075             /* Calculate the angle between reference and actual rotation group atom. */
3076             angle(erg, xj_xc, yj0_yc0, &alpha, &weight); /* angle in rad, weighted */
3077             erg->angle_v += alpha * weight;
3078             erg->weight_v += weight;
3079         }
3080
3081         PRINT_FORCE_J
3082
3083     } /* end of loop over local rotation group atoms */
3084     erg->V = 0.5 * erg->rotg->k * V;
3085 }
3086
3087
3088 /* Precalculate the inner sum for the radial motion 2 forces */
3089 static void radial_motion2_precalc_inner_sum(const gmx_enfrotgrp* erg, rvec innersumvec)
3090 {
3091     int  i;
3092     rvec xi_xc; /* xj - xc */
3093     rvec tmpvec, tmpvec2;
3094     real fac;
3095     rvec ri, si;
3096     real siri;
3097     rvec v_xi_xc; /* v x (xj - u) */
3098     real psii, psiistar;
3099     real wi;  /* Mass-weighting of the positions */
3100     real N_M; /* N/M */
3101     rvec sumvec;
3102
3103     N_M = erg->rotg->nat * erg->invmass;
3104
3105     /* Loop over the collective set of positions */
3106     clear_rvec(sumvec);
3107     for (i = 0; i < erg->rotg->nat; i++)
3108     {
3109         /* Mass-weighting */
3110         wi = N_M * erg->mc[i];
3111
3112         rvec_sub(erg->xc[i], erg->xc_center, xi_xc); /* xi_xc = xi-xc         */
3113
3114         /* Calculate ri. Note that xc_ref_center has already been subtracted from
3115          * x_ref in init_rot_group.*/
3116         mvmul(erg->rotmat, erg->rotg->x_ref[i], ri); /* ri  = Omega.(yi0-yc0) */
3117
3118         cprod(erg->vec, xi_xc, v_xi_xc); /* v_xi_xc = v x (xi-u)  */
3119
3120         fac = norm2(v_xi_xc);
3121         /*                                 *                      1           */
3122         psiistar = 1.0 / (fac + erg->rotg->eps); /* psiistar = --------------------- */
3123         /*            |v x (xi-xc)|^2 + eps */
3124
3125         psii = gmx::invsqrt(fac); /*                 1                */
3126                                   /*  psii    = -------------         */
3127                                   /*            |v x (xi-xc)|         */
3128
3129         svmul(psii, v_xi_xc, si); /*  si = psii * (v x (xi-xc) )     */
3130
3131         siri = iprod(si, ri); /* siri = si.ri           */
3132
3133         svmul(psiistar / psii, ri, tmpvec);
3134         svmul(psiistar * psiistar / (psii * psii * psii) * siri, si, tmpvec2);
3135         rvec_dec(tmpvec, tmpvec2);
3136         cprod(tmpvec, erg->vec, tmpvec2);
3137
3138         svmul(wi * siri, tmpvec2, tmpvec);
3139
3140         rvec_inc(sumvec, tmpvec);
3141     }
3142     svmul(erg->rotg->k * erg->invmass, sumvec, innersumvec);
3143 }
3144
3145
3146 /* Calculate the radial motion 2 potential and forces */
3147 static void do_radial_motion2(gmx_enfrotgrp*                 erg,
3148                               gmx::ArrayRef<const gmx::RVec> coords, /* The positions */
3149                               const matrix box,     /* The simulation box                         */
3150                               gmx_bool bOutstepRot, /* Output to main rotation output file        */
3151                               gmx_bool bOutstepSlab) /* Output per-slab data */
3152 {
3153     rvec     xj;      /* Position */
3154     real     alpha;   /* a single angle between an actual and a reference position */
3155     real     weight;  /* single weight for a single angle */
3156     rvec     xj_u;    /* xj - u */
3157     rvec     yj0_yc0; /* yj0 -yc0 */
3158     rvec     tmpvec, tmpvec2;
3159     real     fac, fit_fac, fac2, Vpart = 0.0;
3160     rvec     rj, fit_rj, sj;
3161     real     sjrj;
3162     rvec     v_xj_u; /* v x (xj - u) */
3163     real     psij, psijstar;
3164     real     mj, wj; /* For mass-weighting of the positions */
3165     real     N_M;    /* N/M */
3166     gmx_bool bPF;
3167     rvec     innersumvec;
3168     gmx_bool bCalcPotFit;
3169
3170     bPF         = erg->rotg->eType == EnforcedRotationGroupType::Rm2pf;
3171     bCalcPotFit = (bOutstepRot || bOutstepSlab) && (RotationGroupFitting::Pot == erg->rotg->eFittype);
3172
3173     clear_rvec(yj0_yc0); /* Make the compiler happy */
3174
3175     clear_rvec(innersumvec);
3176     if (bPF)
3177     {
3178         /* For the pivot-free variant we have to use the current center of
3179          * mass of the rotation group instead of the pivot u */
3180         get_center(erg->xc, erg->mc, erg->rotg->nat, erg->xc_center);
3181
3182         /* Also, we precalculate the second term of the forces that is identical
3183          * (up to the weight factor mj) for all forces */
3184         radial_motion2_precalc_inner_sum(erg, innersumvec);
3185     }
3186
3187     N_M = erg->rotg->nat * erg->invmass;
3188
3189     /* Each process calculates the forces on its local atoms */
3190     const auto& localRotationGroupIndex      = erg->atomSet->localIndex();
3191     const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
3192     for (gmx::index j = 0; j < localRotationGroupIndex.ssize(); j++)
3193     {
3194         if (bPF)
3195         {
3196             /* Local index of a rotation group atom  */
3197             int ii = localRotationGroupIndex[j];
3198             /* Position of this atom in the collective array */
3199             int iigrp = collectiveRotationGroupIndex[j];
3200             /* Mass-weighting */
3201             mj = erg->mc[iigrp];
3202
3203             /* Current position of this atom: x[ii] */
3204             copy_rvec(coords[ii], xj);
3205
3206             /* Shift this atom such that it is near its reference */
3207             shift_single_coord(box, xj, erg->xc_shifts[iigrp]);
3208
3209             /* The (unrotated) reference position is yj0. yc0 has already
3210              * been subtracted in init_rot_group */
3211             copy_rvec(erg->rotg->x_ref[iigrp], yj0_yc0); /* yj0_yc0 = yj0 - yc0  */
3212
3213             /* Calculate Omega.(yj0-yc0) */
3214             mvmul(erg->rotmat, yj0_yc0, rj); /* rj = Omega.(yj0-yc0)  */
3215         }
3216         else
3217         {
3218             mj = erg->m_loc[j];
3219             copy_rvec(erg->x_loc_pbc[j], xj);
3220             copy_rvec(erg->xr_loc[j], rj); /* rj = Omega.(yj0-u)    */
3221         }
3222         /* Mass-weighting */
3223         wj = N_M * mj;
3224
3225         /* Calculate (xj-u) resp. (xj-xc) */
3226         rvec_sub(xj, erg->xc_center, xj_u); /* xj_u = xj-u           */
3227
3228         cprod(erg->vec, xj_u, v_xj_u); /* v_xj_u = v x (xj-u)   */
3229
3230         fac = norm2(v_xj_u);
3231         /*                                      *                      1           */
3232         psijstar = 1.0 / (fac + erg->rotg->eps); /*  psistar = --------------------  */
3233         /*                                      *            |v x (xj-u)|^2 + eps  */
3234
3235         psij = gmx::invsqrt(fac); /*                 1                */
3236                                   /*  psij    = ------------          */
3237                                   /*            |v x (xj-u)|          */
3238
3239         svmul(psij, v_xj_u, sj); /*  sj = psij * (v x (xj-u) )       */
3240
3241         fac  = iprod(v_xj_u, rj); /* fac = (v x (xj-u)).rj */
3242         fac2 = fac * fac;
3243
3244         sjrj = iprod(sj, rj); /* sjrj = sj.rj          */
3245
3246         svmul(psijstar / psij, rj, tmpvec);
3247         svmul(psijstar * psijstar / (psij * psij * psij) * sjrj, sj, tmpvec2);
3248         rvec_dec(tmpvec, tmpvec2);
3249         cprod(tmpvec, erg->vec, tmpvec2);
3250
3251         /* Store the additional force so that it can be added to the force
3252          * array after the normal forces have been evaluated */
3253         svmul(-erg->rotg->k * wj * sjrj, tmpvec2, tmpvec);
3254         svmul(mj, innersumvec, tmpvec2); /* This is != 0 only for the pivot-free variant */
3255
3256         rvec_add(tmpvec2, tmpvec, erg->f_rot_loc[j]);
3257         Vpart += wj * psijstar * fac2;
3258
3259         /* If requested, also calculate the potential for a set of angles
3260          * near the current reference angle */
3261         if (bCalcPotFit)
3262         {
3263             for (int ifit = 0; ifit < erg->rotg->PotAngle_nstep; ifit++)
3264             {
3265                 if (bPF)
3266                 {
3267                     mvmul(erg->PotAngleFit->rotmat[ifit], yj0_yc0, fit_rj); /* fit_rj = Omega.(yj0-yc0) */
3268                 }
3269                 else
3270                 {
3271                     /* Position of this atom in the collective array */
3272                     int iigrp = collectiveRotationGroupIndex[j];
3273                     /* Rotate with the alternative angle. Like rotate_local_reference(),
3274                      * just for a single local atom */
3275                     mvmul(erg->PotAngleFit->rotmat[ifit], erg->rotg->x_ref[iigrp], fit_rj); /* fit_rj = Omega*(yj0-u) */
3276                 }
3277                 fit_fac = iprod(v_xj_u, fit_rj); /* fac = (v x (xj-u)).fit_rj */
3278                 /* Add to the rotation potential for this angle: */
3279                 erg->PotAngleFit->V[ifit] += 0.5 * erg->rotg->k * wj * psijstar * fit_fac * fit_fac;
3280             }
3281         }
3282
3283         if (bOutstepRot)
3284         {
3285             /* Add to the torque of this rotation group */
3286             erg->torque_v += torque(erg->vec, erg->f_rot_loc[j], xj, erg->xc_center);
3287
3288             /* Calculate the angle between reference and actual rotation group atom. */
3289             angle(erg, xj_u, rj, &alpha, &weight); /* angle in rad, weighted */
3290             erg->angle_v += alpha * weight;
3291             erg->weight_v += weight;
3292         }
3293
3294         PRINT_FORCE_J
3295
3296     } /* end of loop over local rotation group atoms */
3297     erg->V = 0.5 * erg->rotg->k * Vpart;
3298 }
3299
3300
3301 /* Determine the smallest and largest position vector (with respect to the
3302  * rotation vector) for the reference group */
3303 static void get_firstlast_atom_ref(const gmx_enfrotgrp* erg, int* firstindex, int* lastindex)
3304 {
3305     int  i;
3306     real xcproj;           /* The projection of a reference position on the
3307                               rotation vector */
3308     real minproj, maxproj; /* Smallest and largest projection on v */
3309
3310     /* Start with some value */
3311     minproj = iprod(erg->rotg->x_ref[0], erg->vec);
3312     maxproj = minproj;
3313
3314     /* This is just to ensure that it still works if all the atoms of the
3315      * reference structure are situated in a plane perpendicular to the rotation
3316      * vector */
3317     *firstindex = 0;
3318     *lastindex  = erg->rotg->nat - 1;
3319
3320     /* Loop over all atoms of the reference group,
3321      * project them on the rotation vector to find the extremes */
3322     for (i = 0; i < erg->rotg->nat; i++)
3323     {
3324         xcproj = iprod(erg->rotg->x_ref[i], erg->vec);
3325         if (xcproj < minproj)
3326         {
3327             minproj     = xcproj;
3328             *firstindex = i;
3329         }
3330         if (xcproj > maxproj)
3331         {
3332             maxproj    = xcproj;
3333             *lastindex = i;
3334         }
3335     }
3336 }
3337
3338
3339 /* Allocate memory for the slabs */
3340 static void allocate_slabs(gmx_enfrotgrp* erg, FILE* fplog, gmx_bool bVerbose)
3341 {
3342     /* More slabs than are defined for the reference are never needed */
3343     int nslabs = erg->slab_last_ref - erg->slab_first_ref + 1;
3344
3345     /* Remember how many we allocated */
3346     erg->nslabs_alloc = nslabs;
3347
3348     if ((nullptr != fplog) && bVerbose)
3349     {
3350         fprintf(fplog,
3351                 "%s allocating memory to store data for %d slabs (rotation group %d).\n",
3352                 RotStr.c_str(),
3353                 nslabs,
3354                 erg->groupIndex);
3355     }
3356     snew(erg->slab_center, nslabs);
3357     snew(erg->slab_center_ref, nslabs);
3358     snew(erg->slab_weights, nslabs);
3359     snew(erg->slab_torque_v, nslabs);
3360     snew(erg->slab_data, nslabs);
3361     snew(erg->gn_atom, nslabs);
3362     snew(erg->gn_slabind, nslabs);
3363     snew(erg->slab_innersumvec, nslabs);
3364     for (int i = 0; i < nslabs; i++)
3365     {
3366         snew(erg->slab_data[i].x, erg->rotg->nat);
3367         snew(erg->slab_data[i].ref, erg->rotg->nat);
3368         snew(erg->slab_data[i].weight, erg->rotg->nat);
3369     }
3370     snew(erg->xc_ref_sorted, erg->rotg->nat);
3371     snew(erg->xc_sortind, erg->rotg->nat);
3372     snew(erg->firstatom, nslabs);
3373     snew(erg->lastatom, nslabs);
3374 }
3375
3376
3377 /* From the extreme positions of the reference group, determine the first
3378  * and last slab of the reference. We can never have more slabs in the real
3379  * simulation than calculated here for the reference.
3380  */
3381 static void get_firstlast_slab_ref(gmx_enfrotgrp* erg, real mc[], int ref_firstindex, int ref_lastindex)
3382 {
3383     rvec dummy;
3384
3385     int first = get_first_slab(erg, erg->rotg->x_ref[ref_firstindex]);
3386     int last  = get_last_slab(erg, erg->rotg->x_ref[ref_lastindex]);
3387
3388     while (get_slab_weight(first, erg, erg->rotg->x_ref, mc, &dummy) > WEIGHT_MIN)
3389     {
3390         first--;
3391     }
3392     erg->slab_first_ref = first + 1;
3393     while (get_slab_weight(last, erg, erg->rotg->x_ref, mc, &dummy) > WEIGHT_MIN)
3394     {
3395         last++;
3396     }
3397     erg->slab_last_ref = last - 1;
3398 }
3399
3400
3401 /* Special version of copy_rvec:
3402  * During the copy procedure of xcurr to b, the correct PBC image is chosen
3403  * such that the copied vector ends up near its reference position xref */
3404 static inline void copy_correct_pbc_image(const rvec xcurr, /* copy vector xcurr ... */
3405                                           rvec       b, /* ... to b ...                         */
3406                                           const rvec xref, /* choosing the PBC image such that b ends up near xref */
3407                                           const matrix box,
3408                                           int          npbcdim)
3409 {
3410     rvec dx;
3411     int  d, m;
3412     ivec shift;
3413
3414
3415     /* Shortest PBC distance between the atom and its reference */
3416     rvec_sub(xcurr, xref, dx);
3417
3418     /* Determine the shift for this atom */
3419     clear_ivec(shift);
3420     for (m = npbcdim - 1; m >= 0; m--)
3421     {
3422         while (dx[m] < -0.5 * box[m][m])
3423         {
3424             for (d = 0; d < DIM; d++)
3425             {
3426                 dx[d] += box[m][d];
3427             }
3428             shift[m]++;
3429         }
3430         while (dx[m] >= 0.5 * box[m][m])
3431         {
3432             for (d = 0; d < DIM; d++)
3433             {
3434                 dx[d] -= box[m][d];
3435             }
3436             shift[m]--;
3437         }
3438     }
3439
3440     /* Apply the shift to the position */
3441     copy_rvec(xcurr, b);
3442     shift_single_coord(box, b, shift);
3443 }
3444
3445
3446 static void init_rot_group(FILE*             fplog,
3447                            const t_commrec*  cr,
3448                            gmx_enfrotgrp*    erg,
3449                            rvec*             x,
3450                            const gmx_mtop_t& mtop,
3451                            gmx_bool          bVerbose,
3452                            FILE*             out_slabs,
3453                            const matrix      box,
3454                            t_inputrec*       ir,
3455                            gmx_bool          bOutputCenters)
3456 {
3457     rvec            coord, xref, *xdum;
3458     gmx_bool        bFlex, bColl;
3459     int             ref_firstindex, ref_lastindex;
3460     real            mass, totalmass;
3461     real            start = 0.0;
3462     double          t_start;
3463     const t_rotgrp* rotg = erg->rotg;
3464
3465
3466     /* Do we have a flexible axis? */
3467     bFlex = ISFLEX(rotg);
3468     /* Do we use a global set of coordinates? */
3469     bColl = ISCOLL(rotg);
3470
3471     /* Allocate space for collective coordinates if needed */
3472     if (bColl)
3473     {
3474         snew(erg->xc, erg->rotg->nat);
3475         snew(erg->xc_shifts, erg->rotg->nat);
3476         snew(erg->xc_eshifts, erg->rotg->nat);
3477         snew(erg->xc_old, erg->rotg->nat);
3478
3479         if (erg->rotg->eFittype == RotationGroupFitting::Norm)
3480         {
3481             snew(erg->xc_ref_length, erg->rotg->nat); /* in case fit type NORM is chosen */
3482             snew(erg->xc_norm, erg->rotg->nat);
3483         }
3484     }
3485     else
3486     {
3487         snew(erg->xr_loc, erg->rotg->nat);
3488         snew(erg->x_loc_pbc, erg->rotg->nat);
3489     }
3490
3491     copy_rvec(erg->rotg->inputVec, erg->vec);
3492     snew(erg->f_rot_loc, erg->rotg->nat);
3493
3494     /* Make space for the calculation of the potential at other angles (used
3495      * for fitting only) */
3496     if (RotationGroupFitting::Pot == erg->rotg->eFittype)
3497     {
3498         snew(erg->PotAngleFit, 1);
3499         snew(erg->PotAngleFit->degangle, erg->rotg->PotAngle_nstep);
3500         snew(erg->PotAngleFit->V, erg->rotg->PotAngle_nstep);
3501         snew(erg->PotAngleFit->rotmat, erg->rotg->PotAngle_nstep);
3502
3503         /* Get the set of angles around the reference angle */
3504         start = -0.5 * (erg->rotg->PotAngle_nstep - 1) * erg->rotg->PotAngle_step;
3505         for (int i = 0; i < erg->rotg->PotAngle_nstep; i++)
3506         {
3507             erg->PotAngleFit->degangle[i] = start + i * erg->rotg->PotAngle_step;
3508         }
3509     }
3510     else
3511     {
3512         erg->PotAngleFit = nullptr;
3513     }
3514
3515     /* Copy the masses so that the center can be determined. For all types of
3516      * enforced rotation, we store the masses in the erg->mc array. */
3517     snew(erg->mc, erg->rotg->nat);
3518     if (bFlex)
3519     {
3520         snew(erg->mc_sorted, erg->rotg->nat);
3521     }
3522     if (!bColl)
3523     {
3524         snew(erg->m_loc, erg->rotg->nat);
3525     }
3526     totalmass = 0.0;
3527     int molb  = 0;
3528     for (int i = 0; i < erg->rotg->nat; i++)
3529     {
3530         if (erg->rotg->bMassW)
3531         {
3532             mass = mtopGetAtomMass(mtop, erg->rotg->ind[i], &molb);
3533         }
3534         else
3535         {
3536             mass = 1.0;
3537         }
3538         erg->mc[i] = mass;
3539         totalmass += mass;
3540     }
3541     erg->invmass = 1.0 / totalmass;
3542
3543     /* Set xc_ref_center for any rotation potential */
3544     if ((erg->rotg->eType == EnforcedRotationGroupType::Iso)
3545         || (erg->rotg->eType == EnforcedRotationGroupType::Pm)
3546         || (erg->rotg->eType == EnforcedRotationGroupType::Rm)
3547         || (erg->rotg->eType == EnforcedRotationGroupType::Rm2))
3548     {
3549         /* Set the pivot point for the fixed, stationary-axis potentials. This
3550          * won't change during the simulation */
3551         copy_rvec(erg->rotg->pivot, erg->xc_ref_center);
3552         copy_rvec(erg->rotg->pivot, erg->xc_center);
3553     }
3554     else
3555     {
3556         /* Center of the reference positions */
3557         get_center(erg->rotg->x_ref, erg->mc, erg->rotg->nat, erg->xc_ref_center);
3558
3559         /* Center of the actual positions */
3560         if (MASTER(cr))
3561         {
3562             snew(xdum, erg->rotg->nat);
3563             for (int i = 0; i < erg->rotg->nat; i++)
3564             {
3565                 int ii = erg->rotg->ind[i];
3566                 copy_rvec(x[ii], xdum[i]);
3567             }
3568             get_center(xdum, erg->mc, erg->rotg->nat, erg->xc_center);
3569             sfree(xdum);
3570         }
3571 #if GMX_MPI
3572         if (PAR(cr))
3573         {
3574             gmx_bcast(sizeof(erg->xc_center), erg->xc_center, cr->mpi_comm_mygroup);
3575         }
3576 #endif
3577     }
3578
3579     if (bColl)
3580     {
3581         /* Save the original (whole) set of positions in xc_old such that at later
3582          * steps the rotation group can always be made whole again. If the simulation is
3583          * restarted, we compute the starting reference positions (given the time)
3584          * and assume that the correct PBC image of each position is the one nearest
3585          * to the current reference */
3586         if (MASTER(cr))
3587         {
3588             /* Calculate the rotation matrix for this angle: */
3589             t_start       = ir->init_t + ir->init_step * ir->delta_t;
3590             erg->degangle = erg->rotg->rate * t_start;
3591             calc_rotmat(erg->vec, erg->degangle, erg->rotmat);
3592
3593             for (int i = 0; i < erg->rotg->nat; i++)
3594             {
3595                 int ii = erg->rotg->ind[i];
3596
3597                 /* Subtract pivot, rotate, and add pivot again. This will yield the
3598                  * reference position for time t */
3599                 rvec_sub(erg->rotg->x_ref[i], erg->xc_ref_center, coord);
3600                 mvmul(erg->rotmat, coord, xref);
3601                 rvec_inc(xref, erg->xc_ref_center);
3602
3603                 copy_correct_pbc_image(x[ii], erg->xc_old[i], xref, box, 3);
3604             }
3605         }
3606 #if GMX_MPI
3607         if (PAR(cr))
3608         {
3609             gmx_bcast(erg->rotg->nat * sizeof(erg->xc_old[0]), erg->xc_old, cr->mpi_comm_mygroup);
3610         }
3611 #endif
3612     }
3613
3614     if ((erg->rotg->eType != EnforcedRotationGroupType::Flex)
3615         && (erg->rotg->eType != EnforcedRotationGroupType::Flex2))
3616     {
3617         /* Put the reference positions into origin: */
3618         for (int i = 0; i < erg->rotg->nat; i++)
3619         {
3620             rvec_dec(erg->rotg->x_ref[i], erg->xc_ref_center);
3621         }
3622     }
3623
3624     /* Enforced rotation with flexible axis */
3625     if (bFlex)
3626     {
3627         /* Calculate maximum beta value from minimum gaussian (performance opt.) */
3628         erg->max_beta = calc_beta_max(erg->rotg->min_gaussian, erg->rotg->slab_dist);
3629
3630         /* Determine the smallest and largest coordinate with respect to the rotation vector */
3631         get_firstlast_atom_ref(erg, &ref_firstindex, &ref_lastindex);
3632
3633         /* From the extreme positions of the reference group, determine the first
3634          * and last slab of the reference. */
3635         get_firstlast_slab_ref(erg, erg->mc, ref_firstindex, ref_lastindex);
3636
3637         /* Allocate memory for the slabs */
3638         allocate_slabs(erg, fplog, bVerbose);
3639
3640         /* Flexible rotation: determine the reference centers for the rest of the simulation */
3641         erg->slab_first = erg->slab_first_ref;
3642         erg->slab_last  = erg->slab_last_ref;
3643         get_slab_centers(erg, erg->rotg->x_ref, erg->mc, -1, out_slabs, bOutputCenters, TRUE);
3644
3645         /* Length of each x_rotref vector from center (needed if fit routine NORM is chosen): */
3646         if (erg->rotg->eFittype == RotationGroupFitting::Norm)
3647         {
3648             for (int i = 0; i < erg->rotg->nat; i++)
3649             {
3650                 rvec_sub(erg->rotg->x_ref[i], erg->xc_ref_center, coord);
3651                 erg->xc_ref_length[i] = norm(coord);
3652             }
3653         }
3654     }
3655 }
3656
3657 /* Calculate the size of the MPI buffer needed in reduce_output() */
3658 static int calc_mpi_bufsize(const gmx_enfrot* er)
3659
3660 {
3661     int count_total = 0;
3662     for (int g = 0; g < er->rot->ngrp; g++)
3663     {
3664         const t_rotgrp*      rotg = &er->rot->grp[g];
3665         const gmx_enfrotgrp* erg  = &er->enfrotgrp[g];
3666
3667         /* Count the items that are transferred for this group: */
3668         int count_group = 4; /* V, torque, angle, weight */
3669
3670         /* Add the maximum number of slabs for flexible groups */
3671         if (ISFLEX(rotg))
3672         {
3673             count_group += erg->slab_last_ref - erg->slab_first_ref + 1;
3674         }
3675
3676         /* Add space for the potentials at different angles: */
3677         if (RotationGroupFitting::Pot == erg->rotg->eFittype)
3678         {
3679             count_group += erg->rotg->PotAngle_nstep;
3680         }
3681
3682         /* Add to the total number: */
3683         count_total += count_group;
3684     }
3685
3686     return count_total;
3687 }
3688
3689
3690 std::unique_ptr<gmx::EnforcedRotation> init_rot(FILE*                       fplog,
3691                                                 t_inputrec*                 ir,
3692                                                 int                         nfile,
3693                                                 const t_filenm              fnm[],
3694                                                 const t_commrec*            cr,
3695                                                 gmx::LocalAtomSetManager*   atomSets,
3696                                                 const t_state*              globalState,
3697                                                 const gmx_mtop_t&           mtop,
3698                                                 const gmx_output_env_t*     oenv,
3699                                                 const gmx::MdrunOptions&    mdrunOptions,
3700                                                 const gmx::StartingBehavior startingBehavior)
3701 {
3702     int   nat_max = 0;       /* Size of biggest rotation group */
3703     rvec* x_pbc   = nullptr; /* Space for the pbc-correct atom positions */
3704
3705     if (MASTER(cr) && mdrunOptions.verbose)
3706     {
3707         fprintf(stdout, "%s Initializing ...\n", RotStr.c_str());
3708     }
3709
3710     auto        enforcedRotation = std::make_unique<gmx::EnforcedRotation>();
3711     gmx_enfrot* er               = enforcedRotation->getLegacyEnfrot();
3712     // TODO When this module implements IMdpOptions, the ownership will become more clear.
3713     er->rot                  = ir->rot;
3714     er->restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
3715
3716     /* When appending, skip first output to avoid duplicate entries in the data files */
3717     er->bOut = er->restartWithAppending;
3718
3719     if (MASTER(cr) && er->bOut)
3720     {
3721         please_cite(fplog, "Kutzner2011");
3722     }
3723
3724     /* Output every step for reruns */
3725     if (mdrunOptions.rerun)
3726     {
3727         if (nullptr != fplog)
3728         {
3729             fprintf(fplog, "%s rerun - will write rotation output every available step.\n", RotStr.c_str());
3730         }
3731         er->nstrout = 1;
3732         er->nstsout = 1;
3733     }
3734     else
3735     {
3736         er->nstrout = er->rot->nstrout;
3737         er->nstsout = er->rot->nstsout;
3738     }
3739
3740     er->out_slabs = nullptr;
3741     if (MASTER(cr) && HaveFlexibleGroups(er->rot))
3742     {
3743         er->out_slabs = open_slab_out(opt2fn("-rs", nfile, fnm), er);
3744     }
3745
3746     if (MASTER(cr))
3747     {
3748         /* Remove pbc, make molecule whole.
3749          * When ir->bContinuation=TRUE this has already been done, but ok. */
3750         snew(x_pbc, mtop.natoms);
3751         copy_rvecn(globalState->x.rvec_array(), x_pbc, 0, mtop.natoms);
3752         do_pbc_first_mtop(nullptr, ir->pbcType, globalState->box, &mtop, x_pbc);
3753         /* All molecules will be whole now, but not necessarily in the home box.
3754          * Additionally, if a rotation group consists of more than one molecule
3755          * (e.g. two strands of DNA), each one of them can end up in a different
3756          * periodic box. This is taken care of in init_rot_group.  */
3757     }
3758
3759     /* Allocate space for the per-rotation-group data: */
3760     er->enfrotgrp.resize(er->rot->ngrp);
3761     int groupIndex = 0;
3762     for (auto& ergRef : er->enfrotgrp)
3763     {
3764         gmx_enfrotgrp* erg = &ergRef;
3765         erg->rotg          = &er->rot->grp[groupIndex];
3766         erg->atomSet       = std::make_unique<gmx::LocalAtomSet>(
3767                 atomSets->add({ erg->rotg->ind, erg->rotg->ind + erg->rotg->nat }));
3768         erg->groupIndex = groupIndex;
3769
3770         if (nullptr != fplog)
3771         {
3772             fprintf(fplog,
3773                     "%s group %d type '%s'\n",
3774                     RotStr.c_str(),
3775                     groupIndex,
3776                     enumValueToString(erg->rotg->eType));
3777         }
3778
3779         if (erg->rotg->nat > 0)
3780         {
3781             nat_max = std::max(nat_max, erg->rotg->nat);
3782
3783             init_rot_group(fplog,
3784                            cr,
3785                            erg,
3786                            x_pbc,
3787                            mtop,
3788                            mdrunOptions.verbose,
3789                            er->out_slabs,
3790                            MASTER(cr) ? globalState->box : nullptr,
3791                            ir,
3792                            !er->restartWithAppending); /* Do not output the reference centers
3793                                                         * again if we are appending */
3794         }
3795         ++groupIndex;
3796     }
3797
3798     /* Allocate space for enforced rotation buffer variables */
3799     er->bufsize = nat_max;
3800     snew(er->data, nat_max);
3801     snew(er->xbuf, nat_max);
3802     snew(er->mbuf, nat_max);
3803
3804     /* Buffers for MPI reducing torques, angles, weights (for each group), and V */
3805     if (PAR(cr))
3806     {
3807         er->mpi_bufsize = calc_mpi_bufsize(er) + 100; /* larger to catch errors */
3808         snew(er->mpi_inbuf, er->mpi_bufsize);
3809         snew(er->mpi_outbuf, er->mpi_bufsize);
3810     }
3811     else
3812     {
3813         er->mpi_bufsize = 0;
3814         er->mpi_inbuf   = nullptr;
3815         er->mpi_outbuf  = nullptr;
3816     }
3817
3818     /* Only do I/O on the MASTER */
3819     er->out_angles = nullptr;
3820     er->out_rot    = nullptr;
3821     er->out_torque = nullptr;
3822     if (MASTER(cr))
3823     {
3824         er->out_rot = open_rot_out(opt2fn("-ro", nfile, fnm), oenv, er);
3825
3826         if (er->nstsout > 0)
3827         {
3828             if (HaveFlexibleGroups(er->rot) || HavePotFitGroups(er->rot))
3829             {
3830                 er->out_angles = open_angles_out(opt2fn("-ra", nfile, fnm), er);
3831             }
3832             if (HaveFlexibleGroups(er->rot))
3833             {
3834                 er->out_torque = open_torque_out(opt2fn("-rt", nfile, fnm), er);
3835             }
3836         }
3837
3838         sfree(x_pbc);
3839     }
3840     return enforcedRotation;
3841 }
3842
3843 /* Rotate the local reference positions and store them in
3844  * erg->xr_loc[0...(nat_loc-1)]
3845  *
3846  * Note that we already subtracted u or y_c from the reference positions
3847  * in init_rot_group().
3848  */
3849 static void rotate_local_reference(gmx_enfrotgrp* erg)
3850 {
3851     const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
3852     for (size_t i = 0; i < erg->atomSet->numAtomsLocal(); i++)
3853     {
3854         /* Index of this rotation group atom with respect to the whole rotation group */
3855         int ii = collectiveRotationGroupIndex[i];
3856         /* Rotate */
3857         mvmul(erg->rotmat, erg->rotg->x_ref[ii], erg->xr_loc[i]);
3858     }
3859 }
3860
3861
3862 /* Select the PBC representation for each local x position and store that
3863  * for later usage. We assume the right PBC image of an x is the one nearest to
3864  * its rotated reference */
3865 static void choose_pbc_image(gmx::ArrayRef<const gmx::RVec> coords, gmx_enfrotgrp* erg, const matrix box, int npbcdim)
3866 {
3867     const auto& localRotationGroupIndex = erg->atomSet->localIndex();
3868     for (gmx::index i = 0; i < localRotationGroupIndex.ssize(); i++)
3869     {
3870         /* Index of a rotation group atom  */
3871         int ii = localRotationGroupIndex[i];
3872
3873         /* Get the correctly rotated reference position. The pivot was already
3874          * subtracted in init_rot_group() from the reference positions. Also,
3875          * the reference positions have already been rotated in
3876          * rotate_local_reference(). For the current reference position we thus
3877          * only need to add the pivot again. */
3878         rvec xref;
3879         copy_rvec(erg->xr_loc[i], xref);
3880         rvec_inc(xref, erg->xc_ref_center);
3881
3882         copy_correct_pbc_image(coords[ii], erg->x_loc_pbc[i], xref, box, npbcdim);
3883     }
3884 }
3885
3886
3887 void do_rotation(const t_commrec*               cr,
3888                  gmx_enfrot*                    er,
3889                  const matrix                   box,
3890                  gmx::ArrayRef<const gmx::RVec> coords,
3891                  real                           t,
3892                  int64_t                        step,
3893                  gmx_bool                       bNS)
3894 {
3895     gmx_bool    outstep_slab, outstep_rot;
3896     gmx_bool    bColl;
3897     rvec        transvec;
3898     gmx_potfit* fit = nullptr; /* For fit type 'potential' determine the fit
3899                                     angle via the potential minimum            */
3900
3901 #ifdef TAKETIME
3902     double t0;
3903 #endif
3904
3905     /* When to output in main rotation output file */
3906     outstep_rot = do_per_step(step, er->nstrout) && er->bOut;
3907     /* When to output per-slab data */
3908     outstep_slab = do_per_step(step, er->nstsout) && er->bOut;
3909
3910     /* Output time into rotation output file */
3911     if (outstep_rot && MASTER(cr))
3912     {
3913         fprintf(er->out_rot, "%12.3e", t);
3914     }
3915
3916     /**************************************************************************/
3917     /* First do ALL the communication! */
3918     for (auto& ergRef : er->enfrotgrp)
3919     {
3920         gmx_enfrotgrp*  erg  = &ergRef;
3921         const t_rotgrp* rotg = erg->rotg;
3922
3923         /* Do we use a collective (global) set of coordinates? */
3924         bColl = ISCOLL(rotg);
3925
3926         /* Calculate the rotation matrix for this angle: */
3927         erg->degangle = rotg->rate * t;
3928         calc_rotmat(erg->vec, erg->degangle, erg->rotmat);
3929
3930         if (bColl)
3931         {
3932             /* Transfer the rotation group's positions such that every node has
3933              * all of them. Every node contributes its local positions x and stores
3934              * it in the collective erg->xc array. */
3935             communicate_group_positions(cr,
3936                                         erg->xc,
3937                                         erg->xc_shifts,
3938                                         erg->xc_eshifts,
3939                                         bNS,
3940                                         as_rvec_array(coords.data()),
3941                                         rotg->nat,
3942                                         erg->atomSet->numAtomsLocal(),
3943                                         erg->atomSet->localIndex().data(),
3944                                         erg->atomSet->collectiveIndex().data(),
3945                                         erg->xc_old,
3946                                         box);
3947         }
3948         else
3949         {
3950             /* Fill the local masses array;
3951              * this array changes in DD/neighborsearching steps */
3952             if (bNS)
3953             {
3954                 const auto& collectiveRotationGroupIndex = erg->atomSet->collectiveIndex();
3955                 for (gmx::index i = 0; i < collectiveRotationGroupIndex.ssize(); i++)
3956                 {
3957                     /* Index of local atom w.r.t. the collective rotation group */
3958                     int ii        = collectiveRotationGroupIndex[i];
3959                     erg->m_loc[i] = erg->mc[ii];
3960                 }
3961             }
3962
3963             /* Calculate Omega*(y_i-y_c) for the local positions */
3964             rotate_local_reference(erg);
3965
3966             /* Choose the nearest PBC images of the group atoms with respect
3967              * to the rotated reference positions */
3968             choose_pbc_image(coords, erg, box, 3);
3969
3970             /* Get the center of the rotation group */
3971             if ((rotg->eType == EnforcedRotationGroupType::Isopf)
3972                 || (rotg->eType == EnforcedRotationGroupType::Pmpf))
3973             {
3974                 get_center_comm(
3975                         cr, erg->x_loc_pbc, erg->m_loc, erg->atomSet->numAtomsLocal(), rotg->nat, erg->xc_center);
3976             }
3977         }
3978
3979     } /* End of loop over rotation groups */
3980
3981     /**************************************************************************/
3982     /* Done communicating, we can start to count cycles for the load balancing now ... */
3983     if (DOMAINDECOMP(cr))
3984     {
3985         ddReopenBalanceRegionCpu(cr->dd);
3986     }
3987
3988 #ifdef TAKETIME
3989     t0 = MPI_Wtime();
3990 #endif
3991
3992     for (auto& ergRef : er->enfrotgrp)
3993     {
3994         gmx_enfrotgrp*  erg  = &ergRef;
3995         const t_rotgrp* rotg = erg->rotg;
3996
3997         if (outstep_rot && MASTER(cr))
3998         {
3999             fprintf(er->out_rot, "%12.4f", erg->degangle);
4000         }
4001
4002         /* Calculate angles and rotation matrices for potential fitting: */
4003         if ((outstep_rot || outstep_slab) && (RotationGroupFitting::Pot == rotg->eFittype))
4004         {
4005             fit = erg->PotAngleFit;
4006             for (int i = 0; i < rotg->PotAngle_nstep; i++)
4007             {
4008                 calc_rotmat(erg->vec, erg->degangle + fit->degangle[i], fit->rotmat[i]);
4009
4010                 /* Clear value from last step */
4011                 erg->PotAngleFit->V[i] = 0.0;
4012             }
4013         }
4014
4015         /* Clear values from last time step */
4016         erg->V        = 0.0;
4017         erg->torque_v = 0.0;
4018         erg->angle_v  = 0.0;
4019         erg->weight_v = 0.0;
4020
4021         switch (rotg->eType)
4022         {
4023             case EnforcedRotationGroupType::Iso:
4024             case EnforcedRotationGroupType::Isopf:
4025             case EnforcedRotationGroupType::Pm:
4026             case EnforcedRotationGroupType::Pmpf: do_fixed(erg, outstep_rot, outstep_slab); break;
4027             case EnforcedRotationGroupType::Rm:
4028                 do_radial_motion(erg, outstep_rot, outstep_slab);
4029                 break;
4030             case EnforcedRotationGroupType::Rmpf:
4031                 do_radial_motion_pf(erg, coords, box, outstep_rot, outstep_slab);
4032                 break;
4033             case EnforcedRotationGroupType::Rm2:
4034             case EnforcedRotationGroupType::Rm2pf:
4035                 do_radial_motion2(erg, coords, box, outstep_rot, outstep_slab);
4036                 break;
4037             case EnforcedRotationGroupType::Flext:
4038             case EnforcedRotationGroupType::Flex2t:
4039                 /* Subtract the center of the rotation group from the collective positions array
4040                  * Also store the center in erg->xc_center since it needs to be subtracted
4041                  * in the low level routines from the local coordinates as well */
4042                 get_center(erg->xc, erg->mc, rotg->nat, erg->xc_center);
4043                 svmul(-1.0, erg->xc_center, transvec);
4044                 translate_x(erg->xc, rotg->nat, transvec);
4045                 do_flexible(MASTER(cr), er, erg, coords, box, t, outstep_rot, outstep_slab);
4046                 break;
4047             case EnforcedRotationGroupType::Flex:
4048             case EnforcedRotationGroupType::Flex2:
4049                 /* Do NOT subtract the center of mass in the low level routines! */
4050                 clear_rvec(erg->xc_center);
4051                 do_flexible(MASTER(cr), er, erg, coords, box, t, outstep_rot, outstep_slab);
4052                 break;
4053             default: gmx_fatal(FARGS, "No such rotation potential.");
4054         }
4055     }
4056
4057 #ifdef TAKETIME
4058     if (MASTER(cr))
4059     {
4060         fprintf(stderr, "%s calculation (step %d) took %g seconds.\n", RotStr, step, MPI_Wtime() - t0);
4061     }
4062 #endif
4063 }