src/mdlib/pull_rotation.c

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