src/gromacs/mdlib/pull_rotation.c

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