src/mdlib/pull_rotation.c

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