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