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