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