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