2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2008, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "pull_rotation.h"
49 #include "gromacs/commandline/filenm.h"
50 #include "gromacs/domdec/dlbtiming.h"
51 #include "gromacs/domdec/domdec_struct.h"
52 #include "gromacs/domdec/ga2la.h"
53 #include "gromacs/fileio/gmxfio.h"
54 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/gmxlib/network.h"
56 #include "gromacs/linearalgebra/nrjac.h"
57 #include "gromacs/math/functions.h"
58 #include "gromacs/math/utilities.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/mdlib/groupcoord.h"
61 #include "gromacs/mdlib/mdrun.h"
62 #include "gromacs/mdlib/sim_util.h"
63 #include "gromacs/mdtypes/inputrec.h"
64 #include "gromacs/mdtypes/md_enums.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/timing/cyclecounter.h"
67 #include "gromacs/timing/wallcycle.h"
68 #include "gromacs/topology/mtop_lookup.h"
69 #include "gromacs/topology/mtop_util.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/pleasecite.h"
72 #include "gromacs/utility/qsort_threadsafe.h"
73 #include "gromacs/utility/smalloc.h"
75 static char const *RotStr
= {"Enforced rotation:"};
77 /* Set the minimum weight for the determination of the slab centers */
78 #define WEIGHT_MIN (10*GMX_FLOAT_MIN)
80 /* Helper structure for sorting positions along rotation vector */
82 real xcproj
; /* Projection of xc on the rotation vector */
83 int ind
; /* Index of xc */
85 rvec x
; /* Position */
86 rvec x_ref
; /* Reference position */
90 /* Enforced rotation / flexible: determine the angle of each slab */
91 typedef struct gmx_slabdata
93 int nat
; /* Number of atoms belonging to this slab */
94 rvec
*x
; /* The positions belonging to this slab. In
95 general, this should be all positions of the
96 whole rotation group, but we leave those away
97 that have a small enough weight */
98 rvec
*ref
; /* Same for reference */
99 real
*weight
; /* The weight for each atom */
103 /* Helper structure for potential fitting */
104 typedef struct gmx_potfit
106 real
*degangle
; /* Set of angles for which the potential is
107 calculated. The optimum fit is determined as
108 the angle for with the potential is minimal */
109 real
*V
; /* Potential for the different angles */
110 matrix
*rotmat
; /* Rotation matrix corresponding to the angles */
114 /* Enforced rotation data for all groups */
115 typedef struct gmx_enfrot
117 FILE *out_rot
; /* Output file for rotation data */
118 FILE *out_torque
; /* Output file for torque data */
119 FILE *out_angles
; /* Output file for slab angles for flexible type */
120 FILE *out_slabs
; /* Output file for slab centers */
121 int bufsize
; /* Allocation size of buf */
122 rvec
*xbuf
; /* Coordinate buffer variable for sorting */
123 real
*mbuf
; /* Masses buffer variable for sorting */
124 sort_along_vec_t
*data
; /* Buffer variable needed for position sorting */
125 real
*mpi_inbuf
; /* MPI buffer */
126 real
*mpi_outbuf
; /* MPI buffer */
127 int mpi_bufsize
; /* Allocation size of in & outbuf */
128 unsigned long Flags
; /* mdrun flags */
129 gmx_bool bOut
; /* Used to skip first output when appending to
130 * avoid duplicate entries in rotation outfiles */
134 /* Global enforced rotation data for a single rotation group */
135 typedef struct gmx_enfrotgrp
137 real degangle
; /* Rotation angle in degrees */
138 matrix rotmat
; /* Rotation matrix */
139 int *ind_loc
; /* Local rotation indices */
140 int nat_loc
; /* Number of local group atoms */
141 int nalloc_loc
; /* Allocation size for ind_loc and weight_loc */
143 real V
; /* Rotation potential for this rotation group */
144 rvec
*f_rot_loc
; /* Array to store the forces on the local atoms
145 resulting from enforced rotation potential */
147 /* Collective coordinates for the whole rotation group */
148 real
*xc_ref_length
; /* Length of each x_rotref vector after x_rotref
149 has been put into origin */
150 int *xc_ref_ind
; /* Position of each local atom in the collective
152 rvec xc_center
; /* Center of the rotation group positions, may
154 rvec xc_ref_center
; /* dito, for the reference positions */
155 rvec
*xc
; /* Current (collective) positions */
156 ivec
*xc_shifts
; /* Current (collective) shifts */
157 ivec
*xc_eshifts
; /* Extra shifts since last DD step */
158 rvec
*xc_old
; /* Old (collective) positions */
159 rvec
*xc_norm
; /* Normalized form of the current positions */
160 rvec
*xc_ref_sorted
; /* Reference positions (sorted in the same order
161 as xc when sorted) */
162 int *xc_sortind
; /* Where is a position found after sorting? */
163 real
*mc
; /* Collective masses */
165 real invmass
; /* one over the total mass of the rotation group */
167 real torque_v
; /* Torque in the direction of rotation vector */
168 real angle_v
; /* Actual angle of the whole rotation group */
169 /* Fixed rotation only */
170 real weight_v
; /* Weights for angle determination */
171 rvec
*xr_loc
; /* Local reference coords, correctly rotated */
172 rvec
*x_loc_pbc
; /* Local current coords, correct PBC image */
173 real
*m_loc
; /* Masses of the current local atoms */
175 /* Flexible rotation only */
176 int nslabs_alloc
; /* For this many slabs memory is allocated */
177 int slab_first
; /* Lowermost slab for that the calculation needs
178 to be performed at a given time step */
179 int slab_last
; /* Uppermost slab ... */
180 int slab_first_ref
; /* First slab for which ref. center is stored */
181 int slab_last_ref
; /* Last ... */
182 int slab_buffer
; /* Slab buffer region around reference slabs */
183 int *firstatom
; /* First relevant atom for a slab */
184 int *lastatom
; /* Last relevant atom for a slab */
185 rvec
*slab_center
; /* Gaussian-weighted slab center */
186 rvec
*slab_center_ref
; /* Gaussian-weighted slab center for the
187 reference positions */
188 real
*slab_weights
; /* Sum of gaussian weights in a slab */
189 real
*slab_torque_v
; /* Torque T = r x f for each slab. */
190 /* torque_v = m.v = angular momentum in the
192 real max_beta
; /* min_gaussian from inputrec->rotgrp is the
193 minimum value the gaussian must have so that
194 the force is actually evaluated max_beta is
195 just another way to put it */
196 real
*gn_atom
; /* Precalculated gaussians for a single atom */
197 int *gn_slabind
; /* Tells to which slab each precalculated gaussian
199 rvec
*slab_innersumvec
; /* Inner sum of the flexible2 potential per slab;
200 this is precalculated for optimization reasons */
201 t_gmx_slabdata
*slab_data
; /* Holds atom positions and gaussian weights
202 of atoms belonging to a slab */
204 /* For potential fits with varying angle: */
205 t_gmx_potfit
*PotAngleFit
; /* Used for fit type 'potential' */
209 /* Activate output of forces for correctness checks */
210 /* #define PRINT_FORCES */
212 #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]);
213 #define PRINT_POT_TAU if (MASTER(cr)) { \
214 fprintf(stderr, "potential = %15.8f\n" "torque = %15.8f\n", erg->V, erg->torque_v); \
217 #define PRINT_FORCE_J
218 #define PRINT_POT_TAU
221 /* Shortcuts for often used queries */
222 #define ISFLEX(rg) ( (rg->eType == erotgFLEX) || (rg->eType == erotgFLEXT) || (rg->eType == erotgFLEX2) || (rg->eType == erotgFLEX2T) )
223 #define ISCOLL(rg) ( (rg->eType == erotgFLEX) || (rg->eType == erotgFLEXT) || (rg->eType == erotgFLEX2) || (rg->eType == erotgFLEX2T) || (rg->eType == erotgRMPF) || (rg->eType == erotgRM2PF) )
226 /* Does any of the rotation groups use slab decomposition? */
227 static gmx_bool
HaveFlexibleGroups(t_rot
*rot
)
233 for (g
= 0; g
< rot
->ngrp
; g
++)
246 /* Is for any group the fit angle determined by finding the minimum of the
247 * rotation potential? */
248 static gmx_bool
HavePotFitGroups(t_rot
*rot
)
254 for (g
= 0; g
< rot
->ngrp
; g
++)
257 if (erotgFitPOT
== rotg
->eFittype
)
267 static double** allocate_square_matrix(int dim
)
270 double** mat
= nullptr;
274 for (i
= 0; i
< dim
; i
++)
283 static void free_square_matrix(double** mat
, int dim
)
288 for (i
= 0; i
< dim
; i
++)
296 /* Return the angle for which the potential is minimal */
297 static real
get_fitangle(t_rotgrp
*rotg
, gmx_enfrotgrp_t erg
)
300 real fitangle
= -999.9;
301 real pot_min
= GMX_FLOAT_MAX
;
305 fit
= erg
->PotAngleFit
;
307 for (i
= 0; i
< rotg
->PotAngle_nstep
; i
++)
309 if (fit
->V
[i
] < pot_min
)
312 fitangle
= fit
->degangle
[i
];
320 /* Reduce potential angle fit data for this group at this time step? */
321 static gmx_inline gmx_bool
bPotAngle(t_rot
*rot
, t_rotgrp
*rotg
, gmx_int64_t step
)
323 return ( (erotgFitPOT
== rotg
->eFittype
) && (do_per_step(step
, rot
->nstsout
) || do_per_step(step
, rot
->nstrout
)) );
326 /* Reduce slab torqe data for this group at this time step? */
327 static gmx_inline gmx_bool
bSlabTau(t_rot
*rot
, t_rotgrp
*rotg
, gmx_int64_t step
)
329 return ( (ISFLEX(rotg
)) && do_per_step(step
, rot
->nstsout
) );
332 /* Output rotation energy, torques, etc. for each rotation group */
333 static void reduce_output(t_commrec
*cr
, t_rot
*rot
, real t
, gmx_int64_t step
)
335 int g
, i
, islab
, nslabs
= 0;
336 int count
; /* MPI element counter */
338 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
339 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
346 /* Fill the MPI buffer with stuff to reduce. If items are added for reduction
347 * here, the MPI buffer size has to be enlarged also in calc_mpi_bufsize() */
351 for (g
= 0; g
< rot
->ngrp
; g
++)
354 erg
= rotg
->enfrotgrp
;
355 nslabs
= erg
->slab_last
- erg
->slab_first
+ 1;
356 er
->mpi_inbuf
[count
++] = erg
->V
;
357 er
->mpi_inbuf
[count
++] = erg
->torque_v
;
358 er
->mpi_inbuf
[count
++] = erg
->angle_v
;
359 er
->mpi_inbuf
[count
++] = erg
->weight_v
; /* weights are not needed for flex types, but this is just a single value */
361 if (bPotAngle(rot
, rotg
, step
))
363 for (i
= 0; i
< rotg
->PotAngle_nstep
; i
++)
365 er
->mpi_inbuf
[count
++] = erg
->PotAngleFit
->V
[i
];
368 if (bSlabTau(rot
, rotg
, step
))
370 for (i
= 0; i
< nslabs
; i
++)
372 er
->mpi_inbuf
[count
++] = erg
->slab_torque_v
[i
];
376 if (count
> er
->mpi_bufsize
)
378 gmx_fatal(FARGS
, "%s MPI buffer overflow, please report this error.", RotStr
);
382 MPI_Reduce(er
->mpi_inbuf
, er
->mpi_outbuf
, count
, GMX_MPI_REAL
, MPI_SUM
, MASTERRANK(cr
), cr
->mpi_comm_mygroup
);
385 /* Copy back the reduced data from the buffer on the master */
389 for (g
= 0; g
< rot
->ngrp
; g
++)
392 erg
= rotg
->enfrotgrp
;
393 nslabs
= erg
->slab_last
- erg
->slab_first
+ 1;
394 erg
->V
= er
->mpi_outbuf
[count
++];
395 erg
->torque_v
= er
->mpi_outbuf
[count
++];
396 erg
->angle_v
= er
->mpi_outbuf
[count
++];
397 erg
->weight_v
= er
->mpi_outbuf
[count
++];
399 if (bPotAngle(rot
, rotg
, step
))
401 for (i
= 0; i
< rotg
->PotAngle_nstep
; i
++)
403 erg
->PotAngleFit
->V
[i
] = er
->mpi_outbuf
[count
++];
406 if (bSlabTau(rot
, rotg
, step
))
408 for (i
= 0; i
< nslabs
; i
++)
410 erg
->slab_torque_v
[i
] = er
->mpi_outbuf
[count
++];
420 /* Angle and torque for each rotation group */
421 for (g
= 0; g
< rot
->ngrp
; g
++)
424 bFlex
= ISFLEX(rotg
);
426 erg
= rotg
->enfrotgrp
;
428 /* Output to main rotation output file: */
429 if (do_per_step(step
, rot
->nstrout
) )
431 if (erotgFitPOT
== rotg
->eFittype
)
433 fitangle
= get_fitangle(rotg
, erg
);
439 fitangle
= erg
->angle_v
; /* RMSD fit angle */
443 fitangle
= (erg
->angle_v
/erg
->weight_v
)*180.0*M_1_PI
;
446 fprintf(er
->out_rot
, "%12.4f", fitangle
);
447 fprintf(er
->out_rot
, "%12.3e", erg
->torque_v
);
448 fprintf(er
->out_rot
, "%12.3e", erg
->V
);
451 if (do_per_step(step
, rot
->nstsout
) )
453 /* Output to torque log file: */
456 fprintf(er
->out_torque
, "%12.3e%6d", t
, g
);
457 for (i
= erg
->slab_first
; i
<= erg
->slab_last
; i
++)
459 islab
= i
- erg
->slab_first
; /* slab index */
460 /* Only output if enough weight is in slab */
461 if (erg
->slab_weights
[islab
] > rotg
->min_gaussian
)
463 fprintf(er
->out_torque
, "%6d%12.3e", i
, erg
->slab_torque_v
[islab
]);
466 fprintf(er
->out_torque
, "\n");
469 /* Output to angles log file: */
470 if (erotgFitPOT
== rotg
->eFittype
)
472 fprintf(er
->out_angles
, "%12.3e%6d%12.4f", t
, g
, erg
->degangle
);
473 /* Output energies at a set of angles around the reference angle */
474 for (i
= 0; i
< rotg
->PotAngle_nstep
; i
++)
476 fprintf(er
->out_angles
, "%12.3e", erg
->PotAngleFit
->V
[i
]);
478 fprintf(er
->out_angles
, "\n");
482 if (do_per_step(step
, rot
->nstrout
) )
484 fprintf(er
->out_rot
, "\n");
490 /* Add the forces from enforced rotation potential to the local forces.
491 * Should be called after the SR forces have been evaluated */
492 extern real
add_rot_forces(t_rot
*rot
, rvec f
[], t_commrec
*cr
, gmx_int64_t step
, real t
)
496 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
497 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
498 real Vrot
= 0.0; /* If more than one rotation group is present, Vrot
499 assembles the local parts from all groups */
504 /* Loop over enforced rotation groups (usually 1, though)
505 * Apply the forces from rotation potentials */
506 for (g
= 0; g
< rot
->ngrp
; g
++)
509 erg
= rotg
->enfrotgrp
;
510 Vrot
+= erg
->V
; /* add the local parts from the nodes */
511 for (l
= 0; l
< erg
->nat_loc
; l
++)
513 /* Get the right index of the local force */
514 ii
= erg
->ind_loc
[l
];
516 rvec_inc(f
[ii
], erg
->f_rot_loc
[l
]);
520 /* Reduce energy,torque, angles etc. to get the sum values (per rotation group)
521 * on the master and output these values to file. */
522 if ( (do_per_step(step
, rot
->nstrout
) || do_per_step(step
, rot
->nstsout
)) && er
->bOut
)
524 reduce_output(cr
, rot
, t
, step
);
527 /* When appending, er->bOut is FALSE the first time to avoid duplicate entries */
536 /* The Gaussian norm is chosen such that the sum of the gaussian functions
537 * over the slabs is approximately 1.0 everywhere */
538 #define GAUSS_NORM 0.569917543430618
541 /* Calculate the maximum beta that leads to a gaussian larger min_gaussian,
542 * also does some checks
544 static double calc_beta_max(real min_gaussian
, real slab_dist
)
550 /* Actually the next two checks are already made in grompp */
553 gmx_fatal(FARGS
, "Slab distance of flexible rotation groups must be >=0 !");
555 if (min_gaussian
<= 0)
557 gmx_fatal(FARGS
, "Cutoff value for Gaussian must be > 0. (You requested %f)");
560 /* Define the sigma value */
561 sigma
= 0.7*slab_dist
;
563 /* Calculate the argument for the logarithm and check that the log() result is negative or 0 */
564 arg
= min_gaussian
/GAUSS_NORM
;
567 gmx_fatal(FARGS
, "min_gaussian of flexible rotation groups must be <%g", GAUSS_NORM
);
570 return std::sqrt(-2.0*sigma
*sigma
*log(min_gaussian
/GAUSS_NORM
));
574 static gmx_inline real
calc_beta(rvec curr_x
, t_rotgrp
*rotg
, int n
)
576 return iprod(curr_x
, rotg
->vec
) - rotg
->slab_dist
* n
;
580 static gmx_inline real
gaussian_weight(rvec curr_x
, t_rotgrp
*rotg
, int n
)
582 const real norm
= GAUSS_NORM
;
586 /* Define the sigma value */
587 sigma
= 0.7*rotg
->slab_dist
;
588 /* Calculate the Gaussian value of slab n for position curr_x */
589 return norm
* exp( -0.5 * gmx::square( calc_beta(curr_x
, rotg
, n
)/sigma
) );
593 /* Returns the weight in a single slab, also calculates the Gaussian- and mass-
594 * weighted sum of positions for that slab */
595 static real
get_slab_weight(int j
, t_rotgrp
*rotg
, rvec xc
[], real mc
[], rvec
*x_weighted_sum
)
597 rvec curr_x
; /* The position of an atom */
598 rvec curr_x_weighted
; /* The gaussian-weighted position */
599 real gaussian
; /* A single gaussian weight */
600 real wgauss
; /* gaussian times current mass */
601 real slabweight
= 0.0; /* The sum of weights in the slab */
605 clear_rvec(*x_weighted_sum
);
607 /* Loop over all atoms in the rotation group */
608 for (i
= 0; i
< rotg
->nat
; i
++)
610 copy_rvec(xc
[i
], curr_x
);
611 gaussian
= gaussian_weight(curr_x
, rotg
, j
);
612 wgauss
= gaussian
* mc
[i
];
613 svmul(wgauss
, curr_x
, curr_x_weighted
);
614 rvec_add(*x_weighted_sum
, curr_x_weighted
, *x_weighted_sum
);
615 slabweight
+= wgauss
;
616 } /* END of loop over rotation group atoms */
622 static void get_slab_centers(
623 t_rotgrp
*rotg
, /* The rotation group information */
624 rvec
*xc
, /* The rotation group positions; will
625 typically be enfrotgrp->xc, but at first call
626 it is enfrotgrp->xc_ref */
627 real
*mc
, /* The masses of the rotation group atoms */
628 int g
, /* The number of the rotation group */
629 real time
, /* Used for output only */
630 FILE *out_slabs
, /* For outputting center per slab information */
631 gmx_bool bOutStep
, /* Is this an output step? */
632 gmx_bool bReference
) /* If this routine is called from
633 init_rot_group we need to store
634 the reference slab centers */
638 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
641 erg
= rotg
->enfrotgrp
;
643 /* Loop over slabs */
644 for (j
= erg
->slab_first
; j
<= erg
->slab_last
; j
++)
646 islab
= j
- erg
->slab_first
;
647 erg
->slab_weights
[islab
] = get_slab_weight(j
, rotg
, xc
, mc
, &erg
->slab_center
[islab
]);
649 /* We can do the calculations ONLY if there is weight in the slab! */
650 if (erg
->slab_weights
[islab
] > WEIGHT_MIN
)
652 svmul(1.0/erg
->slab_weights
[islab
], erg
->slab_center
[islab
], erg
->slab_center
[islab
]);
656 /* We need to check this here, since we divide through slab_weights
657 * in the flexible low-level routines! */
658 gmx_fatal(FARGS
, "Not enough weight in slab %d. Slab center cannot be determined!", j
);
661 /* At first time step: save the centers of the reference structure */
664 copy_rvec(erg
->slab_center
[islab
], erg
->slab_center_ref
[islab
]);
666 } /* END of loop over slabs */
668 /* Output on the master */
669 if ( (nullptr != out_slabs
) && bOutStep
)
671 fprintf(out_slabs
, "%12.3e%6d", time
, g
);
672 for (j
= erg
->slab_first
; j
<= erg
->slab_last
; j
++)
674 islab
= j
- erg
->slab_first
;
675 fprintf(out_slabs
, "%6d%12.3e%12.3e%12.3e",
676 j
, erg
->slab_center
[islab
][XX
], erg
->slab_center
[islab
][YY
], erg
->slab_center
[islab
][ZZ
]);
678 fprintf(out_slabs
, "\n");
683 static void calc_rotmat(
685 real degangle
, /* Angle alpha of rotation at time t in degrees */
686 matrix rotmat
) /* Rotation matrix */
688 real radangle
; /* Rotation angle in radians */
689 real cosa
; /* cosine alpha */
690 real sina
; /* sine alpha */
691 real OMcosa
; /* 1 - cos(alpha) */
692 real dumxy
, dumxz
, dumyz
; /* save computations */
693 rvec rot_vec
; /* Rotate around rot_vec ... */
696 radangle
= degangle
* M_PI
/180.0;
697 copy_rvec(vec
, rot_vec
);
699 /* Precompute some variables: */
700 cosa
= cos(radangle
);
701 sina
= sin(radangle
);
703 dumxy
= rot_vec
[XX
]*rot_vec
[YY
]*OMcosa
;
704 dumxz
= rot_vec
[XX
]*rot_vec
[ZZ
]*OMcosa
;
705 dumyz
= rot_vec
[YY
]*rot_vec
[ZZ
]*OMcosa
;
707 /* Construct the rotation matrix for this rotation group: */
709 rotmat
[XX
][XX
] = cosa
+ rot_vec
[XX
]*rot_vec
[XX
]*OMcosa
;
710 rotmat
[YY
][XX
] = dumxy
+ rot_vec
[ZZ
]*sina
;
711 rotmat
[ZZ
][XX
] = dumxz
- rot_vec
[YY
]*sina
;
713 rotmat
[XX
][YY
] = dumxy
- rot_vec
[ZZ
]*sina
;
714 rotmat
[YY
][YY
] = cosa
+ rot_vec
[YY
]*rot_vec
[YY
]*OMcosa
;
715 rotmat
[ZZ
][YY
] = dumyz
+ rot_vec
[XX
]*sina
;
717 rotmat
[XX
][ZZ
] = dumxz
+ rot_vec
[YY
]*sina
;
718 rotmat
[YY
][ZZ
] = dumyz
- rot_vec
[XX
]*sina
;
719 rotmat
[ZZ
][ZZ
] = cosa
+ rot_vec
[ZZ
]*rot_vec
[ZZ
]*OMcosa
;
724 for (iii
= 0; iii
< 3; iii
++)
726 for (jjj
= 0; jjj
< 3; jjj
++)
728 fprintf(stderr
, " %10.8f ", rotmat
[iii
][jjj
]);
730 fprintf(stderr
, "\n");
736 /* Calculates torque on the rotation axis tau = position x force */
737 static gmx_inline real
torque(
738 rvec rotvec
, /* rotation vector; MUST be normalized! */
739 rvec force
, /* force */
740 rvec x
, /* position of atom on which the force acts */
741 rvec pivot
) /* pivot point of rotation axis */
746 /* Subtract offset */
747 rvec_sub(x
, pivot
, vectmp
);
749 /* position x force */
750 cprod(vectmp
, force
, tau
);
752 /* Return the part of the torque which is parallel to the rotation vector */
753 return iprod(tau
, rotvec
);
757 /* Right-aligned output of value with standard width */
758 static void print_aligned(FILE *fp
, char const *str
)
760 fprintf(fp
, "%12s", str
);
764 /* Right-aligned output of value with standard short width */
765 static void print_aligned_short(FILE *fp
, char const *str
)
767 fprintf(fp
, "%6s", str
);
771 static FILE *open_output_file(const char *fn
, int steps
, const char what
[])
776 fp
= gmx_ffopen(fn
, "w");
778 fprintf(fp
, "# Output of %s is written in intervals of %d time step%s.\n#\n",
779 what
, steps
, steps
> 1 ? "s" : "");
785 /* Open output file for slab center data. Call on master only */
786 static FILE *open_slab_out(const char *fn
, t_rot
*rot
)
793 if (rot
->enfrot
->Flags
& MD_APPENDFILES
)
795 fp
= gmx_fio_fopen(fn
, "a");
799 fp
= open_output_file(fn
, rot
->nstsout
, "gaussian weighted slab centers");
801 for (g
= 0; g
< rot
->ngrp
; g
++)
806 fprintf(fp
, "# Rotation group %d (%s), slab distance %f nm, %s.\n",
807 g
, erotg_names
[rotg
->eType
], rotg
->slab_dist
,
808 rotg
->bMassW
? "centers of mass" : "geometrical centers");
812 fprintf(fp
, "# Reference centers are listed first (t=-1).\n");
813 fprintf(fp
, "# The following columns have the syntax:\n");
815 print_aligned_short(fp
, "t");
816 print_aligned_short(fp
, "grp");
817 /* Print legend for the first two entries only ... */
818 for (i
= 0; i
< 2; i
++)
820 print_aligned_short(fp
, "slab");
821 print_aligned(fp
, "X center");
822 print_aligned(fp
, "Y center");
823 print_aligned(fp
, "Z center");
825 fprintf(fp
, " ...\n");
833 /* Adds 'buf' to 'str' */
834 static void add_to_string(char **str
, char *buf
)
839 len
= strlen(*str
) + strlen(buf
) + 1;
845 static void add_to_string_aligned(char **str
, char *buf
)
847 char buf_aligned
[STRLEN
];
849 sprintf(buf_aligned
, "%12s", buf
);
850 add_to_string(str
, buf_aligned
);
854 /* Open output file and print some general information about the rotation groups.
855 * Call on master only */
856 static FILE *open_rot_out(const char *fn
, t_rot
*rot
, const gmx_output_env_t
*oenv
)
861 const char **setname
;
862 char buf
[50], buf2
[75];
863 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
865 char *LegendStr
= nullptr;
868 if (rot
->enfrot
->Flags
& MD_APPENDFILES
)
870 fp
= gmx_fio_fopen(fn
, "a");
874 fp
= xvgropen(fn
, "Rotation angles and energy", "Time (ps)", "angles (degrees) and energies (kJ/mol)", oenv
);
875 fprintf(fp
, "# Output of enforced rotation data is written in intervals of %d time step%s.\n#\n", rot
->nstrout
, rot
->nstrout
> 1 ? "s" : "");
876 fprintf(fp
, "# The scalar tau is the torque (kJ/mol) in the direction of the rotation vector v.\n");
877 fprintf(fp
, "# To obtain the vectorial torque, multiply tau with the group's rot-vec.\n");
878 fprintf(fp
, "# For flexible groups, tau(t,n) from all slabs n have been summed in a single value tau(t) here.\n");
879 fprintf(fp
, "# The torques tau(t,n) are found in the rottorque.log (-rt) output file\n");
881 for (g
= 0; g
< rot
->ngrp
; g
++)
884 erg
= rotg
->enfrotgrp
;
885 bFlex
= ISFLEX(rotg
);
888 fprintf(fp
, "# ROTATION GROUP %d, potential type '%s':\n", g
, erotg_names
[rotg
->eType
]);
889 fprintf(fp
, "# rot-massw%d %s\n", g
, yesno_names
[rotg
->bMassW
]);
890 fprintf(fp
, "# rot-vec%d %12.5e %12.5e %12.5e\n", g
, rotg
->vec
[XX
], rotg
->vec
[YY
], rotg
->vec
[ZZ
]);
891 fprintf(fp
, "# rot-rate%d %12.5e degrees/ps\n", g
, rotg
->rate
);
892 fprintf(fp
, "# rot-k%d %12.5e kJ/(mol*nm^2)\n", g
, rotg
->k
);
893 if (rotg
->eType
== erotgISO
|| rotg
->eType
== erotgPM
|| rotg
->eType
== erotgRM
|| rotg
->eType
== erotgRM2
)
895 fprintf(fp
, "# rot-pivot%d %12.5e %12.5e %12.5e nm\n", g
, rotg
->pivot
[XX
], rotg
->pivot
[YY
], rotg
->pivot
[ZZ
]);
900 fprintf(fp
, "# rot-slab-distance%d %f nm\n", g
, rotg
->slab_dist
);
901 fprintf(fp
, "# rot-min-gaussian%d %12.5e\n", g
, rotg
->min_gaussian
);
904 /* Output the centers of the rotation groups for the pivot-free potentials */
905 if ((rotg
->eType
== erotgISOPF
) || (rotg
->eType
== erotgPMPF
) || (rotg
->eType
== erotgRMPF
) || (rotg
->eType
== erotgRM2PF
906 || (rotg
->eType
== erotgFLEXT
) || (rotg
->eType
== erotgFLEX2T
)) )
908 fprintf(fp
, "# ref. grp. %d center %12.5e %12.5e %12.5e\n", g
,
909 erg
->xc_ref_center
[XX
], erg
->xc_ref_center
[YY
], erg
->xc_ref_center
[ZZ
]);
911 fprintf(fp
, "# grp. %d init.center %12.5e %12.5e %12.5e\n", g
,
912 erg
->xc_center
[XX
], erg
->xc_center
[YY
], erg
->xc_center
[ZZ
]);
915 if ( (rotg
->eType
== erotgRM2
) || (rotg
->eType
== erotgFLEX2
) || (rotg
->eType
== erotgFLEX2T
) )
917 fprintf(fp
, "# rot-eps%d %12.5e nm^2\n", g
, rotg
->eps
);
919 if (erotgFitPOT
== rotg
->eFittype
)
922 fprintf(fp
, "# theta_fit%d is determined by first evaluating the potential for %d angles around theta_ref%d.\n",
923 g
, rotg
->PotAngle_nstep
, g
);
924 fprintf(fp
, "# The fit angle is the one with the smallest potential. It is given as the deviation\n");
925 fprintf(fp
, "# from the reference angle, i.e. if theta_ref=X and theta_fit=Y, then the angle with\n");
926 fprintf(fp
, "# minimal value of the potential is X+Y. Angular resolution is %g degrees.\n", rotg
->PotAngle_step
);
930 /* Print a nice legend */
933 sprintf(buf
, "# %6s", "time");
934 add_to_string_aligned(&LegendStr
, buf
);
937 snew(setname
, 4*rot
->ngrp
);
939 for (g
= 0; g
< rot
->ngrp
; g
++)
941 sprintf(buf
, "theta_ref%d", g
);
942 add_to_string_aligned(&LegendStr
, buf
);
944 sprintf(buf2
, "%s (degrees)", buf
);
945 setname
[nsets
] = gmx_strdup(buf2
);
948 for (g
= 0; g
< rot
->ngrp
; g
++)
951 bFlex
= ISFLEX(rotg
);
953 /* For flexible axis rotation we use RMSD fitting to determine the
954 * actual angle of the rotation group */
955 if (bFlex
|| erotgFitPOT
== rotg
->eFittype
)
957 sprintf(buf
, "theta_fit%d", g
);
961 sprintf(buf
, "theta_av%d", g
);
963 add_to_string_aligned(&LegendStr
, buf
);
964 sprintf(buf2
, "%s (degrees)", buf
);
965 setname
[nsets
] = gmx_strdup(buf2
);
968 sprintf(buf
, "tau%d", g
);
969 add_to_string_aligned(&LegendStr
, buf
);
970 sprintf(buf2
, "%s (kJ/mol)", buf
);
971 setname
[nsets
] = gmx_strdup(buf2
);
974 sprintf(buf
, "energy%d", g
);
975 add_to_string_aligned(&LegendStr
, buf
);
976 sprintf(buf2
, "%s (kJ/mol)", buf
);
977 setname
[nsets
] = gmx_strdup(buf2
);
984 xvgr_legend(fp
, nsets
, setname
, oenv
);
988 fprintf(fp
, "#\n# Legend for the following data columns:\n");
989 fprintf(fp
, "%s\n", LegendStr
);
999 /* Call on master only */
1000 static FILE *open_angles_out(const char *fn
, t_rot
*rot
)
1005 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1009 if (rot
->enfrot
->Flags
& MD_APPENDFILES
)
1011 fp
= gmx_fio_fopen(fn
, "a");
1015 /* Open output file and write some information about it's structure: */
1016 fp
= open_output_file(fn
, rot
->nstsout
, "rotation group angles");
1017 fprintf(fp
, "# All angles given in degrees, time in ps.\n");
1018 for (g
= 0; g
< rot
->ngrp
; g
++)
1020 rotg
= &rot
->grp
[g
];
1021 erg
= rotg
->enfrotgrp
;
1023 /* Output for this group happens only if potential type is flexible or
1024 * if fit type is potential! */
1025 if (ISFLEX(rotg
) || (erotgFitPOT
== rotg
->eFittype
) )
1029 sprintf(buf
, " slab distance %f nm, ", rotg
->slab_dist
);
1036 fprintf(fp
, "#\n# ROTATION GROUP %d '%s',%s fit type '%s'.\n",
1037 g
, erotg_names
[rotg
->eType
], buf
, erotg_fitnames
[rotg
->eFittype
]);
1039 /* Special type of fitting using the potential minimum. This is
1040 * done for the whole group only, not for the individual slabs. */
1041 if (erotgFitPOT
== rotg
->eFittype
)
1043 fprintf(fp
, "# To obtain theta_fit%d, the potential is evaluated for %d angles around theta_ref%d\n", g
, rotg
->PotAngle_nstep
, g
);
1044 fprintf(fp
, "# The fit angle in the rotation standard outfile is the one with minimal energy E(theta_fit) [kJ/mol].\n");
1048 fprintf(fp
, "# Legend for the group %d data columns:\n", g
);
1050 print_aligned_short(fp
, "time");
1051 print_aligned_short(fp
, "grp");
1052 print_aligned(fp
, "theta_ref");
1054 if (erotgFitPOT
== rotg
->eFittype
)
1056 /* Output the set of angles around the reference angle */
1057 for (i
= 0; i
< rotg
->PotAngle_nstep
; i
++)
1059 sprintf(buf
, "E(%g)", erg
->PotAngleFit
->degangle
[i
]);
1060 print_aligned(fp
, buf
);
1065 /* Output fit angle for each slab */
1066 print_aligned_short(fp
, "slab");
1067 print_aligned_short(fp
, "atoms");
1068 print_aligned(fp
, "theta_fit");
1069 print_aligned_short(fp
, "slab");
1070 print_aligned_short(fp
, "atoms");
1071 print_aligned(fp
, "theta_fit");
1072 fprintf(fp
, " ...");
1084 /* Open torque output file and write some information about it's structure.
1085 * Call on master only */
1086 static FILE *open_torque_out(const char *fn
, t_rot
*rot
)
1093 if (rot
->enfrot
->Flags
& MD_APPENDFILES
)
1095 fp
= gmx_fio_fopen(fn
, "a");
1099 fp
= open_output_file(fn
, rot
->nstsout
, "torques");
1101 for (g
= 0; g
< rot
->ngrp
; g
++)
1103 rotg
= &rot
->grp
[g
];
1106 fprintf(fp
, "# Rotation group %d (%s), slab distance %f nm.\n", g
, erotg_names
[rotg
->eType
], rotg
->slab_dist
);
1107 fprintf(fp
, "# The scalar tau is the torque (kJ/mol) in the direction of the rotation vector.\n");
1108 fprintf(fp
, "# To obtain the vectorial torque, multiply tau with\n");
1109 fprintf(fp
, "# rot-vec%d %10.3e %10.3e %10.3e\n", g
, rotg
->vec
[XX
], rotg
->vec
[YY
], rotg
->vec
[ZZ
]);
1113 fprintf(fp
, "# Legend for the following data columns: (tau=torque for that slab):\n");
1115 print_aligned_short(fp
, "t");
1116 print_aligned_short(fp
, "grp");
1117 print_aligned_short(fp
, "slab");
1118 print_aligned(fp
, "tau");
1119 print_aligned_short(fp
, "slab");
1120 print_aligned(fp
, "tau");
1121 fprintf(fp
, " ...\n");
1129 static void swap_val(double* vec
, int i
, int j
)
1131 double tmp
= vec
[j
];
1139 static void swap_col(double **mat
, int i
, int j
)
1141 double tmp
[3] = {mat
[0][j
], mat
[1][j
], mat
[2][j
]};
1144 mat
[0][j
] = mat
[0][i
];
1145 mat
[1][j
] = mat
[1][i
];
1146 mat
[2][j
] = mat
[2][i
];
1154 /* Eigenvectors are stored in columns of eigen_vec */
1155 static void diagonalize_symmetric(
1163 jacobi(matrix
, 3, eigenval
, eigen_vec
, &n_rot
);
1165 /* sort in ascending order */
1166 if (eigenval
[0] > eigenval
[1])
1168 swap_val(eigenval
, 0, 1);
1169 swap_col(eigen_vec
, 0, 1);
1171 if (eigenval
[1] > eigenval
[2])
1173 swap_val(eigenval
, 1, 2);
1174 swap_col(eigen_vec
, 1, 2);
1176 if (eigenval
[0] > eigenval
[1])
1178 swap_val(eigenval
, 0, 1);
1179 swap_col(eigen_vec
, 0, 1);
1184 static void align_with_z(
1185 rvec
* s
, /* Structure to align */
1190 rvec zet
= {0.0, 0.0, 1.0};
1191 rvec rot_axis
= {0.0, 0.0, 0.0};
1192 rvec
*rotated_str
= nullptr;
1198 snew(rotated_str
, natoms
);
1200 /* Normalize the axis */
1201 ooanorm
= 1.0/norm(axis
);
1202 svmul(ooanorm
, axis
, axis
);
1204 /* Calculate the angle for the fitting procedure */
1205 cprod(axis
, zet
, rot_axis
);
1206 angle
= acos(axis
[2]);
1212 /* Calculate the rotation matrix */
1213 calc_rotmat(rot_axis
, angle
*180.0/M_PI
, rotmat
);
1215 /* Apply the rotation matrix to s */
1216 for (i
= 0; i
< natoms
; i
++)
1218 for (j
= 0; j
< 3; j
++)
1220 for (k
= 0; k
< 3; k
++)
1222 rotated_str
[i
][j
] += rotmat
[j
][k
]*s
[i
][k
];
1227 /* Rewrite the rotated structure to s */
1228 for (i
= 0; i
< natoms
; i
++)
1230 for (j
= 0; j
< 3; j
++)
1232 s
[i
][j
] = rotated_str
[i
][j
];
1240 static void calc_correl_matrix(rvec
* Xstr
, rvec
* Ystr
, double** Rmat
, int natoms
)
1245 for (i
= 0; i
< 3; i
++)
1247 for (j
= 0; j
< 3; j
++)
1253 for (i
= 0; i
< 3; i
++)
1255 for (j
= 0; j
< 3; j
++)
1257 for (k
= 0; k
< natoms
; k
++)
1259 Rmat
[i
][j
] += Ystr
[k
][i
] * Xstr
[k
][j
];
1266 static void weigh_coords(rvec
* str
, real
* weight
, int natoms
)
1271 for (i
= 0; i
< natoms
; i
++)
1273 for (j
= 0; j
< 3; j
++)
1275 str
[i
][j
] *= std::sqrt(weight
[i
]);
1281 static real
opt_angle_analytic(
1291 rvec
*ref_s_1
= nullptr;
1292 rvec
*act_s_1
= nullptr;
1294 double **Rmat
, **RtR
, **eigvec
;
1296 double V
[3][3], WS
[3][3];
1297 double rot_matrix
[3][3];
1301 /* Do not change the original coordinates */
1302 snew(ref_s_1
, natoms
);
1303 snew(act_s_1
, natoms
);
1304 for (i
= 0; i
< natoms
; i
++)
1306 copy_rvec(ref_s
[i
], ref_s_1
[i
]);
1307 copy_rvec(act_s
[i
], act_s_1
[i
]);
1310 /* Translate the structures to the origin */
1311 shift
[XX
] = -ref_com
[XX
];
1312 shift
[YY
] = -ref_com
[YY
];
1313 shift
[ZZ
] = -ref_com
[ZZ
];
1314 translate_x(ref_s_1
, natoms
, shift
);
1316 shift
[XX
] = -act_com
[XX
];
1317 shift
[YY
] = -act_com
[YY
];
1318 shift
[ZZ
] = -act_com
[ZZ
];
1319 translate_x(act_s_1
, natoms
, shift
);
1321 /* Align rotation axis with z */
1322 align_with_z(ref_s_1
, natoms
, axis
);
1323 align_with_z(act_s_1
, natoms
, axis
);
1325 /* Correlation matrix */
1326 Rmat
= allocate_square_matrix(3);
1328 for (i
= 0; i
< natoms
; i
++)
1330 ref_s_1
[i
][2] = 0.0;
1331 act_s_1
[i
][2] = 0.0;
1334 /* Weight positions with sqrt(weight) */
1335 if (nullptr != weight
)
1337 weigh_coords(ref_s_1
, weight
, natoms
);
1338 weigh_coords(act_s_1
, weight
, natoms
);
1341 /* Calculate correlation matrices R=YXt (X=ref_s; Y=act_s) */
1342 calc_correl_matrix(ref_s_1
, act_s_1
, Rmat
, natoms
);
1345 RtR
= allocate_square_matrix(3);
1346 for (i
= 0; i
< 3; i
++)
1348 for (j
= 0; j
< 3; j
++)
1350 for (k
= 0; k
< 3; k
++)
1352 RtR
[i
][j
] += Rmat
[k
][i
] * Rmat
[k
][j
];
1356 /* Diagonalize RtR */
1358 for (i
= 0; i
< 3; i
++)
1363 diagonalize_symmetric(RtR
, eigvec
, eigval
);
1364 swap_col(eigvec
, 0, 1);
1365 swap_col(eigvec
, 1, 2);
1366 swap_val(eigval
, 0, 1);
1367 swap_val(eigval
, 1, 2);
1370 for (i
= 0; i
< 3; i
++)
1372 for (j
= 0; j
< 3; j
++)
1379 for (i
= 0; i
< 2; i
++)
1381 for (j
= 0; j
< 2; j
++)
1383 WS
[i
][j
] = eigvec
[i
][j
] / std::sqrt(eigval
[j
]);
1387 for (i
= 0; i
< 3; i
++)
1389 for (j
= 0; j
< 3; j
++)
1391 for (k
= 0; k
< 3; k
++)
1393 V
[i
][j
] += Rmat
[i
][k
]*WS
[k
][j
];
1397 free_square_matrix(Rmat
, 3);
1399 /* Calculate optimal rotation matrix */
1400 for (i
= 0; i
< 3; i
++)
1402 for (j
= 0; j
< 3; j
++)
1404 rot_matrix
[i
][j
] = 0.0;
1408 for (i
= 0; i
< 3; i
++)
1410 for (j
= 0; j
< 3; j
++)
1412 for (k
= 0; k
< 3; k
++)
1414 rot_matrix
[i
][j
] += eigvec
[i
][k
]*V
[j
][k
];
1418 rot_matrix
[2][2] = 1.0;
1420 /* In some cases abs(rot_matrix[0][0]) can be slighly larger
1421 * than unity due to numerical inacurracies. To be able to calculate
1422 * the acos function, we put these values back in range. */
1423 if (rot_matrix
[0][0] > 1.0)
1425 rot_matrix
[0][0] = 1.0;
1427 else if (rot_matrix
[0][0] < -1.0)
1429 rot_matrix
[0][0] = -1.0;
1432 /* Determine the optimal rotation angle: */
1433 opt_angle
= (-1.0)*acos(rot_matrix
[0][0])*180.0/M_PI
;
1434 if (rot_matrix
[0][1] < 0.0)
1436 opt_angle
= (-1.0)*opt_angle
;
1439 /* Give back some memory */
1440 free_square_matrix(RtR
, 3);
1443 for (i
= 0; i
< 3; i
++)
1449 return (real
) opt_angle
;
1453 /* Determine angle of the group by RMSD fit to the reference */
1454 /* Not parallelized, call this routine only on the master */
1455 static real
flex_fit_angle(t_rotgrp
*rotg
)
1458 rvec
*fitcoords
= nullptr;
1459 rvec center
; /* Center of positions passed to the fit routine */
1460 real fitangle
; /* Angle of the rotation group derived by fitting */
1463 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1466 erg
= rotg
->enfrotgrp
;
1468 /* Get the center of the rotation group.
1469 * Note, again, erg->xc has been sorted in do_flexible */
1470 get_center(erg
->xc
, erg
->mc_sorted
, rotg
->nat
, center
);
1472 /* === Determine the optimal fit angle for the rotation group === */
1473 if (rotg
->eFittype
== erotgFitNORM
)
1475 /* Normalize every position to it's reference length */
1476 for (i
= 0; i
< rotg
->nat
; i
++)
1478 /* Put the center of the positions into the origin */
1479 rvec_sub(erg
->xc
[i
], center
, coord
);
1480 /* Determine the scaling factor for the length: */
1481 scal
= erg
->xc_ref_length
[erg
->xc_sortind
[i
]] / norm(coord
);
1482 /* Get position, multiply with the scaling factor and save */
1483 svmul(scal
, coord
, erg
->xc_norm
[i
]);
1485 fitcoords
= erg
->xc_norm
;
1489 fitcoords
= erg
->xc
;
1491 /* From the point of view of the current positions, the reference has rotated
1492 * backwards. Since we output the angle relative to the fixed reference,
1493 * we need the minus sign. */
1494 fitangle
= -opt_angle_analytic(erg
->xc_ref_sorted
, fitcoords
, erg
->mc_sorted
,
1495 rotg
->nat
, erg
->xc_ref_center
, center
, rotg
->vec
);
1501 /* Determine actual angle of each slab by RMSD fit to the reference */
1502 /* Not parallelized, call this routine only on the master */
1503 static void flex_fit_angle_perslab(
1510 int i
, l
, n
, islab
, ind
;
1512 rvec act_center
; /* Center of actual positions that are passed to the fit routine */
1513 rvec ref_center
; /* Same for the reference positions */
1514 real fitangle
; /* Angle of a slab derived from an RMSD fit to
1515 * the reference structure at t=0 */
1517 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1518 real OOm_av
; /* 1/average_mass of a rotation group atom */
1519 real m_rel
; /* Relative mass of a rotation group atom */
1522 erg
= rotg
->enfrotgrp
;
1524 /* Average mass of a rotation group atom: */
1525 OOm_av
= erg
->invmass
*rotg
->nat
;
1527 /**********************************/
1528 /* First collect the data we need */
1529 /**********************************/
1531 /* Collect the data for the individual slabs */
1532 for (n
= erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1534 islab
= n
- erg
->slab_first
; /* slab index */
1535 sd
= &(rotg
->enfrotgrp
->slab_data
[islab
]);
1536 sd
->nat
= erg
->lastatom
[islab
]-erg
->firstatom
[islab
]+1;
1539 /* Loop over the relevant atoms in the slab */
1540 for (l
= erg
->firstatom
[islab
]; l
<= erg
->lastatom
[islab
]; l
++)
1542 /* Current position of this atom: x[ii][XX/YY/ZZ] */
1543 copy_rvec(erg
->xc
[l
], curr_x
);
1545 /* The (unrotated) reference position of this atom is copied to ref_x.
1546 * Beware, the xc coords have been sorted in do_flexible */
1547 copy_rvec(erg
->xc_ref_sorted
[l
], ref_x
);
1549 /* Save data for doing angular RMSD fit later */
1550 /* Save the current atom position */
1551 copy_rvec(curr_x
, sd
->x
[ind
]);
1552 /* Save the corresponding reference position */
1553 copy_rvec(ref_x
, sd
->ref
[ind
]);
1555 /* Maybe also mass-weighting was requested. If yes, additionally
1556 * multiply the weights with the relative mass of the atom. If not,
1557 * multiply with unity. */
1558 m_rel
= erg
->mc_sorted
[l
]*OOm_av
;
1560 /* Save the weight for this atom in this slab */
1561 sd
->weight
[ind
] = gaussian_weight(curr_x
, rotg
, n
) * m_rel
;
1563 /* Next atom in this slab */
1568 /******************************/
1569 /* Now do the fit calculation */
1570 /******************************/
1572 fprintf(fp
, "%12.3e%6d%12.3f", t
, g
, degangle
);
1574 /* === Now do RMSD fitting for each slab === */
1575 /* We require at least SLAB_MIN_ATOMS in a slab, such that the fit makes sense. */
1576 #define SLAB_MIN_ATOMS 4
1578 for (n
= erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1580 islab
= n
- erg
->slab_first
; /* slab index */
1581 sd
= &(rotg
->enfrotgrp
->slab_data
[islab
]);
1582 if (sd
->nat
>= SLAB_MIN_ATOMS
)
1584 /* Get the center of the slabs reference and current positions */
1585 get_center(sd
->ref
, sd
->weight
, sd
->nat
, ref_center
);
1586 get_center(sd
->x
, sd
->weight
, sd
->nat
, act_center
);
1587 if (rotg
->eFittype
== erotgFitNORM
)
1589 /* Normalize every position to it's reference length
1590 * prior to performing the fit */
1591 for (i
= 0; i
< sd
->nat
; i
++) /* Center */
1593 rvec_dec(sd
->ref
[i
], ref_center
);
1594 rvec_dec(sd
->x
[i
], act_center
);
1595 /* Normalize x_i such that it gets the same length as ref_i */
1596 svmul( norm(sd
->ref
[i
])/norm(sd
->x
[i
]), sd
->x
[i
], sd
->x
[i
] );
1598 /* We already subtracted the centers */
1599 clear_rvec(ref_center
);
1600 clear_rvec(act_center
);
1602 fitangle
= -opt_angle_analytic(sd
->ref
, sd
->x
, sd
->weight
, sd
->nat
,
1603 ref_center
, act_center
, rotg
->vec
);
1604 fprintf(fp
, "%6d%6d%12.3f", n
, sd
->nat
, fitangle
);
1609 #undef SLAB_MIN_ATOMS
1613 /* Shift x with is */
1614 static gmx_inline
void shift_single_coord(matrix box
, rvec x
, const ivec is
)
1625 x
[XX
] += tx
*box
[XX
][XX
]+ty
*box
[YY
][XX
]+tz
*box
[ZZ
][XX
];
1626 x
[YY
] += ty
*box
[YY
][YY
]+tz
*box
[ZZ
][YY
];
1627 x
[ZZ
] += tz
*box
[ZZ
][ZZ
];
1631 x
[XX
] += tx
*box
[XX
][XX
];
1632 x
[YY
] += ty
*box
[YY
][YY
];
1633 x
[ZZ
] += tz
*box
[ZZ
][ZZ
];
1638 /* Determine the 'home' slab of this atom which is the
1639 * slab with the highest Gaussian weight of all */
1640 #define round(a) (int)(a+0.5)
1641 static gmx_inline
int get_homeslab(
1642 rvec curr_x
, /* The position for which the home slab shall be determined */
1643 rvec rotvec
, /* The rotation vector */
1644 real slabdist
) /* The slab distance */
1649 /* The distance of the atom to the coordinate center (where the
1650 * slab with index 0) is */
1651 dist
= iprod(rotvec
, curr_x
);
1653 return round(dist
/ slabdist
);
1657 /* For a local atom determine the relevant slabs, i.e. slabs in
1658 * which the gaussian is larger than min_gaussian
1660 static int get_single_atom_gaussians(
1667 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1670 erg
= rotg
->enfrotgrp
;
1672 /* Determine the 'home' slab of this atom: */
1673 homeslab
= get_homeslab(curr_x
, rotg
->vec
, rotg
->slab_dist
);
1675 /* First determine the weight in the atoms home slab: */
1676 g
= gaussian_weight(curr_x
, rotg
, homeslab
);
1678 erg
->gn_atom
[count
] = g
;
1679 erg
->gn_slabind
[count
] = homeslab
;
1683 /* Determine the max slab */
1685 while (g
> rotg
->min_gaussian
)
1688 g
= gaussian_weight(curr_x
, rotg
, slab
);
1689 erg
->gn_slabind
[count
] = slab
;
1690 erg
->gn_atom
[count
] = g
;
1695 /* Determine the min slab */
1700 g
= gaussian_weight(curr_x
, rotg
, slab
);
1701 erg
->gn_slabind
[count
] = slab
;
1702 erg
->gn_atom
[count
] = g
;
1705 while (g
> rotg
->min_gaussian
);
1712 static void flex2_precalc_inner_sum(t_rotgrp
*rotg
)
1715 rvec xi
; /* positions in the i-sum */
1716 rvec xcn
, ycn
; /* the current and the reference slab centers */
1719 rvec rin
; /* Helper variables */
1722 real OOpsii
, OOpsiistar
;
1723 real sin_rin
; /* s_ii.r_ii */
1724 rvec s_in
, tmpvec
, tmpvec2
;
1725 real mi
, wi
; /* Mass-weighting of the positions */
1727 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1730 erg
= rotg
->enfrotgrp
;
1731 N_M
= rotg
->nat
* erg
->invmass
;
1733 /* Loop over all slabs that contain something */
1734 for (n
= erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1736 islab
= n
- erg
->slab_first
; /* slab index */
1738 /* The current center of this slab is saved in xcn: */
1739 copy_rvec(erg
->slab_center
[islab
], xcn
);
1740 /* ... and the reference center in ycn: */
1741 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
1743 /*** D. Calculate the whole inner sum used for second and third sum */
1744 /* For slab n, we need to loop over all atoms i again. Since we sorted
1745 * the atoms with respect to the rotation vector, we know that it is sufficient
1746 * to calculate from firstatom to lastatom only. All other contributions will
1748 clear_rvec(innersumvec
);
1749 for (i
= erg
->firstatom
[islab
]; i
<= erg
->lastatom
[islab
]; i
++)
1751 /* Coordinate xi of this atom */
1752 copy_rvec(erg
->xc
[i
], xi
);
1755 gaussian_xi
= gaussian_weight(xi
, rotg
, n
);
1756 mi
= erg
->mc_sorted
[i
]; /* need the sorted mass here */
1760 copy_rvec(erg
->xc_ref_sorted
[i
], yi0
); /* Reference position yi0 */
1761 rvec_sub(yi0
, ycn
, tmpvec2
); /* tmpvec2 = yi0 - ycn */
1762 mvmul(erg
->rotmat
, tmpvec2
, rin
); /* rin = Omega.(yi0 - ycn) */
1764 /* Calculate psi_i* and sin */
1765 rvec_sub(xi
, xcn
, tmpvec2
); /* tmpvec2 = xi - xcn */
1766 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x (xi - xcn) */
1767 OOpsiistar
= norm2(tmpvec
)+rotg
->eps
; /* OOpsii* = 1/psii* = |v x (xi-xcn)|^2 + eps */
1768 OOpsii
= norm(tmpvec
); /* OOpsii = 1 / psii = |v x (xi - xcn)| */
1770 /* * v x (xi - xcn) */
1771 unitv(tmpvec
, s_in
); /* sin = ---------------- */
1772 /* |v x (xi - xcn)| */
1774 sin_rin
= iprod(s_in
, rin
); /* sin_rin = sin . rin */
1776 /* Now the whole sum */
1777 fac
= OOpsii
/OOpsiistar
;
1778 svmul(fac
, rin
, tmpvec
);
1779 fac2
= fac
*fac
*OOpsii
;
1780 svmul(fac2
*sin_rin
, s_in
, tmpvec2
);
1781 rvec_dec(tmpvec
, tmpvec2
);
1783 svmul(wi
*gaussian_xi
*sin_rin
, tmpvec
, tmpvec2
);
1785 rvec_inc(innersumvec
, tmpvec2
);
1786 } /* now we have the inner sum, used both for sum2 and sum3 */
1788 /* Save it to be used in do_flex2_lowlevel */
1789 copy_rvec(innersumvec
, erg
->slab_innersumvec
[islab
]);
1790 } /* END of loop over slabs */
1794 static void flex_precalc_inner_sum(t_rotgrp
*rotg
)
1797 rvec xi
; /* position */
1798 rvec xcn
, ycn
; /* the current and the reference slab centers */
1799 rvec qin
, rin
; /* q_i^n and r_i^n */
1802 rvec innersumvec
; /* Inner part of sum_n2 */
1803 real gaussian_xi
; /* Gaussian weight gn(xi) */
1804 real mi
, wi
; /* Mass-weighting of the positions */
1807 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1810 erg
= rotg
->enfrotgrp
;
1811 N_M
= rotg
->nat
* erg
->invmass
;
1813 /* Loop over all slabs that contain something */
1814 for (n
= erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1816 islab
= n
- erg
->slab_first
; /* slab index */
1818 /* The current center of this slab is saved in xcn: */
1819 copy_rvec(erg
->slab_center
[islab
], xcn
);
1820 /* ... and the reference center in ycn: */
1821 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
1823 /* For slab n, we need to loop over all atoms i again. Since we sorted
1824 * the atoms with respect to the rotation vector, we know that it is sufficient
1825 * to calculate from firstatom to lastatom only. All other contributions will
1827 clear_rvec(innersumvec
);
1828 for (i
= erg
->firstatom
[islab
]; i
<= erg
->lastatom
[islab
]; i
++)
1830 /* Coordinate xi of this atom */
1831 copy_rvec(erg
->xc
[i
], xi
);
1834 gaussian_xi
= gaussian_weight(xi
, rotg
, n
);
1835 mi
= erg
->mc_sorted
[i
]; /* need the sorted mass here */
1838 /* Calculate rin and qin */
1839 rvec_sub(erg
->xc_ref_sorted
[i
], ycn
, tmpvec
); /* tmpvec = yi0-ycn */
1840 mvmul(erg
->rotmat
, tmpvec
, rin
); /* rin = Omega.(yi0 - ycn) */
1841 cprod(rotg
->vec
, rin
, tmpvec
); /* tmpvec = v x Omega*(yi0-ycn) */
1843 /* * v x Omega*(yi0-ycn) */
1844 unitv(tmpvec
, qin
); /* qin = --------------------- */
1845 /* |v x Omega*(yi0-ycn)| */
1848 rvec_sub(xi
, xcn
, tmpvec
); /* tmpvec = xi-xcn */
1849 bin
= iprod(qin
, tmpvec
); /* bin = qin*(xi-xcn) */
1851 svmul(wi
*gaussian_xi
*bin
, qin
, tmpvec
);
1853 /* Add this contribution to the inner sum: */
1854 rvec_add(innersumvec
, tmpvec
, innersumvec
);
1855 } /* now we have the inner sum vector S^n for this slab */
1856 /* Save it to be used in do_flex_lowlevel */
1857 copy_rvec(innersumvec
, erg
->slab_innersumvec
[islab
]);
1862 static real
do_flex2_lowlevel(
1864 real sigma
, /* The Gaussian width sigma */
1866 gmx_bool bOutstepRot
,
1867 gmx_bool bOutstepSlab
,
1870 int count
, ic
, ii
, j
, m
, n
, islab
, iigrp
, ifit
;
1871 rvec xj
; /* position in the i-sum */
1872 rvec yj0
; /* the reference position in the j-sum */
1873 rvec xcn
, ycn
; /* the current and the reference slab centers */
1874 real V
; /* This node's part of the rotation pot. energy */
1875 real gaussian_xj
; /* Gaussian weight */
1878 real numerator
, fit_numerator
;
1879 rvec rjn
, fit_rjn
; /* Helper variables */
1882 real OOpsij
, OOpsijstar
;
1883 real OOsigma2
; /* 1/(sigma^2) */
1886 rvec sjn
, tmpvec
, tmpvec2
, yj0_ycn
;
1887 rvec sum1vec_part
, sum1vec
, sum2vec_part
, sum2vec
, sum3vec
, sum4vec
, innersumvec
;
1889 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1890 real mj
, wj
; /* Mass-weighting of the positions */
1892 real Wjn
; /* g_n(x_j) m_j / Mjn */
1893 gmx_bool bCalcPotFit
;
1895 /* To calculate the torque per slab */
1896 rvec slab_force
; /* Single force from slab n on one atom */
1897 rvec slab_sum1vec_part
;
1898 real slab_sum3part
, slab_sum4part
;
1899 rvec slab_sum1vec
, slab_sum2vec
, slab_sum3vec
, slab_sum4vec
;
1902 erg
= rotg
->enfrotgrp
;
1904 /* Pre-calculate the inner sums, so that we do not have to calculate
1905 * them again for every atom */
1906 flex2_precalc_inner_sum(rotg
);
1908 bCalcPotFit
= (bOutstepRot
|| bOutstepSlab
) && (erotgFitPOT
== rotg
->eFittype
);
1910 /********************************************************/
1911 /* Main loop over all local atoms of the rotation group */
1912 /********************************************************/
1913 N_M
= rotg
->nat
* erg
->invmass
;
1915 OOsigma2
= 1.0 / (sigma
*sigma
);
1916 for (j
= 0; j
< erg
->nat_loc
; j
++)
1918 /* Local index of a rotation group atom */
1919 ii
= erg
->ind_loc
[j
];
1920 /* Position of this atom in the collective array */
1921 iigrp
= erg
->xc_ref_ind
[j
];
1922 /* Mass-weighting */
1923 mj
= erg
->mc
[iigrp
]; /* need the unsorted mass here */
1926 /* Current position of this atom: x[ii][XX/YY/ZZ]
1927 * Note that erg->xc_center contains the center of mass in case the flex2-t
1928 * potential was chosen. For the flex2 potential erg->xc_center must be
1930 rvec_sub(x
[ii
], erg
->xc_center
, xj
);
1932 /* Shift this atom such that it is near its reference */
1933 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
1935 /* Determine the slabs to loop over, i.e. the ones with contributions
1936 * larger than min_gaussian */
1937 count
= get_single_atom_gaussians(xj
, rotg
);
1939 clear_rvec(sum1vec_part
);
1940 clear_rvec(sum2vec_part
);
1943 /* Loop over the relevant slabs for this atom */
1944 for (ic
= 0; ic
< count
; ic
++)
1946 n
= erg
->gn_slabind
[ic
];
1948 /* Get the precomputed Gaussian value of curr_slab for curr_x */
1949 gaussian_xj
= erg
->gn_atom
[ic
];
1951 islab
= n
- erg
->slab_first
; /* slab index */
1953 /* The (unrotated) reference position of this atom is copied to yj0: */
1954 copy_rvec(rotg
->x_ref
[iigrp
], yj0
);
1956 beta
= calc_beta(xj
, rotg
, n
);
1958 /* The current center of this slab is saved in xcn: */
1959 copy_rvec(erg
->slab_center
[islab
], xcn
);
1960 /* ... and the reference center in ycn: */
1961 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
1963 rvec_sub(yj0
, ycn
, yj0_ycn
); /* yj0_ycn = yj0 - ycn */
1966 mvmul(erg
->rotmat
, yj0_ycn
, rjn
); /* rjn = Omega.(yj0 - ycn) */
1968 /* Subtract the slab center from xj */
1969 rvec_sub(xj
, xcn
, tmpvec2
); /* tmpvec2 = xj - xcn */
1971 /* In rare cases, when an atom position coincides with a slab center
1972 * (tmpvec2 == 0) we cannot compute the vector product for sjn.
1973 * However, since the atom is located directly on the pivot, this
1974 * slab's contribution to the force on that atom will be zero
1975 * anyway. Therefore, we directly move on to the next slab. */
1976 if (gmx_numzero(norm(tmpvec2
))) /* 0 == norm(xj - xcn) */
1982 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x (xj - xcn) */
1984 OOpsijstar
= norm2(tmpvec
)+rotg
->eps
; /* OOpsij* = 1/psij* = |v x (xj-xcn)|^2 + eps */
1986 numerator
= gmx::square(iprod(tmpvec
, rjn
));
1988 /*********************************/
1989 /* Add to the rotation potential */
1990 /*********************************/
1991 V
+= 0.5*rotg
->k
*wj
*gaussian_xj
*numerator
/OOpsijstar
;
1993 /* If requested, also calculate the potential for a set of angles
1994 * near the current reference angle */
1997 for (ifit
= 0; ifit
< rotg
->PotAngle_nstep
; ifit
++)
1999 mvmul(erg
->PotAngleFit
->rotmat
[ifit
], yj0_ycn
, fit_rjn
);
2000 fit_numerator
= gmx::square(iprod(tmpvec
, fit_rjn
));
2001 erg
->PotAngleFit
->V
[ifit
] += 0.5*rotg
->k
*wj
*gaussian_xj
*fit_numerator
/OOpsijstar
;
2005 /*************************************/
2006 /* Now calculate the force on atom j */
2007 /*************************************/
2009 OOpsij
= norm(tmpvec
); /* OOpsij = 1 / psij = |v x (xj - xcn)| */
2011 /* * v x (xj - xcn) */
2012 unitv(tmpvec
, sjn
); /* sjn = ---------------- */
2013 /* |v x (xj - xcn)| */
2015 sjn_rjn
= iprod(sjn
, rjn
); /* sjn_rjn = sjn . rjn */
2018 /*** A. Calculate the first of the four sum terms: ****************/
2019 fac
= OOpsij
/OOpsijstar
;
2020 svmul(fac
, rjn
, tmpvec
);
2021 fac2
= fac
*fac
*OOpsij
;
2022 svmul(fac2
*sjn_rjn
, sjn
, tmpvec2
);
2023 rvec_dec(tmpvec
, tmpvec2
);
2024 fac2
= wj
*gaussian_xj
; /* also needed for sum4 */
2025 svmul(fac2
*sjn_rjn
, tmpvec
, slab_sum1vec_part
);
2026 /********************/
2027 /*** Add to sum1: ***/
2028 /********************/
2029 rvec_inc(sum1vec_part
, slab_sum1vec_part
); /* sum1 still needs to vector multiplied with v */
2031 /*** B. Calculate the forth of the four sum terms: ****************/
2032 betasigpsi
= beta
*OOsigma2
*OOpsij
; /* this is also needed for sum3 */
2033 /********************/
2034 /*** Add to sum4: ***/
2035 /********************/
2036 slab_sum4part
= fac2
*betasigpsi
*fac
*sjn_rjn
*sjn_rjn
; /* Note that fac is still valid from above */
2037 sum4
+= slab_sum4part
;
2039 /*** C. Calculate Wjn for second and third sum */
2040 /* Note that we can safely divide by slab_weights since we check in
2041 * get_slab_centers that it is non-zero. */
2042 Wjn
= gaussian_xj
*mj
/erg
->slab_weights
[islab
];
2044 /* We already have precalculated the inner sum for slab n */
2045 copy_rvec(erg
->slab_innersumvec
[islab
], innersumvec
);
2047 /* Weigh the inner sum vector with Wjn */
2048 svmul(Wjn
, innersumvec
, innersumvec
);
2050 /*** E. Calculate the second of the four sum terms: */
2051 /********************/
2052 /*** Add to sum2: ***/
2053 /********************/
2054 rvec_inc(sum2vec_part
, innersumvec
); /* sum2 still needs to be vector crossproduct'ed with v */
2056 /*** F. Calculate the third of the four sum terms: */
2057 slab_sum3part
= betasigpsi
* iprod(sjn
, innersumvec
);
2058 sum3
+= slab_sum3part
; /* still needs to be multiplied with v */
2060 /*** G. Calculate the torque on the local slab's axis: */
2064 cprod(slab_sum1vec_part
, rotg
->vec
, slab_sum1vec
);
2066 cprod(innersumvec
, rotg
->vec
, slab_sum2vec
);
2068 svmul(slab_sum3part
, rotg
->vec
, slab_sum3vec
);
2070 svmul(slab_sum4part
, rotg
->vec
, slab_sum4vec
);
2072 /* The force on atom ii from slab n only: */
2073 for (m
= 0; m
< DIM
; m
++)
2075 slab_force
[m
] = rotg
->k
* (-slab_sum1vec
[m
] + slab_sum2vec
[m
] - slab_sum3vec
[m
] + 0.5*slab_sum4vec
[m
]);
2078 erg
->slab_torque_v
[islab
] += torque(rotg
->vec
, slab_force
, xj
, xcn
);
2080 } /* END of loop over slabs */
2082 /* Construct the four individual parts of the vector sum: */
2083 cprod(sum1vec_part
, rotg
->vec
, sum1vec
); /* sum1vec = { } x v */
2084 cprod(sum2vec_part
, rotg
->vec
, sum2vec
); /* sum2vec = { } x v */
2085 svmul(sum3
, rotg
->vec
, sum3vec
); /* sum3vec = { } . v */
2086 svmul(sum4
, rotg
->vec
, sum4vec
); /* sum4vec = { } . v */
2088 /* Store the additional force so that it can be added to the force
2089 * array after the normal forces have been evaluated */
2090 for (m
= 0; m
< DIM
; m
++)
2092 erg
->f_rot_loc
[j
][m
] = rotg
->k
* (-sum1vec
[m
] + sum2vec
[m
] - sum3vec
[m
] + 0.5*sum4vec
[m
]);
2096 fprintf(stderr
, "sum1: %15.8f %15.8f %15.8f\n", -rotg
->k
*sum1vec
[XX
], -rotg
->k
*sum1vec
[YY
], -rotg
->k
*sum1vec
[ZZ
]);
2097 fprintf(stderr
, "sum2: %15.8f %15.8f %15.8f\n", rotg
->k
*sum2vec
[XX
], rotg
->k
*sum2vec
[YY
], rotg
->k
*sum2vec
[ZZ
]);
2098 fprintf(stderr
, "sum3: %15.8f %15.8f %15.8f\n", -rotg
->k
*sum3vec
[XX
], -rotg
->k
*sum3vec
[YY
], -rotg
->k
*sum3vec
[ZZ
]);
2099 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
]);
2104 } /* END of loop over local atoms */
2110 static real
do_flex_lowlevel(
2112 real sigma
, /* The Gaussian width sigma */
2114 gmx_bool bOutstepRot
,
2115 gmx_bool bOutstepSlab
,
2118 int count
, ic
, ifit
, ii
, j
, m
, n
, islab
, iigrp
;
2119 rvec xj
, yj0
; /* current and reference position */
2120 rvec xcn
, ycn
; /* the current and the reference slab centers */
2121 rvec yj0_ycn
; /* yj0 - ycn */
2122 rvec xj_xcn
; /* xj - xcn */
2123 rvec qjn
, fit_qjn
; /* q_i^n */
2124 rvec sum_n1
, sum_n2
; /* Two contributions to the rotation force */
2125 rvec innersumvec
; /* Inner part of sum_n2 */
2127 rvec force_n
; /* Single force from slab n on one atom */
2128 rvec force_n1
, force_n2
; /* First and second part of force_n */
2129 rvec tmpvec
, tmpvec2
, tmp_f
; /* Helper variables */
2130 real V
; /* The rotation potential energy */
2131 real OOsigma2
; /* 1/(sigma^2) */
2132 real beta
; /* beta_n(xj) */
2133 real bjn
, fit_bjn
; /* b_j^n */
2134 real gaussian_xj
; /* Gaussian weight gn(xj) */
2135 real betan_xj_sigma2
;
2136 real mj
, wj
; /* Mass-weighting of the positions */
2138 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2139 gmx_bool bCalcPotFit
;
2142 erg
= rotg
->enfrotgrp
;
2144 /* Pre-calculate the inner sums, so that we do not have to calculate
2145 * them again for every atom */
2146 flex_precalc_inner_sum(rotg
);
2148 bCalcPotFit
= (bOutstepRot
|| bOutstepSlab
) && (erotgFitPOT
== rotg
->eFittype
);
2150 /********************************************************/
2151 /* Main loop over all local atoms of the rotation group */
2152 /********************************************************/
2153 OOsigma2
= 1.0/(sigma
*sigma
);
2154 N_M
= rotg
->nat
* erg
->invmass
;
2156 for (j
= 0; j
< erg
->nat_loc
; j
++)
2158 /* Local index of a rotation group atom */
2159 ii
= erg
->ind_loc
[j
];
2160 /* Position of this atom in the collective array */
2161 iigrp
= erg
->xc_ref_ind
[j
];
2162 /* Mass-weighting */
2163 mj
= erg
->mc
[iigrp
]; /* need the unsorted mass here */
2166 /* Current position of this atom: x[ii][XX/YY/ZZ]
2167 * Note that erg->xc_center contains the center of mass in case the flex-t
2168 * potential was chosen. For the flex potential erg->xc_center must be
2170 rvec_sub(x
[ii
], erg
->xc_center
, xj
);
2172 /* Shift this atom such that it is near its reference */
2173 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
2175 /* Determine the slabs to loop over, i.e. the ones with contributions
2176 * larger than min_gaussian */
2177 count
= get_single_atom_gaussians(xj
, rotg
);
2182 /* Loop over the relevant slabs for this atom */
2183 for (ic
= 0; ic
< count
; ic
++)
2185 n
= erg
->gn_slabind
[ic
];
2187 /* Get the precomputed Gaussian for xj in slab n */
2188 gaussian_xj
= erg
->gn_atom
[ic
];
2190 islab
= n
- erg
->slab_first
; /* slab index */
2192 /* The (unrotated) reference position of this atom is saved in yj0: */
2193 copy_rvec(rotg
->x_ref
[iigrp
], yj0
);
2195 beta
= calc_beta(xj
, rotg
, n
);
2197 /* The current center of this slab is saved in xcn: */
2198 copy_rvec(erg
->slab_center
[islab
], xcn
);
2199 /* ... and the reference center in ycn: */
2200 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
2202 rvec_sub(yj0
, ycn
, yj0_ycn
); /* yj0_ycn = yj0 - ycn */
2204 /* In rare cases, when an atom position coincides with a reference slab
2205 * center (yj0_ycn == 0) we cannot compute the normal vector qjn.
2206 * However, since the atom is located directly on the pivot, this
2207 * slab's contribution to the force on that atom will be zero
2208 * anyway. Therefore, we directly move on to the next slab. */
2209 if (gmx_numzero(norm(yj0_ycn
))) /* 0 == norm(yj0 - ycn) */
2215 mvmul(erg
->rotmat
, yj0_ycn
, tmpvec2
); /* tmpvec2= Omega.(yj0-ycn) */
2217 /* Subtract the slab center from xj */
2218 rvec_sub(xj
, xcn
, xj_xcn
); /* xj_xcn = xj - xcn */
2221 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec= v x Omega.(yj0-ycn) */
2223 /* * v x Omega.(yj0-ycn) */
2224 unitv(tmpvec
, qjn
); /* qjn = --------------------- */
2225 /* |v x Omega.(yj0-ycn)| */
2227 bjn
= iprod(qjn
, xj_xcn
); /* bjn = qjn * (xj - xcn) */
2229 /*********************************/
2230 /* Add to the rotation potential */
2231 /*********************************/
2232 V
+= 0.5*rotg
->k
*wj
*gaussian_xj
*gmx::square(bjn
);
2234 /* If requested, also calculate the potential for a set of angles
2235 * near the current reference angle */
2238 for (ifit
= 0; ifit
< rotg
->PotAngle_nstep
; ifit
++)
2240 /* As above calculate Omega.(yj0-ycn), now for the other angles */
2241 mvmul(erg
->PotAngleFit
->rotmat
[ifit
], yj0_ycn
, tmpvec2
); /* tmpvec2= Omega.(yj0-ycn) */
2242 /* As above calculate qjn */
2243 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec= v x Omega.(yj0-ycn) */
2244 /* * v x Omega.(yj0-ycn) */
2245 unitv(tmpvec
, fit_qjn
); /* fit_qjn = --------------------- */
2246 /* |v x Omega.(yj0-ycn)| */
2247 fit_bjn
= iprod(fit_qjn
, xj_xcn
); /* fit_bjn = fit_qjn * (xj - xcn) */
2248 /* Add to the rotation potential for this angle */
2249 erg
->PotAngleFit
->V
[ifit
] += 0.5*rotg
->k
*wj
*gaussian_xj
*gmx::square(fit_bjn
);
2253 /****************************************************************/
2254 /* sum_n1 will typically be the main contribution to the force: */
2255 /****************************************************************/
2256 betan_xj_sigma2
= beta
*OOsigma2
; /* beta_n(xj)/sigma^2 */
2258 /* The next lines calculate
2259 * qjn - (bjn*beta(xj)/(2sigma^2))v */
2260 svmul(bjn
*0.5*betan_xj_sigma2
, rotg
->vec
, tmpvec2
);
2261 rvec_sub(qjn
, tmpvec2
, tmpvec
);
2263 /* Multiply with gn(xj)*bjn: */
2264 svmul(gaussian_xj
*bjn
, tmpvec
, tmpvec2
);
2267 rvec_inc(sum_n1
, tmpvec2
);
2269 /* We already have precalculated the Sn term for slab n */
2270 copy_rvec(erg
->slab_innersumvec
[islab
], s_n
);
2272 svmul(betan_xj_sigma2
*iprod(s_n
, xj_xcn
), rotg
->vec
, tmpvec
); /* tmpvec = ---------- s_n (xj-xcn) */
2275 rvec_sub(s_n
, tmpvec
, innersumvec
);
2277 /* We can safely divide by slab_weights since we check in get_slab_centers
2278 * that it is non-zero. */
2279 svmul(gaussian_xj
/erg
->slab_weights
[islab
], innersumvec
, innersumvec
);
2281 rvec_add(sum_n2
, innersumvec
, sum_n2
);
2283 /* Calculate the torque: */
2286 /* The force on atom ii from slab n only: */
2287 svmul(-rotg
->k
*wj
, tmpvec2
, force_n1
); /* part 1 */
2288 svmul( rotg
->k
*mj
, innersumvec
, force_n2
); /* part 2 */
2289 rvec_add(force_n1
, force_n2
, force_n
);
2290 erg
->slab_torque_v
[islab
] += torque(rotg
->vec
, force_n
, xj
, xcn
);
2292 } /* END of loop over slabs */
2294 /* Put both contributions together: */
2295 svmul(wj
, sum_n1
, sum_n1
);
2296 svmul(mj
, sum_n2
, sum_n2
);
2297 rvec_sub(sum_n2
, sum_n1
, tmp_f
); /* F = -grad V */
2299 /* Store the additional force so that it can be added to the force
2300 * array after the normal forces have been evaluated */
2301 for (m
= 0; m
< DIM
; m
++)
2303 erg
->f_rot_loc
[j
][m
] = rotg
->k
*tmp_f
[m
];
2308 } /* END of loop over local atoms */
2314 static void print_coordinates(t_rotgrp
*rotg
, rvec x
[], matrix box
, int step
)
2318 static char buf
[STRLEN
];
2319 static gmx_bool bFirst
= 1;
2324 sprintf(buf
, "coords%d.txt", cr
->nodeid
);
2325 fp
= fopen(buf
, "w");
2329 fprintf(fp
, "\nStep %d\n", step
);
2330 fprintf(fp
, "box: %f %f %f %f %f %f %f %f %f\n",
2331 box
[XX
][XX
], box
[XX
][YY
], box
[XX
][ZZ
],
2332 box
[YY
][XX
], box
[YY
][YY
], box
[YY
][ZZ
],
2333 box
[ZZ
][XX
], box
[ZZ
][ZZ
], box
[ZZ
][ZZ
]);
2334 for (i
= 0; i
< rotg
->nat
; i
++)
2336 fprintf(fp
, "%4d %f %f %f\n", i
,
2337 erg
->xc
[i
][XX
], erg
->xc
[i
][YY
], erg
->xc
[i
][ZZ
]);
2345 static int projection_compare(const void *a
, const void *b
)
2347 sort_along_vec_t
*xca
, *xcb
;
2350 xca
= (sort_along_vec_t
*)a
;
2351 xcb
= (sort_along_vec_t
*)b
;
2353 if (xca
->xcproj
< xcb
->xcproj
)
2357 else if (xca
->xcproj
> xcb
->xcproj
)
2368 static void sort_collective_coordinates(
2369 t_rotgrp
*rotg
, /* Rotation group */
2370 sort_along_vec_t
*data
) /* Buffer for sorting the positions */
2373 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2376 erg
= rotg
->enfrotgrp
;
2378 /* The projection of the position vector on the rotation vector is
2379 * the relevant value for sorting. Fill the 'data' structure */
2380 for (i
= 0; i
< rotg
->nat
; i
++)
2382 data
[i
].xcproj
= iprod(erg
->xc
[i
], rotg
->vec
); /* sort criterium */
2383 data
[i
].m
= erg
->mc
[i
];
2385 copy_rvec(erg
->xc
[i
], data
[i
].x
);
2386 copy_rvec(rotg
->x_ref
[i
], data
[i
].x_ref
);
2388 /* Sort the 'data' structure */
2389 gmx_qsort(data
, rotg
->nat
, sizeof(sort_along_vec_t
), projection_compare
);
2391 /* Copy back the sorted values */
2392 for (i
= 0; i
< rotg
->nat
; i
++)
2394 copy_rvec(data
[i
].x
, erg
->xc
[i
] );
2395 copy_rvec(data
[i
].x_ref
, erg
->xc_ref_sorted
[i
]);
2396 erg
->mc_sorted
[i
] = data
[i
].m
;
2397 erg
->xc_sortind
[i
] = data
[i
].ind
;
2402 /* For each slab, get the first and the last index of the sorted atom
2404 static void get_firstlast_atom_per_slab(t_rotgrp
*rotg
)
2408 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2411 erg
= rotg
->enfrotgrp
;
2413 /* Find the first atom that needs to enter the calculation for each slab */
2414 n
= erg
->slab_first
; /* slab */
2415 i
= 0; /* start with the first atom */
2418 /* Find the first atom that significantly contributes to this slab */
2419 do /* move forward in position until a large enough beta is found */
2421 beta
= calc_beta(erg
->xc
[i
], rotg
, n
);
2424 while ((beta
< -erg
->max_beta
) && (i
< rotg
->nat
));
2426 islab
= n
- erg
->slab_first
; /* slab index */
2427 erg
->firstatom
[islab
] = i
;
2428 /* Proceed to the next slab */
2431 while (n
<= erg
->slab_last
);
2433 /* Find the last atom for each slab */
2434 n
= erg
->slab_last
; /* start with last slab */
2435 i
= rotg
->nat
-1; /* start with the last atom */
2438 do /* move backward in position until a large enough beta is found */
2440 beta
= calc_beta(erg
->xc
[i
], rotg
, n
);
2443 while ((beta
> erg
->max_beta
) && (i
> -1));
2445 islab
= n
- erg
->slab_first
; /* slab index */
2446 erg
->lastatom
[islab
] = i
;
2447 /* Proceed to the next slab */
2450 while (n
>= erg
->slab_first
);
2454 /* Determine the very first and very last slab that needs to be considered
2455 * For the first slab that needs to be considered, we have to find the smallest
2458 * x_first * v - n*Delta_x <= beta_max
2460 * slab index n, slab distance Delta_x, rotation vector v. For the last slab we
2461 * have to find the largest n that obeys
2463 * x_last * v - n*Delta_x >= -beta_max
2466 static gmx_inline
int get_first_slab(
2467 t_rotgrp
*rotg
, /* The rotation group (inputrec data) */
2468 real max_beta
, /* The max_beta value, instead of min_gaussian */
2469 rvec firstatom
) /* First atom after sorting along the rotation vector v */
2471 /* Find the first slab for the first atom */
2472 return static_cast<int>(ceil(static_cast<double>((iprod(firstatom
, rotg
->vec
) - max_beta
)/rotg
->slab_dist
)));
2476 static gmx_inline
int get_last_slab(
2477 t_rotgrp
*rotg
, /* The rotation group (inputrec data) */
2478 real max_beta
, /* The max_beta value, instead of min_gaussian */
2479 rvec lastatom
) /* Last atom along v */
2481 /* Find the last slab for the last atom */
2482 return static_cast<int>(floor(static_cast<double>((iprod(lastatom
, rotg
->vec
) + max_beta
)/rotg
->slab_dist
)));
2486 static void get_firstlast_slab_check(
2487 t_rotgrp
*rotg
, /* The rotation group (inputrec data) */
2488 t_gmx_enfrotgrp
*erg
, /* The rotation group (data only accessible in this file) */
2489 rvec firstatom
, /* First atom after sorting along the rotation vector v */
2490 rvec lastatom
) /* Last atom along v */
2492 erg
->slab_first
= get_first_slab(rotg
, erg
->max_beta
, firstatom
);
2493 erg
->slab_last
= get_last_slab(rotg
, erg
->max_beta
, lastatom
);
2495 /* Calculate the slab buffer size, which changes when slab_first changes */
2496 erg
->slab_buffer
= erg
->slab_first
- erg
->slab_first_ref
;
2498 /* Check whether we have reference data to compare against */
2499 if (erg
->slab_first
< erg
->slab_first_ref
)
2501 gmx_fatal(FARGS
, "%s No reference data for first slab (n=%d), unable to proceed.",
2502 RotStr
, erg
->slab_first
);
2505 /* Check whether we have reference data to compare against */
2506 if (erg
->slab_last
> erg
->slab_last_ref
)
2508 gmx_fatal(FARGS
, "%s No reference data for last slab (n=%d), unable to proceed.",
2509 RotStr
, erg
->slab_last
);
2514 /* Enforced rotation with a flexible axis */
2515 static void do_flexible(
2517 gmx_enfrot_t enfrot
, /* Other rotation data */
2518 t_rotgrp
*rotg
, /* The rotation group */
2519 int g
, /* Group number */
2520 rvec x
[], /* The local positions */
2522 double t
, /* Time in picoseconds */
2523 gmx_bool bOutstepRot
, /* Output to main rotation output file */
2524 gmx_bool bOutstepSlab
) /* Output per-slab data */
2527 real sigma
; /* The Gaussian width sigma */
2528 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2531 erg
= rotg
->enfrotgrp
;
2533 /* Define the sigma value */
2534 sigma
= 0.7*rotg
->slab_dist
;
2536 /* Sort the collective coordinates erg->xc along the rotation vector. This is
2537 * an optimization for the inner loop. */
2538 sort_collective_coordinates(rotg
, enfrot
->data
);
2540 /* Determine the first relevant slab for the first atom and the last
2541 * relevant slab for the last atom */
2542 get_firstlast_slab_check(rotg
, erg
, erg
->xc
[0], erg
->xc
[rotg
->nat
-1]);
2544 /* Determine for each slab depending on the min_gaussian cutoff criterium,
2545 * a first and a last atom index inbetween stuff needs to be calculated */
2546 get_firstlast_atom_per_slab(rotg
);
2548 /* Determine the gaussian-weighted center of positions for all slabs */
2549 get_slab_centers(rotg
, erg
->xc
, erg
->mc_sorted
, g
, t
, enfrot
->out_slabs
, bOutstepSlab
, FALSE
);
2551 /* Clear the torque per slab from last time step: */
2552 nslabs
= erg
->slab_last
- erg
->slab_first
+ 1;
2553 for (l
= 0; l
< nslabs
; l
++)
2555 erg
->slab_torque_v
[l
] = 0.0;
2558 /* Call the rotational forces kernel */
2559 if (rotg
->eType
== erotgFLEX
|| rotg
->eType
== erotgFLEXT
)
2561 erg
->V
= do_flex_lowlevel(rotg
, sigma
, x
, bOutstepRot
, bOutstepSlab
, box
);
2563 else if (rotg
->eType
== erotgFLEX2
|| rotg
->eType
== erotgFLEX2T
)
2565 erg
->V
= do_flex2_lowlevel(rotg
, sigma
, x
, bOutstepRot
, bOutstepSlab
, box
);
2569 gmx_fatal(FARGS
, "Unknown flexible rotation type");
2572 /* Determine angle by RMSD fit to the reference - Let's hope this */
2573 /* only happens once in a while, since this is not parallelized! */
2574 if (bMaster
&& (erotgFitPOT
!= rotg
->eFittype
) )
2578 /* Fit angle of the whole rotation group */
2579 erg
->angle_v
= flex_fit_angle(rotg
);
2583 /* Fit angle of each slab */
2584 flex_fit_angle_perslab(g
, rotg
, t
, erg
->degangle
, enfrot
->out_angles
);
2588 /* Lump together the torques from all slabs: */
2589 erg
->torque_v
= 0.0;
2590 for (l
= 0; l
< nslabs
; l
++)
2592 erg
->torque_v
+= erg
->slab_torque_v
[l
];
2597 /* Calculate the angle between reference and actual rotation group atom,
2598 * both projected into a plane perpendicular to the rotation vector: */
2599 static void angle(t_rotgrp
*rotg
,
2603 real
*weight
) /* atoms near the rotation axis should count less than atoms far away */
2605 rvec xp
, xrp
; /* current and reference positions projected on a plane perpendicular to pg->vec */
2609 /* Project x_ref and x into a plane through the origin perpendicular to rot_vec: */
2610 /* Project x_ref: xrp = x_ref - (vec * x_ref) * vec */
2611 svmul(iprod(rotg
->vec
, x_ref
), rotg
->vec
, dum
);
2612 rvec_sub(x_ref
, dum
, xrp
);
2613 /* Project x_act: */
2614 svmul(iprod(rotg
->vec
, x_act
), rotg
->vec
, dum
);
2615 rvec_sub(x_act
, dum
, xp
);
2617 /* Retrieve information about which vector precedes. gmx_angle always
2618 * returns a positive angle. */
2619 cprod(xp
, xrp
, dum
); /* if reference precedes, this is pointing into the same direction as vec */
2621 if (iprod(rotg
->vec
, dum
) >= 0)
2623 *alpha
= -gmx_angle(xrp
, xp
);
2627 *alpha
= +gmx_angle(xrp
, xp
);
2630 /* Also return the weight */
2635 /* Project first vector onto a plane perpendicular to the second vector
2637 * Note that v must be of unit length.
2639 static gmx_inline
void project_onto_plane(rvec dr
, const rvec v
)
2644 svmul(iprod(dr
, v
), v
, tmp
); /* tmp = (dr.v)v */
2645 rvec_dec(dr
, tmp
); /* dr = dr - (dr.v)v */
2649 /* Fixed rotation: The rotation reference group rotates around the v axis. */
2650 /* The atoms of the actual rotation group are attached with imaginary */
2651 /* springs to the reference atoms. */
2652 static void do_fixed(
2653 t_rotgrp
*rotg
, /* The rotation group */
2654 gmx_bool bOutstepRot
, /* Output to main rotation output file */
2655 gmx_bool bOutstepSlab
) /* Output per-slab data */
2659 rvec tmp_f
; /* Force */
2660 real alpha
; /* a single angle between an actual and a reference position */
2661 real weight
; /* single weight for a single angle */
2662 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2663 rvec xi_xc
; /* xi - xc */
2664 gmx_bool bCalcPotFit
;
2667 /* for mass weighting: */
2668 real wi
; /* Mass-weighting of the positions */
2670 real k_wi
; /* k times wi */
2675 erg
= rotg
->enfrotgrp
;
2676 bProject
= (rotg
->eType
== erotgPM
) || (rotg
->eType
== erotgPMPF
);
2677 bCalcPotFit
= (bOutstepRot
|| bOutstepSlab
) && (erotgFitPOT
== rotg
->eFittype
);
2679 N_M
= rotg
->nat
* erg
->invmass
;
2681 /* Each process calculates the forces on its local atoms */
2682 for (j
= 0; j
< erg
->nat_loc
; j
++)
2684 /* Calculate (x_i-x_c) resp. (x_i-u) */
2685 rvec_sub(erg
->x_loc_pbc
[j
], erg
->xc_center
, xi_xc
);
2687 /* Calculate Omega*(y_i-y_c)-(x_i-x_c) */
2688 rvec_sub(erg
->xr_loc
[j
], xi_xc
, dr
);
2692 project_onto_plane(dr
, rotg
->vec
);
2695 /* Mass-weighting */
2696 wi
= N_M
*erg
->m_loc
[j
];
2698 /* Store the additional force so that it can be added to the force
2699 * array after the normal forces have been evaluated */
2701 for (m
= 0; m
< DIM
; m
++)
2703 tmp_f
[m
] = k_wi
*dr
[m
];
2704 erg
->f_rot_loc
[j
][m
] = tmp_f
[m
];
2705 erg
->V
+= 0.5*k_wi
*gmx::square(dr
[m
]);
2708 /* If requested, also calculate the potential for a set of angles
2709 * near the current reference angle */
2712 for (ifit
= 0; ifit
< rotg
->PotAngle_nstep
; ifit
++)
2714 /* Index of this rotation group atom with respect to the whole rotation group */
2715 jj
= erg
->xc_ref_ind
[j
];
2717 /* Rotate with the alternative angle. Like rotate_local_reference(),
2718 * just for a single local atom */
2719 mvmul(erg
->PotAngleFit
->rotmat
[ifit
], rotg
->x_ref
[jj
], fit_xr_loc
); /* fit_xr_loc = Omega*(y_i-y_c) */
2721 /* Calculate Omega*(y_i-y_c)-(x_i-x_c) */
2722 rvec_sub(fit_xr_loc
, xi_xc
, dr
);
2726 project_onto_plane(dr
, rotg
->vec
);
2729 /* Add to the rotation potential for this angle: */
2730 erg
->PotAngleFit
->V
[ifit
] += 0.5*k_wi
*norm2(dr
);
2736 /* Add to the torque of this rotation group */
2737 erg
->torque_v
+= torque(rotg
->vec
, tmp_f
, erg
->x_loc_pbc
[j
], erg
->xc_center
);
2739 /* Calculate the angle between reference and actual rotation group atom. */
2740 angle(rotg
, xi_xc
, erg
->xr_loc
[j
], &alpha
, &weight
); /* angle in rad, weighted */
2741 erg
->angle_v
+= alpha
* weight
;
2742 erg
->weight_v
+= weight
;
2744 /* If you want enforced rotation to contribute to the virial,
2745 * activate the following lines:
2748 Add the rotation contribution to the virial
2749 for(j=0; j<DIM; j++)
2751 vir[j][m] += 0.5*f[ii][j]*dr[m];
2757 } /* end of loop over local rotation group atoms */
2761 /* Calculate the radial motion potential and forces */
2762 static void do_radial_motion(
2763 t_rotgrp
*rotg
, /* The rotation group */
2764 gmx_bool bOutstepRot
, /* Output to main rotation output file */
2765 gmx_bool bOutstepSlab
) /* Output per-slab data */
2768 rvec tmp_f
; /* Force */
2769 real alpha
; /* a single angle between an actual and a reference position */
2770 real weight
; /* single weight for a single angle */
2771 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2772 rvec xj_u
; /* xj - u */
2773 rvec tmpvec
, fit_tmpvec
;
2774 real fac
, fac2
, sum
= 0.0;
2776 gmx_bool bCalcPotFit
;
2778 /* For mass weighting: */
2779 real wj
; /* Mass-weighting of the positions */
2783 erg
= rotg
->enfrotgrp
;
2784 bCalcPotFit
= (bOutstepRot
|| bOutstepSlab
) && (erotgFitPOT
== rotg
->eFittype
);
2786 N_M
= rotg
->nat
* erg
->invmass
;
2788 /* Each process calculates the forces on its local atoms */
2789 for (j
= 0; j
< erg
->nat_loc
; j
++)
2791 /* Calculate (xj-u) */
2792 rvec_sub(erg
->x_loc_pbc
[j
], erg
->xc_center
, xj_u
); /* xj_u = xj-u */
2794 /* Calculate Omega.(yj0-u) */
2795 cprod(rotg
->vec
, erg
->xr_loc
[j
], tmpvec
); /* tmpvec = v x Omega.(yj0-u) */
2797 /* * v x Omega.(yj0-u) */
2798 unitv(tmpvec
, pj
); /* pj = --------------------- */
2799 /* | v x Omega.(yj0-u) | */
2801 fac
= iprod(pj
, xj_u
); /* fac = pj.(xj-u) */
2804 /* Mass-weighting */
2805 wj
= N_M
*erg
->m_loc
[j
];
2807 /* Store the additional force so that it can be added to the force
2808 * array after the normal forces have been evaluated */
2809 svmul(-rotg
->k
*wj
*fac
, pj
, tmp_f
);
2810 copy_rvec(tmp_f
, erg
->f_rot_loc
[j
]);
2813 /* If requested, also calculate the potential for a set of angles
2814 * near the current reference angle */
2817 for (ifit
= 0; ifit
< rotg
->PotAngle_nstep
; ifit
++)
2819 /* Index of this rotation group atom with respect to the whole rotation group */
2820 jj
= erg
->xc_ref_ind
[j
];
2822 /* Rotate with the alternative angle. Like rotate_local_reference(),
2823 * just for a single local atom */
2824 mvmul(erg
->PotAngleFit
->rotmat
[ifit
], rotg
->x_ref
[jj
], fit_tmpvec
); /* fit_tmpvec = Omega*(yj0-u) */
2826 /* Calculate Omega.(yj0-u) */
2827 cprod(rotg
->vec
, fit_tmpvec
, tmpvec
); /* tmpvec = v x Omega.(yj0-u) */
2828 /* * v x Omega.(yj0-u) */
2829 unitv(tmpvec
, pj
); /* pj = --------------------- */
2830 /* | v x Omega.(yj0-u) | */
2832 fac
= iprod(pj
, xj_u
); /* fac = pj.(xj-u) */
2835 /* Add to the rotation potential for this angle: */
2836 erg
->PotAngleFit
->V
[ifit
] += 0.5*rotg
->k
*wj
*fac2
;
2842 /* Add to the torque of this rotation group */
2843 erg
->torque_v
+= torque(rotg
->vec
, tmp_f
, erg
->x_loc_pbc
[j
], erg
->xc_center
);
2845 /* Calculate the angle between reference and actual rotation group atom. */
2846 angle(rotg
, xj_u
, erg
->xr_loc
[j
], &alpha
, &weight
); /* angle in rad, weighted */
2847 erg
->angle_v
+= alpha
* weight
;
2848 erg
->weight_v
+= weight
;
2853 } /* end of loop over local rotation group atoms */
2854 erg
->V
= 0.5*rotg
->k
*sum
;
2858 /* Calculate the radial motion pivot-free potential and forces */
2859 static void do_radial_motion_pf(
2860 t_rotgrp
*rotg
, /* The rotation group */
2861 rvec x
[], /* The positions */
2862 matrix box
, /* The simulation box */
2863 gmx_bool bOutstepRot
, /* Output to main rotation output file */
2864 gmx_bool bOutstepSlab
) /* Output per-slab data */
2866 int i
, ii
, iigrp
, ifit
, j
;
2867 rvec xj
; /* Current position */
2868 rvec xj_xc
; /* xj - xc */
2869 rvec yj0_yc0
; /* yj0 - yc0 */
2870 rvec tmp_f
; /* Force */
2871 real alpha
; /* a single angle between an actual and a reference position */
2872 real weight
; /* single weight for a single angle */
2873 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2874 rvec tmpvec
, tmpvec2
;
2875 rvec innersumvec
; /* Precalculation of the inner sum */
2877 real fac
, fac2
, V
= 0.0;
2879 gmx_bool bCalcPotFit
;
2881 /* For mass weighting: */
2882 real mj
, wi
, wj
; /* Mass-weighting of the positions */
2886 erg
= rotg
->enfrotgrp
;
2887 bCalcPotFit
= (bOutstepRot
|| bOutstepSlab
) && (erotgFitPOT
== rotg
->eFittype
);
2889 N_M
= rotg
->nat
* erg
->invmass
;
2891 /* Get the current center of the rotation group: */
2892 get_center(erg
->xc
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
2894 /* Precalculate Sum_i [ wi qi.(xi-xc) qi ] which is needed for every single j */
2895 clear_rvec(innersumvec
);
2896 for (i
= 0; i
< rotg
->nat
; i
++)
2898 /* Mass-weighting */
2899 wi
= N_M
*erg
->mc
[i
];
2901 /* Calculate qi. Note that xc_ref_center has already been subtracted from
2902 * x_ref in init_rot_group.*/
2903 mvmul(erg
->rotmat
, rotg
->x_ref
[i
], tmpvec
); /* tmpvec = Omega.(yi0-yc0) */
2905 cprod(rotg
->vec
, tmpvec
, tmpvec2
); /* tmpvec2 = v x Omega.(yi0-yc0) */
2907 /* * v x Omega.(yi0-yc0) */
2908 unitv(tmpvec2
, qi
); /* qi = ----------------------- */
2909 /* | v x Omega.(yi0-yc0) | */
2911 rvec_sub(erg
->xc
[i
], erg
->xc_center
, tmpvec
); /* tmpvec = xi-xc */
2913 svmul(wi
*iprod(qi
, tmpvec
), qi
, tmpvec2
);
2915 rvec_inc(innersumvec
, tmpvec2
);
2917 svmul(rotg
->k
*erg
->invmass
, innersumvec
, innersumveckM
);
2919 /* Each process calculates the forces on its local atoms */
2920 for (j
= 0; j
< erg
->nat_loc
; j
++)
2922 /* Local index of a rotation group atom */
2923 ii
= erg
->ind_loc
[j
];
2924 /* Position of this atom in the collective array */
2925 iigrp
= erg
->xc_ref_ind
[j
];
2926 /* Mass-weighting */
2927 mj
= erg
->mc
[iigrp
]; /* need the unsorted mass here */
2930 /* Current position of this atom: x[ii][XX/YY/ZZ] */
2931 copy_rvec(x
[ii
], xj
);
2933 /* Shift this atom such that it is near its reference */
2934 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
2936 /* The (unrotated) reference position is yj0. yc0 has already
2937 * been subtracted in init_rot_group */
2938 copy_rvec(rotg
->x_ref
[iigrp
], yj0_yc0
); /* yj0_yc0 = yj0 - yc0 */
2940 /* Calculate Omega.(yj0-yc0) */
2941 mvmul(erg
->rotmat
, yj0_yc0
, tmpvec2
); /* tmpvec2 = Omega.(yj0 - yc0) */
2943 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x Omega.(yj0-yc0) */
2945 /* * v x Omega.(yj0-yc0) */
2946 unitv(tmpvec
, qj
); /* qj = ----------------------- */
2947 /* | v x Omega.(yj0-yc0) | */
2949 /* Calculate (xj-xc) */
2950 rvec_sub(xj
, erg
->xc_center
, xj_xc
); /* xj_xc = xj-xc */
2952 fac
= iprod(qj
, xj_xc
); /* fac = qj.(xj-xc) */
2955 /* Store the additional force so that it can be added to the force
2956 * array after the normal forces have been evaluated */
2957 svmul(-rotg
->k
*wj
*fac
, qj
, tmp_f
); /* part 1 of force */
2958 svmul(mj
, innersumveckM
, tmpvec
); /* part 2 of force */
2959 rvec_inc(tmp_f
, tmpvec
);
2960 copy_rvec(tmp_f
, erg
->f_rot_loc
[j
]);
2963 /* If requested, also calculate the potential for a set of angles
2964 * near the current reference angle */
2967 for (ifit
= 0; ifit
< rotg
->PotAngle_nstep
; ifit
++)
2969 /* Rotate with the alternative angle. Like rotate_local_reference(),
2970 * just for a single local atom */
2971 mvmul(erg
->PotAngleFit
->rotmat
[ifit
], yj0_yc0
, tmpvec2
); /* tmpvec2 = Omega*(yj0-yc0) */
2973 /* Calculate Omega.(yj0-u) */
2974 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x Omega.(yj0-yc0) */
2975 /* * v x Omega.(yj0-yc0) */
2976 unitv(tmpvec
, qj
); /* qj = ----------------------- */
2977 /* | v x Omega.(yj0-yc0) | */
2979 fac
= iprod(qj
, xj_xc
); /* fac = qj.(xj-xc) */
2982 /* Add to the rotation potential for this angle: */
2983 erg
->PotAngleFit
->V
[ifit
] += 0.5*rotg
->k
*wj
*fac2
;
2989 /* Add to the torque of this rotation group */
2990 erg
->torque_v
+= torque(rotg
->vec
, tmp_f
, xj
, erg
->xc_center
);
2992 /* Calculate the angle between reference and actual rotation group atom. */
2993 angle(rotg
, xj_xc
, yj0_yc0
, &alpha
, &weight
); /* angle in rad, weighted */
2994 erg
->angle_v
+= alpha
* weight
;
2995 erg
->weight_v
+= weight
;
3000 } /* end of loop over local rotation group atoms */
3001 erg
->V
= 0.5*rotg
->k
*V
;
3005 /* Precalculate the inner sum for the radial motion 2 forces */
3006 static void radial_motion2_precalc_inner_sum(t_rotgrp
*rotg
, rvec innersumvec
)
3009 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3010 rvec xi_xc
; /* xj - xc */
3011 rvec tmpvec
, tmpvec2
;
3015 rvec v_xi_xc
; /* v x (xj - u) */
3016 real psii
, psiistar
;
3017 real wi
; /* Mass-weighting of the positions */
3021 erg
= rotg
->enfrotgrp
;
3022 N_M
= rotg
->nat
* erg
->invmass
;
3024 /* Loop over the collective set of positions */
3026 for (i
= 0; i
< rotg
->nat
; i
++)
3028 /* Mass-weighting */
3029 wi
= N_M
*erg
->mc
[i
];
3031 rvec_sub(erg
->xc
[i
], erg
->xc_center
, xi_xc
); /* xi_xc = xi-xc */
3033 /* Calculate ri. Note that xc_ref_center has already been subtracted from
3034 * x_ref in init_rot_group.*/
3035 mvmul(erg
->rotmat
, rotg
->x_ref
[i
], ri
); /* ri = Omega.(yi0-yc0) */
3037 cprod(rotg
->vec
, xi_xc
, v_xi_xc
); /* v_xi_xc = v x (xi-u) */
3039 fac
= norm2(v_xi_xc
);
3041 psiistar
= 1.0/(fac
+ rotg
->eps
); /* psiistar = --------------------- */
3042 /* |v x (xi-xc)|^2 + eps */
3044 psii
= gmx::invsqrt(fac
); /* 1 */
3045 /* psii = ------------- */
3048 svmul(psii
, v_xi_xc
, si
); /* si = psii * (v x (xi-xc) ) */
3050 siri
= iprod(si
, ri
); /* siri = si.ri */
3052 svmul(psiistar
/psii
, ri
, tmpvec
);
3053 svmul(psiistar
*psiistar
/(psii
*psii
*psii
) * siri
, si
, tmpvec2
);
3054 rvec_dec(tmpvec
, tmpvec2
);
3055 cprod(tmpvec
, rotg
->vec
, tmpvec2
);
3057 svmul(wi
*siri
, tmpvec2
, tmpvec
);
3059 rvec_inc(sumvec
, tmpvec
);
3061 svmul(rotg
->k
*erg
->invmass
, sumvec
, innersumvec
);
3065 /* Calculate the radial motion 2 potential and forces */
3066 static void do_radial_motion2(
3067 t_rotgrp
*rotg
, /* The rotation group */
3068 rvec x
[], /* The positions */
3069 matrix box
, /* The simulation box */
3070 gmx_bool bOutstepRot
, /* Output to main rotation output file */
3071 gmx_bool bOutstepSlab
) /* Output per-slab data */
3073 int ii
, iigrp
, ifit
, j
;
3074 rvec xj
; /* Position */
3075 real alpha
; /* a single angle between an actual and a reference position */
3076 real weight
; /* single weight for a single angle */
3077 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3078 rvec xj_u
; /* xj - u */
3079 rvec yj0_yc0
; /* yj0 -yc0 */
3080 rvec tmpvec
, tmpvec2
;
3081 real fac
, fit_fac
, fac2
, Vpart
= 0.0;
3082 rvec rj
, fit_rj
, sj
;
3084 rvec v_xj_u
; /* v x (xj - u) */
3085 real psij
, psijstar
;
3086 real mj
, wj
; /* For mass-weighting of the positions */
3090 gmx_bool bCalcPotFit
;
3093 erg
= rotg
->enfrotgrp
;
3095 bPF
= rotg
->eType
== erotgRM2PF
;
3096 bCalcPotFit
= (bOutstepRot
|| bOutstepSlab
) && (erotgFitPOT
== rotg
->eFittype
);
3099 clear_rvec(yj0_yc0
); /* Make the compiler happy */
3101 clear_rvec(innersumvec
);
3104 /* For the pivot-free variant we have to use the current center of
3105 * mass of the rotation group instead of the pivot u */
3106 get_center(erg
->xc
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
3108 /* Also, we precalculate the second term of the forces that is identical
3109 * (up to the weight factor mj) for all forces */
3110 radial_motion2_precalc_inner_sum(rotg
, innersumvec
);
3113 N_M
= rotg
->nat
* erg
->invmass
;
3115 /* Each process calculates the forces on its local atoms */
3116 for (j
= 0; j
< erg
->nat_loc
; j
++)
3120 /* Local index of a rotation group atom */
3121 ii
= erg
->ind_loc
[j
];
3122 /* Position of this atom in the collective array */
3123 iigrp
= erg
->xc_ref_ind
[j
];
3124 /* Mass-weighting */
3125 mj
= erg
->mc
[iigrp
];
3127 /* Current position of this atom: x[ii] */
3128 copy_rvec(x
[ii
], xj
);
3130 /* Shift this atom such that it is near its reference */
3131 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
3133 /* The (unrotated) reference position is yj0. yc0 has already
3134 * been subtracted in init_rot_group */
3135 copy_rvec(rotg
->x_ref
[iigrp
], yj0_yc0
); /* yj0_yc0 = yj0 - yc0 */
3137 /* Calculate Omega.(yj0-yc0) */
3138 mvmul(erg
->rotmat
, yj0_yc0
, rj
); /* rj = Omega.(yj0-yc0) */
3143 copy_rvec(erg
->x_loc_pbc
[j
], xj
);
3144 copy_rvec(erg
->xr_loc
[j
], rj
); /* rj = Omega.(yj0-u) */
3146 /* Mass-weighting */
3149 /* Calculate (xj-u) resp. (xj-xc) */
3150 rvec_sub(xj
, erg
->xc_center
, xj_u
); /* xj_u = xj-u */
3152 cprod(rotg
->vec
, xj_u
, v_xj_u
); /* v_xj_u = v x (xj-u) */
3154 fac
= norm2(v_xj_u
);
3156 psijstar
= 1.0/(fac
+ rotg
->eps
); /* psistar = -------------------- */
3157 /* |v x (xj-u)|^2 + eps */
3159 psij
= gmx::invsqrt(fac
); /* 1 */
3160 /* psij = ------------ */
3163 svmul(psij
, v_xj_u
, sj
); /* sj = psij * (v x (xj-u) ) */
3165 fac
= iprod(v_xj_u
, rj
); /* fac = (v x (xj-u)).rj */
3168 sjrj
= iprod(sj
, rj
); /* sjrj = sj.rj */
3170 svmul(psijstar
/psij
, rj
, tmpvec
);
3171 svmul(psijstar
*psijstar
/(psij
*psij
*psij
) * sjrj
, sj
, tmpvec2
);
3172 rvec_dec(tmpvec
, tmpvec2
);
3173 cprod(tmpvec
, rotg
->vec
, tmpvec2
);
3175 /* Store the additional force so that it can be added to the force
3176 * array after the normal forces have been evaluated */
3177 svmul(-rotg
->k
*wj
*sjrj
, tmpvec2
, tmpvec
);
3178 svmul(mj
, innersumvec
, tmpvec2
); /* This is != 0 only for the pivot-free variant */
3180 rvec_add(tmpvec2
, tmpvec
, erg
->f_rot_loc
[j
]);
3181 Vpart
+= wj
*psijstar
*fac2
;
3183 /* If requested, also calculate the potential for a set of angles
3184 * near the current reference angle */
3187 for (ifit
= 0; ifit
< rotg
->PotAngle_nstep
; ifit
++)
3191 mvmul(erg
->PotAngleFit
->rotmat
[ifit
], yj0_yc0
, fit_rj
); /* fit_rj = Omega.(yj0-yc0) */
3195 /* Position of this atom in the collective array */
3196 iigrp
= erg
->xc_ref_ind
[j
];
3197 /* Rotate with the alternative angle. Like rotate_local_reference(),
3198 * just for a single local atom */
3199 mvmul(erg
->PotAngleFit
->rotmat
[ifit
], rotg
->x_ref
[iigrp
], fit_rj
); /* fit_rj = Omega*(yj0-u) */
3201 fit_fac
= iprod(v_xj_u
, fit_rj
); /* fac = (v x (xj-u)).fit_rj */
3202 /* Add to the rotation potential for this angle: */
3203 erg
->PotAngleFit
->V
[ifit
] += 0.5*rotg
->k
*wj
*psijstar
*fit_fac
*fit_fac
;
3209 /* Add to the torque of this rotation group */
3210 erg
->torque_v
+= torque(rotg
->vec
, erg
->f_rot_loc
[j
], xj
, erg
->xc_center
);
3212 /* Calculate the angle between reference and actual rotation group atom. */
3213 angle(rotg
, xj_u
, rj
, &alpha
, &weight
); /* angle in rad, weighted */
3214 erg
->angle_v
+= alpha
* weight
;
3215 erg
->weight_v
+= weight
;
3220 } /* end of loop over local rotation group atoms */
3221 erg
->V
= 0.5*rotg
->k
*Vpart
;
3225 /* Determine the smallest and largest position vector (with respect to the
3226 * rotation vector) for the reference group */
3227 static void get_firstlast_atom_ref(
3233 real xcproj
; /* The projection of a reference position on the
3235 real minproj
, maxproj
; /* Smallest and largest projection on v */
3239 /* Start with some value */
3240 minproj
= iprod(rotg
->x_ref
[0], rotg
->vec
);
3243 /* This is just to ensure that it still works if all the atoms of the
3244 * reference structure are situated in a plane perpendicular to the rotation
3247 *lastindex
= rotg
->nat
-1;
3249 /* Loop over all atoms of the reference group,
3250 * project them on the rotation vector to find the extremes */
3251 for (i
= 0; i
< rotg
->nat
; i
++)
3253 xcproj
= iprod(rotg
->x_ref
[i
], rotg
->vec
);
3254 if (xcproj
< minproj
)
3259 if (xcproj
> maxproj
)
3268 /* Allocate memory for the slabs */
3269 static void allocate_slabs(
3275 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3279 erg
= rotg
->enfrotgrp
;
3281 /* More slabs than are defined for the reference are never needed */
3282 nslabs
= erg
->slab_last_ref
- erg
->slab_first_ref
+ 1;
3284 /* Remember how many we allocated */
3285 erg
->nslabs_alloc
= nslabs
;
3287 if ( (nullptr != fplog
) && bVerbose
)
3289 fprintf(fplog
, "%s allocating memory to store data for %d slabs (rotation group %d).\n",
3292 snew(erg
->slab_center
, nslabs
);
3293 snew(erg
->slab_center_ref
, nslabs
);
3294 snew(erg
->slab_weights
, nslabs
);
3295 snew(erg
->slab_torque_v
, nslabs
);
3296 snew(erg
->slab_data
, nslabs
);
3297 snew(erg
->gn_atom
, nslabs
);
3298 snew(erg
->gn_slabind
, nslabs
);
3299 snew(erg
->slab_innersumvec
, nslabs
);
3300 for (i
= 0; i
< nslabs
; i
++)
3302 snew(erg
->slab_data
[i
].x
, rotg
->nat
);
3303 snew(erg
->slab_data
[i
].ref
, rotg
->nat
);
3304 snew(erg
->slab_data
[i
].weight
, rotg
->nat
);
3306 snew(erg
->xc_ref_sorted
, rotg
->nat
);
3307 snew(erg
->xc_sortind
, rotg
->nat
);
3308 snew(erg
->firstatom
, nslabs
);
3309 snew(erg
->lastatom
, nslabs
);
3313 /* From the extreme positions of the reference group, determine the first
3314 * and last slab of the reference. We can never have more slabs in the real
3315 * simulation than calculated here for the reference.
3317 static void get_firstlast_slab_ref(t_rotgrp
*rotg
, real mc
[], int ref_firstindex
, int ref_lastindex
)
3319 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3324 erg
= rotg
->enfrotgrp
;
3325 first
= get_first_slab(rotg
, erg
->max_beta
, rotg
->x_ref
[ref_firstindex
]);
3326 last
= get_last_slab( rotg
, erg
->max_beta
, rotg
->x_ref
[ref_lastindex
]);
3328 while (get_slab_weight(first
, rotg
, rotg
->x_ref
, mc
, &dummy
) > WEIGHT_MIN
)
3332 erg
->slab_first_ref
= first
+1;
3333 while (get_slab_weight(last
, rotg
, rotg
->x_ref
, mc
, &dummy
) > WEIGHT_MIN
)
3337 erg
->slab_last_ref
= last
-1;
3341 /* Special version of copy_rvec:
3342 * During the copy procedure of xcurr to b, the correct PBC image is chosen
3343 * such that the copied vector ends up near its reference position xref */
3344 static gmx_inline
void copy_correct_pbc_image(
3345 const rvec xcurr
, /* copy vector xcurr ... */
3346 rvec b
, /* ... to b ... */
3347 const rvec xref
, /* choosing the PBC image such that b ends up near xref */
3356 /* Shortest PBC distance between the atom and its reference */
3357 rvec_sub(xcurr
, xref
, dx
);
3359 /* Determine the shift for this atom */
3361 for (m
= npbcdim
-1; m
>= 0; m
--)
3363 while (dx
[m
] < -0.5*box
[m
][m
])
3365 for (d
= 0; d
< DIM
; d
++)
3371 while (dx
[m
] >= 0.5*box
[m
][m
])
3373 for (d
= 0; d
< DIM
; d
++)
3381 /* Apply the shift to the position */
3382 copy_rvec(xcurr
, b
);
3383 shift_single_coord(box
, b
, shift
);
3387 static void init_rot_group(FILE *fplog
, t_commrec
*cr
, int g
, t_rotgrp
*rotg
,
3388 rvec
*x
, gmx_mtop_t
*mtop
, gmx_bool bVerbose
, FILE *out_slabs
, matrix box
,
3389 t_inputrec
*ir
, gmx_bool bOutputCenters
)
3392 rvec coord
, xref
, *xdum
;
3393 gmx_bool bFlex
, bColl
;
3394 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3395 int ref_firstindex
, ref_lastindex
;
3396 real mass
, totalmass
;
3401 /* Do we have a flexible axis? */
3402 bFlex
= ISFLEX(rotg
);
3403 /* Do we use a global set of coordinates? */
3404 bColl
= ISCOLL(rotg
);
3406 erg
= rotg
->enfrotgrp
;
3408 /* Allocate space for collective coordinates if needed */
3411 snew(erg
->xc
, rotg
->nat
);
3412 snew(erg
->xc_shifts
, rotg
->nat
);
3413 snew(erg
->xc_eshifts
, rotg
->nat
);
3414 snew(erg
->xc_old
, rotg
->nat
);
3416 if (rotg
->eFittype
== erotgFitNORM
)
3418 snew(erg
->xc_ref_length
, rotg
->nat
); /* in case fit type NORM is chosen */
3419 snew(erg
->xc_norm
, rotg
->nat
);
3424 snew(erg
->xr_loc
, rotg
->nat
);
3425 snew(erg
->x_loc_pbc
, rotg
->nat
);
3428 snew(erg
->f_rot_loc
, rotg
->nat
);
3429 snew(erg
->xc_ref_ind
, rotg
->nat
);
3431 /* Make space for the calculation of the potential at other angles (used
3432 * for fitting only) */
3433 if (erotgFitPOT
== rotg
->eFittype
)
3435 snew(erg
->PotAngleFit
, 1);
3436 snew(erg
->PotAngleFit
->degangle
, rotg
->PotAngle_nstep
);
3437 snew(erg
->PotAngleFit
->V
, rotg
->PotAngle_nstep
);
3438 snew(erg
->PotAngleFit
->rotmat
, rotg
->PotAngle_nstep
);
3440 /* Get the set of angles around the reference angle */
3441 start
= -0.5 * (rotg
->PotAngle_nstep
- 1)*rotg
->PotAngle_step
;
3442 for (i
= 0; i
< rotg
->PotAngle_nstep
; i
++)
3444 erg
->PotAngleFit
->degangle
[i
] = start
+ i
*rotg
->PotAngle_step
;
3449 erg
->PotAngleFit
= nullptr;
3452 /* xc_ref_ind needs to be set to identity in the serial case */
3455 for (i
= 0; i
< rotg
->nat
; i
++)
3457 erg
->xc_ref_ind
[i
] = i
;
3461 /* Copy the masses so that the center can be determined. For all types of
3462 * enforced rotation, we store the masses in the erg->mc array. */
3463 snew(erg
->mc
, rotg
->nat
);
3466 snew(erg
->mc_sorted
, rotg
->nat
);
3470 snew(erg
->m_loc
, rotg
->nat
);
3474 for (i
= 0; i
< rotg
->nat
; i
++)
3478 mass
= mtopGetAtomMass(mtop
, rotg
->ind
[i
], &molb
);
3487 erg
->invmass
= 1.0/totalmass
;
3489 /* Set xc_ref_center for any rotation potential */
3490 if ((rotg
->eType
== erotgISO
) || (rotg
->eType
== erotgPM
) || (rotg
->eType
== erotgRM
) || (rotg
->eType
== erotgRM2
))
3492 /* Set the pivot point for the fixed, stationary-axis potentials. This
3493 * won't change during the simulation */
3494 copy_rvec(rotg
->pivot
, erg
->xc_ref_center
);
3495 copy_rvec(rotg
->pivot
, erg
->xc_center
);
3499 /* Center of the reference positions */
3500 get_center(rotg
->x_ref
, erg
->mc
, rotg
->nat
, erg
->xc_ref_center
);
3502 /* Center of the actual positions */
3505 snew(xdum
, rotg
->nat
);
3506 for (i
= 0; i
< rotg
->nat
; i
++)
3509 copy_rvec(x
[ii
], xdum
[i
]);
3511 get_center(xdum
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
3517 gmx_bcast(sizeof(erg
->xc_center
), erg
->xc_center
, cr
);
3524 /* Save the original (whole) set of positions in xc_old such that at later
3525 * steps the rotation group can always be made whole again. If the simulation is
3526 * restarted, we compute the starting reference positions (given the time)
3527 * and assume that the correct PBC image of each position is the one nearest
3528 * to the current reference */
3531 /* Calculate the rotation matrix for this angle: */
3532 t_start
= ir
->init_t
+ ir
->init_step
*ir
->delta_t
;
3533 erg
->degangle
= rotg
->rate
* t_start
;
3534 calc_rotmat(rotg
->vec
, erg
->degangle
, erg
->rotmat
);
3536 for (i
= 0; i
< rotg
->nat
; i
++)
3540 /* Subtract pivot, rotate, and add pivot again. This will yield the
3541 * reference position for time t */
3542 rvec_sub(rotg
->x_ref
[i
], erg
->xc_ref_center
, coord
);
3543 mvmul(erg
->rotmat
, coord
, xref
);
3544 rvec_inc(xref
, erg
->xc_ref_center
);
3546 copy_correct_pbc_image(x
[ii
], erg
->xc_old
[i
], xref
, box
, 3);
3552 gmx_bcast(rotg
->nat
*sizeof(erg
->xc_old
[0]), erg
->xc_old
, cr
);
3557 if ( (rotg
->eType
!= erotgFLEX
) && (rotg
->eType
!= erotgFLEX2
) )
3559 /* Put the reference positions into origin: */
3560 for (i
= 0; i
< rotg
->nat
; i
++)
3562 rvec_dec(rotg
->x_ref
[i
], erg
->xc_ref_center
);
3566 /* Enforced rotation with flexible axis */
3569 /* Calculate maximum beta value from minimum gaussian (performance opt.) */
3570 erg
->max_beta
= calc_beta_max(rotg
->min_gaussian
, rotg
->slab_dist
);
3572 /* Determine the smallest and largest coordinate with respect to the rotation vector */
3573 get_firstlast_atom_ref(rotg
, &ref_firstindex
, &ref_lastindex
);
3575 /* From the extreme positions of the reference group, determine the first
3576 * and last slab of the reference. */
3577 get_firstlast_slab_ref(rotg
, erg
->mc
, ref_firstindex
, ref_lastindex
);
3579 /* Allocate memory for the slabs */
3580 allocate_slabs(rotg
, fplog
, g
, bVerbose
);
3582 /* Flexible rotation: determine the reference centers for the rest of the simulation */
3583 erg
->slab_first
= erg
->slab_first_ref
;
3584 erg
->slab_last
= erg
->slab_last_ref
;
3585 get_slab_centers(rotg
, rotg
->x_ref
, erg
->mc
, g
, -1, out_slabs
, bOutputCenters
, TRUE
);
3587 /* Length of each x_rotref vector from center (needed if fit routine NORM is chosen): */
3588 if (rotg
->eFittype
== erotgFitNORM
)
3590 for (i
= 0; i
< rotg
->nat
; i
++)
3592 rvec_sub(rotg
->x_ref
[i
], erg
->xc_ref_center
, coord
);
3593 erg
->xc_ref_length
[i
] = norm(coord
);
3600 extern void dd_make_local_rotation_groups(gmx_domdec_t
*dd
, t_rot
*rot
)
3605 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3609 for (g
= 0; g
< rot
->ngrp
; g
++)
3611 rotg
= &rot
->grp
[g
];
3612 erg
= rotg
->enfrotgrp
;
3615 dd_make_local_group_indices(ga2la
, rotg
->nat
, rotg
->ind
,
3616 &erg
->nat_loc
, &erg
->ind_loc
, &erg
->nalloc_loc
, erg
->xc_ref_ind
);
3621 /* Calculate the size of the MPI buffer needed in reduce_output() */
3622 static int calc_mpi_bufsize(t_rot
*rot
)
3625 int count_group
, count_total
;
3627 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3631 for (g
= 0; g
< rot
->ngrp
; g
++)
3633 rotg
= &rot
->grp
[g
];
3634 erg
= rotg
->enfrotgrp
;
3636 /* Count the items that are transferred for this group: */
3637 count_group
= 4; /* V, torque, angle, weight */
3639 /* Add the maximum number of slabs for flexible groups */
3642 count_group
+= erg
->slab_last_ref
- erg
->slab_first_ref
+ 1;
3645 /* Add space for the potentials at different angles: */
3646 if (erotgFitPOT
== rotg
->eFittype
)
3648 count_group
+= rotg
->PotAngle_nstep
;
3651 /* Add to the total number: */
3652 count_total
+= count_group
;
3659 extern void init_rot(FILE *fplog
, t_inputrec
*ir
, int nfile
, const t_filenm fnm
[],
3660 t_commrec
*cr
, rvec
*x
, matrix box
, gmx_mtop_t
*mtop
, const gmx_output_env_t
*oenv
,
3661 gmx_bool bVerbose
, unsigned long Flags
)
3666 int nat_max
= 0; /* Size of biggest rotation group */
3667 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
3668 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3669 rvec
*x_pbc
= nullptr; /* Space for the pbc-correct atom positions */
3672 if (MASTER(cr
) && bVerbose
)
3674 fprintf(stdout
, "%s Initializing ...\n", RotStr
);
3678 snew(rot
->enfrot
, 1);
3682 /* When appending, skip first output to avoid duplicate entries in the data files */
3683 if (er
->Flags
& MD_APPENDFILES
)
3692 if (MASTER(cr
) && er
->bOut
)
3694 please_cite(fplog
, "Kutzner2011");
3697 /* Output every step for reruns */
3698 if (er
->Flags
& MD_RERUN
)
3700 if (nullptr != fplog
)
3702 fprintf(fplog
, "%s rerun - will write rotation output every available step.\n", RotStr
);
3708 er
->out_slabs
= nullptr;
3709 if (MASTER(cr
) && HaveFlexibleGroups(rot
) )
3711 er
->out_slabs
= open_slab_out(opt2fn("-rs", nfile
, fnm
), rot
);
3716 /* Remove pbc, make molecule whole.
3717 * When ir->bContinuation=TRUE this has already been done, but ok. */
3718 snew(x_pbc
, mtop
->natoms
);
3719 copy_rvecn(x
, x_pbc
, 0, mtop
->natoms
);
3720 do_pbc_first_mtop(nullptr, ir
->ePBC
, box
, mtop
, x_pbc
);
3721 /* All molecules will be whole now, but not necessarily in the home box.
3722 * Additionally, if a rotation group consists of more than one molecule
3723 * (e.g. two strands of DNA), each one of them can end up in a different
3724 * periodic box. This is taken care of in init_rot_group. */
3727 for (g
= 0; g
< rot
->ngrp
; g
++)
3729 rotg
= &rot
->grp
[g
];
3731 if (nullptr != fplog
)
3733 fprintf(fplog
, "%s group %d type '%s'\n", RotStr
, g
, erotg_names
[rotg
->eType
]);
3738 /* Allocate space for the rotation group's data: */
3739 snew(rotg
->enfrotgrp
, 1);
3740 erg
= rotg
->enfrotgrp
;
3742 nat_max
= std::max(nat_max
, rotg
->nat
);
3747 erg
->nalloc_loc
= 0;
3748 erg
->ind_loc
= nullptr;
3752 erg
->nat_loc
= rotg
->nat
;
3753 erg
->ind_loc
= rotg
->ind
;
3755 init_rot_group(fplog
, cr
, g
, rotg
, x_pbc
, mtop
, bVerbose
, er
->out_slabs
, box
, ir
,
3756 !(er
->Flags
& MD_APPENDFILES
) ); /* Do not output the reference centers
3757 * again if we are appending */
3761 /* Allocate space for enforced rotation buffer variables */
3762 er
->bufsize
= nat_max
;
3763 snew(er
->data
, nat_max
);
3764 snew(er
->xbuf
, nat_max
);
3765 snew(er
->mbuf
, nat_max
);
3767 /* Buffers for MPI reducing torques, angles, weights (for each group), and V */
3770 er
->mpi_bufsize
= calc_mpi_bufsize(rot
) + 100; /* larger to catch errors */
3771 snew(er
->mpi_inbuf
, er
->mpi_bufsize
);
3772 snew(er
->mpi_outbuf
, er
->mpi_bufsize
);
3776 er
->mpi_bufsize
= 0;
3777 er
->mpi_inbuf
= nullptr;
3778 er
->mpi_outbuf
= nullptr;
3781 /* Only do I/O on the MASTER */
3782 er
->out_angles
= nullptr;
3783 er
->out_rot
= nullptr;
3784 er
->out_torque
= nullptr;
3787 er
->out_rot
= open_rot_out(opt2fn("-ro", nfile
, fnm
), rot
, oenv
);
3789 if (rot
->nstsout
> 0)
3791 if (HaveFlexibleGroups(rot
) || HavePotFitGroups(rot
) )
3793 er
->out_angles
= open_angles_out(opt2fn("-ra", nfile
, fnm
), rot
);
3795 if (HaveFlexibleGroups(rot
) )
3797 er
->out_torque
= open_torque_out(opt2fn("-rt", nfile
, fnm
), rot
);
3806 extern void finish_rot(t_rot
*rot
)
3808 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
3814 gmx_fio_fclose(er
->out_rot
);
3818 gmx_fio_fclose(er
->out_slabs
);
3822 gmx_fio_fclose(er
->out_angles
);
3826 gmx_fio_fclose(er
->out_torque
);
3831 /* Rotate the local reference positions and store them in
3832 * erg->xr_loc[0...(nat_loc-1)]
3834 * Note that we already subtracted u or y_c from the reference positions
3835 * in init_rot_group().
3837 static void rotate_local_reference(t_rotgrp
*rotg
)
3839 gmx_enfrotgrp_t erg
;
3843 erg
= rotg
->enfrotgrp
;
3845 for (i
= 0; i
< erg
->nat_loc
; i
++)
3847 /* Index of this rotation group atom with respect to the whole rotation group */
3848 ii
= erg
->xc_ref_ind
[i
];
3850 mvmul(erg
->rotmat
, rotg
->x_ref
[ii
], erg
->xr_loc
[i
]);
3855 /* Select the PBC representation for each local x position and store that
3856 * for later usage. We assume the right PBC image of an x is the one nearest to
3857 * its rotated reference */
3858 static void choose_pbc_image(rvec x
[], t_rotgrp
*rotg
, matrix box
, int npbcdim
)
3861 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3865 erg
= rotg
->enfrotgrp
;
3867 for (i
= 0; i
< erg
->nat_loc
; i
++)
3869 /* Index of a rotation group atom */
3870 ii
= erg
->ind_loc
[i
];
3872 /* Get the correctly rotated reference position. The pivot was already
3873 * subtracted in init_rot_group() from the reference positions. Also,
3874 * the reference positions have already been rotated in
3875 * rotate_local_reference(). For the current reference position we thus
3876 * only need to add the pivot again. */
3877 copy_rvec(erg
->xr_loc
[i
], xref
);
3878 rvec_inc(xref
, erg
->xc_ref_center
);
3880 copy_correct_pbc_image(x
[ii
], erg
->x_loc_pbc
[i
], xref
, box
, npbcdim
);
3885 extern void do_rotation(
3897 gmx_bool outstep_slab
, outstep_rot
;
3899 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
3900 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3902 t_gmx_potfit
*fit
= nullptr; /* For fit type 'potential' determine the fit
3903 angle via the potential minimum */
3912 /* When to output in main rotation output file */
3913 outstep_rot
= do_per_step(step
, rot
->nstrout
) && er
->bOut
;
3914 /* When to output per-slab data */
3915 outstep_slab
= do_per_step(step
, rot
->nstsout
) && er
->bOut
;
3917 /* Output time into rotation output file */
3918 if (outstep_rot
&& MASTER(cr
))
3920 fprintf(er
->out_rot
, "%12.3e", t
);
3923 /**************************************************************************/
3924 /* First do ALL the communication! */
3925 for (g
= 0; g
< rot
->ngrp
; g
++)
3927 rotg
= &rot
->grp
[g
];
3928 erg
= rotg
->enfrotgrp
;
3930 /* Do we use a collective (global) set of coordinates? */
3931 bColl
= ISCOLL(rotg
);
3933 /* Calculate the rotation matrix for this angle: */
3934 erg
->degangle
= rotg
->rate
* t
;
3935 calc_rotmat(rotg
->vec
, erg
->degangle
, erg
->rotmat
);
3939 /* Transfer the rotation group's positions such that every node has
3940 * all of them. Every node contributes its local positions x and stores
3941 * it in the collective erg->xc array. */
3942 communicate_group_positions(cr
, erg
->xc
, erg
->xc_shifts
, erg
->xc_eshifts
, bNS
,
3943 x
, rotg
->nat
, erg
->nat_loc
, erg
->ind_loc
, erg
->xc_ref_ind
, erg
->xc_old
, box
);
3947 /* Fill the local masses array;
3948 * this array changes in DD/neighborsearching steps */
3951 for (i
= 0; i
< erg
->nat_loc
; i
++)
3953 /* Index of local atom w.r.t. the collective rotation group */
3954 ii
= erg
->xc_ref_ind
[i
];
3955 erg
->m_loc
[i
] = erg
->mc
[ii
];
3959 /* Calculate Omega*(y_i-y_c) for the local positions */
3960 rotate_local_reference(rotg
);
3962 /* Choose the nearest PBC images of the group atoms with respect
3963 * to the rotated reference positions */
3964 choose_pbc_image(x
, rotg
, box
, 3);
3966 /* Get the center of the rotation group */
3967 if ( (rotg
->eType
== erotgISOPF
) || (rotg
->eType
== erotgPMPF
) )
3969 get_center_comm(cr
, erg
->x_loc_pbc
, erg
->m_loc
, erg
->nat_loc
, rotg
->nat
, erg
->xc_center
);
3973 } /* End of loop over rotation groups */
3975 /**************************************************************************/
3976 /* Done communicating, we can start to count cycles for the load balancing now ... */
3977 if (DOMAINDECOMP(cr
))
3979 ddReopenBalanceRegionCpu(cr
->dd
);
3986 for (g
= 0; g
< rot
->ngrp
; g
++)
3988 rotg
= &rot
->grp
[g
];
3989 erg
= rotg
->enfrotgrp
;
3991 if (outstep_rot
&& MASTER(cr
))
3993 fprintf(er
->out_rot
, "%12.4f", erg
->degangle
);
3996 /* Calculate angles and rotation matrices for potential fitting: */
3997 if ( (outstep_rot
|| outstep_slab
) && (erotgFitPOT
== rotg
->eFittype
) )
3999 fit
= erg
->PotAngleFit
;
4000 for (i
= 0; i
< rotg
->PotAngle_nstep
; i
++)
4002 calc_rotmat(rotg
->vec
, erg
->degangle
+ fit
->degangle
[i
], fit
->rotmat
[i
]);
4004 /* Clear value from last step */
4005 erg
->PotAngleFit
->V
[i
] = 0.0;
4009 /* Clear values from last time step */
4011 erg
->torque_v
= 0.0;
4013 erg
->weight_v
= 0.0;
4015 switch (rotg
->eType
)
4021 do_fixed(rotg
, outstep_rot
, outstep_slab
);
4024 do_radial_motion(rotg
, outstep_rot
, outstep_slab
);
4027 do_radial_motion_pf(rotg
, x
, box
, outstep_rot
, outstep_slab
);
4031 do_radial_motion2(rotg
, x
, box
, outstep_rot
, outstep_slab
);
4035 /* Subtract the center of the rotation group from the collective positions array
4036 * Also store the center in erg->xc_center since it needs to be subtracted
4037 * in the low level routines from the local coordinates as well */
4038 get_center(erg
->xc
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
4039 svmul(-1.0, erg
->xc_center
, transvec
);
4040 translate_x(erg
->xc
, rotg
->nat
, transvec
);
4041 do_flexible(MASTER(cr
), er
, rotg
, g
, x
, box
, t
, outstep_rot
, outstep_slab
);
4045 /* Do NOT subtract the center of mass in the low level routines! */
4046 clear_rvec(erg
->xc_center
);
4047 do_flexible(MASTER(cr
), er
, rotg
, g
, x
, box
, t
, outstep_rot
, outstep_slab
);
4050 gmx_fatal(FARGS
, "No such rotation potential.");
4058 fprintf(stderr
, "%s calculation (step %d) took %g seconds.\n", RotStr
, step
, MPI_Wtime()-t0
);