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