src/mdlib/pull_rotation.c

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