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