3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2008, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
32 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
42 #include "gmx_wallcycle.h"
51 #include "mtop_util.h"
55 #include "gmx_ga2la.h"
58 #include "mpelogging.h"
59 #include "groupcoord.h"
60 #include "pull_rotation.h"
64 static char *RotStr
= {"Enforced rotation:"};
67 /* Set the minimum weight for the determination of the slab centers */
68 #define WEIGHT_MIN (10*GMX_FLOAT_MIN)
70 /* Helper structure for sorting positions along rotation vector */
72 real xcproj
; /* Projection of xc on the rotation vector */
73 int ind
; /* Index of xc */
75 rvec x
; /* Position */
76 rvec x_ref
; /* Reference position */
80 /* Enforced rotation / flexible: determine the angle of each slab */
81 typedef struct gmx_slabdata
83 int nat
; /* Number of atoms belonging to this slab */
84 rvec
*x
; /* The positions belonging to this slab. In
85 general, this should be all positions of the
86 whole rotation group, but we leave those away
87 that have a small enough weight */
88 rvec
*ref
; /* Same for reference */
89 real
*weight
; /* The weight for each atom */
93 /* Enforced rotation data for all groups */
94 typedef struct gmx_enfrot
96 FILE *out_rot
; /* Output file for rotation data */
97 FILE *out_torque
; /* Output file for torque data */
98 FILE *out_angles
; /* Output file for slab angles for flexible type */
99 FILE *out_slabs
; /* Output file for slab centers */
100 int bufsize
; /* Allocation size of buf */
101 rvec
*xbuf
; /* Coordinate buffer variable for sorting */
102 real
*mbuf
; /* Masses buffer variable for sorting */
103 sort_along_vec_t
*data
; /* Buffer variable needed for position sorting */
104 real
*mpi_inbuf
; /* MPI buffer */
105 real
*mpi_outbuf
; /* MPI buffer */
106 int mpi_bufsize
; /* Allocation size of in & outbuf */
107 real Vrot
; /* (Local) part of the enf. rotation potential */
111 /* Global enforced rotation data for a single rotation group */
112 typedef struct gmx_enfrotgrp
114 real degangle
; /* Rotation angle in degree */
115 matrix rotmat
; /* Rotation matrix */
116 atom_id
*ind_loc
; /* Local rotation indices */
117 int nat_loc
; /* Number of local group atoms */
118 int nalloc_loc
; /* Allocation size for ind_loc and weight_loc */
120 real V
; /* Rotation potential for this rotation group */
121 rvec
*f_rot_loc
; /* Array to store the forces on the local atoms
122 resulting from enforced rotation potential */
124 /* Collective coordinates for the whole rotation group */
125 real
*xc_ref_length
; /* Length of each x_rotref vector after x_rotref
126 has been put into origin */
127 int *xc_ref_ind
; /* Position of each local atom in the collective
129 rvec xc_center
; /* Center of the rotation group positions, may
131 rvec xc_ref_center
; /* dito, for the reference positions */
132 rvec
*xc
; /* Current (collective) positions */
133 ivec
*xc_shifts
; /* Current (collective) shifts */
134 ivec
*xc_eshifts
; /* Extra shifts since last DD step */
135 rvec
*xc_old
; /* Old (collective) positions */
136 rvec
*xc_norm
; /* Normalized form of the current positions */
137 rvec
*xc_ref_sorted
; /* Reference positions (sorted in the same order
138 as xc when sorted) */
139 int *xc_sortind
; /* Where is a position found after sorting? */
140 real
*mc
; /* Collective masses */
142 real invmass
; /* one over the total mass of the rotation group */
143 /* Fixed rotation only */
144 rvec
*xr_loc
; /* Local reference coords, correctly rotated */
145 rvec
*x_loc_pbc
; /* Local current coords, correct PBC image */
146 real
*m_loc
; /* Masses of the current local atoms */
147 real fix_torque_v
; /* Torque in the direction of rotation vector */
150 /* Flexible rotation only */
151 int nslabs_alloc
; /* For this many slabs memory is allocated */
152 int slab_first
; /* Lowermost slab for that the calculation needs
153 to be performed at a given time step */
154 int slab_last
; /* Uppermost slab ... */
155 int slab_first_ref
; /* First slab for which reference COG is stored */
156 int slab_last_ref
; /* Last ... */
157 int slab_buffer
; /* Slab buffer region around reference slabs */
158 int *firstatom
; /* First relevant atom for a slab */
159 int *lastatom
; /* Last relevant atom for a slab */
160 rvec
*slab_center
; /* Gaussian-weighted slab center (COG) */
161 rvec
*slab_center_ref
; /* Gaussian-weighted slab COG for the
162 reference positions */
163 real
*slab_weights
; /* Sum of gaussian weights in a slab */
164 real
*slab_torque_v
; /* Torque T = r x f for each slab. */
165 /* torque_v = m.v = angular momentum in the
167 real max_beta
; /* min_gaussian from inputrec->rotgrp is the
168 minimum value the gaussian must have so that
169 the force is actually evaluated max_beta is
170 just another way to put it */
171 real
*gn_atom
; /* Precalculated gaussians for a single atom */
172 int *gn_slabind
; /* Tells to which slab each precalculated gaussian
174 rvec
*slab_innersumvec
;/* Inner sum of the flexible2 potential per slab;
175 this is precalculated for optimization reasons */
176 t_gmx_slabdata
*slab_data
; /* Holds atom positions and gaussian weights
177 of atoms belonging to a slab */
181 static double** allocate_square_matrix(int dim
)
195 static void free_square_matrix(double** mat
, int dim
)
200 for (i
=0; i
<dim
; i
++)
206 /* Output rotation energy and torque for each rotation group */
207 static void reduce_output(t_commrec
*cr
, t_rot
*rot
, real t
)
209 int g
,i
,islab
,nslabs
=0;
210 int count
; /* MPI element counter */
212 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
213 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
219 /* Fill the MPI buffer with stuff to reduce: */
223 for (g
=0; g
< rot
->ngrp
; g
++)
227 nslabs
= erg
->slab_last
- erg
->slab_first
+ 1;
228 er
->mpi_inbuf
[count
++] = erg
->V
;
239 er
->mpi_inbuf
[count
++] = erg
->fix_torque_v
;
240 er
->mpi_inbuf
[count
++] = erg
->fix_angles_v
;
241 er
->mpi_inbuf
[count
++] = erg
->fix_weight_v
;
247 /* (Re-)allocate memory for MPI buffer: */
248 if (er
->mpi_bufsize
< count
+nslabs
)
250 er
->mpi_bufsize
= count
+nslabs
;
251 srenew(er
->mpi_inbuf
, er
->mpi_bufsize
);
252 srenew(er
->mpi_outbuf
, er
->mpi_bufsize
);
254 for (i
=0; i
<nslabs
; i
++)
255 er
->mpi_inbuf
[count
++] = erg
->slab_torque_v
[i
];
262 MPI_Reduce(er
->mpi_inbuf
, er
->mpi_outbuf
, count
, GMX_MPI_REAL
, MPI_SUM
, MASTERRANK(cr
), cr
->mpi_comm_mygroup
);
264 /* Copy back the reduced data from the buffer on the master */
268 for (g
=0; g
< rot
->ngrp
; g
++)
272 nslabs
= erg
->slab_last
- erg
->slab_first
+ 1;
273 erg
->V
= er
->mpi_outbuf
[count
++];
284 erg
->fix_torque_v
= er
->mpi_outbuf
[count
++];
285 erg
->fix_angles_v
= er
->mpi_outbuf
[count
++];
286 erg
->fix_weight_v
= er
->mpi_outbuf
[count
++];
292 for (i
=0; i
<nslabs
; i
++)
293 erg
->slab_torque_v
[i
] = er
->mpi_outbuf
[count
++];
305 /* Av. angle and total torque for each rotation group */
306 for (g
=0; g
< rot
->ngrp
; g
++)
309 bFlex
= ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
310 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) );
314 /* Output to main rotation log file: */
317 fprintf(er
->out_rot
, "%12.4f%12.3e",
318 (erg
->fix_angles_v
/erg
->fix_weight_v
)*180.0*M_1_PI
,
321 fprintf(er
->out_rot
, "%12.3e", erg
->V
);
323 /* Output to torque log file: */
326 fprintf(er
->out_torque
, "%12.3e%6d", t
, g
);
327 for (i
=erg
->slab_first
; i
<=erg
->slab_last
; i
++)
329 islab
= i
- erg
->slab_first
; /* slab index */
330 /* Only output if enough weight is in slab */
331 if (erg
->slab_weights
[islab
] > rotg
->min_gaussian
)
332 fprintf(er
->out_torque
, "%6d%12.3e", i
, erg
->slab_torque_v
[islab
]);
334 fprintf(er
->out_torque
, "\n");
337 fprintf(er
->out_rot
, "\n");
342 /* Add the forces from enforced rotation potential to the local forces.
343 * Should be called after the SR forces have been evaluated */
344 extern real
add_rot_forces(t_rot
*rot
, rvec f
[], t_commrec
*cr
, int step
, real t
)
348 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
349 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
354 GMX_MPE_LOG(ev_add_rot_forces_start
);
356 /* Reduce energy,torque, angles etc. to get the sum values (per rotation group)
357 * on the master and output these values to file. */
358 if (do_per_step(step
, rot
->nsttout
))
359 reduce_output(cr
, rot
, t
);
361 /* Total rotation potential is the sum over all rotation groups */
364 /* Loop over enforced rotation groups (usually 1, though)
365 * Apply the forces from rotation potentials */
366 for (g
=0; g
<rot
->ngrp
; g
++)
371 for (l
=0; l
<erg
->nat_loc
; l
++)
373 /* Get the right index of the local force */
374 ii
= erg
->ind_loc
[l
];
376 rvec_inc(f
[ii
],erg
->f_rot_loc
[l
]);
380 GMX_MPE_LOG(ev_add_rot_forces_finish
);
382 return (MASTER(cr
)? er
->Vrot
: 0.0);
386 /* Calculate the maximum beta that leads to a gaussian larger min_gaussian,
387 * also does some checks
389 static double calc_beta_max(real min_gaussian
, real slab_dist
)
391 const double norm
= 0.5698457353514458216; /* = 1/1.7548609 */
396 /* Actually the next two checks are already made in grompp */
398 gmx_fatal(FARGS
, "Slab distance of flexible rotation groups must be >=0 !");
399 if (min_gaussian
<= 0)
400 gmx_fatal(FARGS
, "Cutoff value for Gaussian must be > 0. (You requested %f)");
402 /* Define the sigma value */
403 sigma
= 0.7*slab_dist
;
405 /* Calculate the argument for the logarithm and check that the log() result is negative or 0 */
406 arg
= min_gaussian
/norm
;
408 gmx_fatal(FARGS
, "min_gaussian of flexible rotation groups must be <%g", norm
);
410 return sqrt(-2.0*sigma
*sigma
*log(min_gaussian
/norm
));
414 static inline real
calc_beta(rvec curr_x
, t_rotgrp
*rotg
, int n
)
416 return iprod(curr_x
, rotg
->vec
) - rotg
->slab_dist
* n
;
420 static inline real
gaussian_weight(rvec curr_x
, t_rotgrp
*rotg
, int n
)
422 /* norm is chosen such that the sum of the gaussians
423 * over the slabs is approximately 1.0 everywhere */
424 /* a previously used value was norm = 0.5698457353514458216 = 1/1.7548609 */
425 const real norm
= 0.569917543430618; /* = 1/1.7546397922417 */
429 /* Define the sigma value */
430 sigma
= 0.7*rotg
->slab_dist
;
431 /* Calculate the Gaussian value of slab n for position curr_x */
432 return norm
* exp( -0.5 * sqr( calc_beta(curr_x
, rotg
, n
)/sigma
) );
436 /* Returns the weight in a single slab, also calculates the Gaussian- and mass-
437 * weighted sum of positions for that slab */
438 static real
get_slab_weight(int j
, t_rotgrp
*rotg
, rvec xc
[], real mc
[], rvec
*x_weighted_sum
)
440 rvec curr_x
; /* The position of an atom */
441 rvec curr_x_weighted
; /* The gaussian-weighted position */
442 real gaussian
; /* A single gaussian weight */
443 real wgauss
; /* gaussian times current mass */
444 real slabweight
= 0.0; /* The sum of weights in the slab */
446 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
450 clear_rvec(*x_weighted_sum
);
453 islab
= j
- erg
->slab_first
;
455 /* Loop over all atoms in the rotation group */
456 for (i
=0; i
<rotg
->nat
; i
++)
458 copy_rvec(xc
[i
], curr_x
);
459 gaussian
= gaussian_weight(curr_x
, rotg
, j
);
460 wgauss
= gaussian
* mc
[i
];
461 svmul(wgauss
, curr_x
, curr_x_weighted
);
462 rvec_add(*x_weighted_sum
, curr_x_weighted
, *x_weighted_sum
);
463 slabweight
+= wgauss
;
464 } /* END of loop over rotation group atoms */
470 static void get_slab_centers(
471 t_rotgrp
*rotg
, /* The rotation group information */
472 rvec
*xc
, /* The rotation group positions; will
473 typically be enfrotgrp->xc, but at first call
474 it is enfrotgrp->xc_ref */
475 real
*mc
, /* The masses of the rotation group atoms */
476 t_commrec
*cr
, /* Communication record */
477 int g
, /* The number of the rotation group */
478 real time
, /* Used for output only */
479 FILE *out_slabs
, /* For outputting center per slab information */
480 gmx_bool bOutStep
, /* Is this an output step? */
481 gmx_bool bReference
) /* If this routine is called from
482 init_rot_group we need to store
483 the reference slab centers */
486 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
491 /* Loop over slabs */
492 for (j
= erg
->slab_first
; j
<= erg
->slab_last
; j
++)
494 islab
= j
- erg
->slab_first
;
495 erg
->slab_weights
[islab
] = get_slab_weight(j
, rotg
, xc
, mc
, &erg
->slab_center
[islab
]);
497 /* We can do the calculations ONLY if there is weight in the slab! */
498 if (erg
->slab_weights
[islab
] > WEIGHT_MIN
)
500 svmul(1.0/erg
->slab_weights
[islab
], erg
->slab_center
[islab
], erg
->slab_center
[islab
]);
504 /* We need to check this here, since we divide through slab_weights
505 * in the flexible low-level routines! */
506 gmx_fatal(FARGS
, "Not enough weight in slab %d. Slab center cannot be determined!", j
);
509 /* At first time step: save the centers of the reference structure */
511 copy_rvec(erg
->slab_center
[islab
], erg
->slab_center_ref
[islab
]);
512 } /* END of loop over slabs */
514 /* Output on the master */
515 if (MASTER(cr
) && bOutStep
)
517 fprintf(out_slabs
, "%12.3e%6d", time
, g
);
518 for (j
= erg
->slab_first
; j
<= erg
->slab_last
; j
++)
520 islab
= j
- erg
->slab_first
;
521 fprintf(out_slabs
, "%6d%12.3e%12.3e%12.3e",
522 j
,erg
->slab_center
[islab
][XX
],erg
->slab_center
[islab
][YY
],erg
->slab_center
[islab
][ZZ
]);
524 fprintf(out_slabs
, "\n");
529 static void calc_rotmat(
531 real degangle
, /* Angle alpha of rotation at time t in degrees */
532 matrix rotmat
) /* Rotation matrix */
534 real radangle
; /* Rotation angle in radians */
535 real cosa
; /* cosine alpha */
536 real sina
; /* sine alpha */
537 real OMcosa
; /* 1 - cos(alpha) */
538 real dumxy
, dumxz
, dumyz
; /* save computations */
539 rvec rot_vec
; /* Rotate around rot_vec ... */
542 radangle
= degangle
* M_PI
/180.0;
543 copy_rvec(vec
, rot_vec
);
545 /* Precompute some variables: */
546 cosa
= cos(radangle
);
547 sina
= sin(radangle
);
549 dumxy
= rot_vec
[XX
]*rot_vec
[YY
]*OMcosa
;
550 dumxz
= rot_vec
[XX
]*rot_vec
[ZZ
]*OMcosa
;
551 dumyz
= rot_vec
[YY
]*rot_vec
[ZZ
]*OMcosa
;
553 /* Construct the rotation matrix for this rotation group: */
555 rotmat
[XX
][XX
] = cosa
+ rot_vec
[XX
]*rot_vec
[XX
]*OMcosa
;
556 rotmat
[YY
][XX
] = dumxy
+ rot_vec
[ZZ
]*sina
;
557 rotmat
[ZZ
][XX
] = dumxz
- rot_vec
[YY
]*sina
;
559 rotmat
[XX
][YY
] = dumxy
- rot_vec
[ZZ
]*sina
;
560 rotmat
[YY
][YY
] = cosa
+ rot_vec
[YY
]*rot_vec
[YY
]*OMcosa
;
561 rotmat
[ZZ
][YY
] = dumyz
+ rot_vec
[XX
]*sina
;
563 rotmat
[XX
][ZZ
] = dumxz
+ rot_vec
[YY
]*sina
;
564 rotmat
[YY
][ZZ
] = dumyz
- rot_vec
[XX
]*sina
;
565 rotmat
[ZZ
][ZZ
] = cosa
+ rot_vec
[ZZ
]*rot_vec
[ZZ
]*OMcosa
;
570 for (iii
=0; iii
<3; iii
++) {
571 for (jjj
=0; jjj
<3; jjj
++)
572 fprintf(stderr
, " %10.8f ", rotmat
[iii
][jjj
]);
573 fprintf(stderr
, "\n");
579 /* Calculates torque on the rotation axis tau = position x force */
580 static inline real
torque(
581 rvec rotvec
, /* rotation vector; MUST be normalized! */
582 rvec force
, /* force */
583 rvec x
, /* position of atom on which the force acts */
584 rvec pivot
) /* pivot point of rotation axis */
589 /* Subtract offset */
590 rvec_sub(x
,pivot
,vectmp
);
592 /* position x force */
593 cprod(vectmp
, force
, tau
);
595 /* Return the part of the torque which is parallel to the rotation vector */
596 return iprod(tau
, rotvec
);
600 /* Right-aligned output of value with standard width */
601 static void print_aligned(FILE *fp
, char *str
)
603 fprintf(fp
, "%12s", str
);
607 /* Right-aligned output of value with standard short width */
608 static void print_aligned_short(FILE *fp
, char *str
)
610 fprintf(fp
, "%6s", str
);
614 /* Right-aligned output of value with standard width */
615 static void print_aligned_group(FILE *fp
, char *str
, int g
)
620 sprintf(sbuf
, "%s%d", str
, g
);
621 fprintf(fp
, "%12s", sbuf
);
625 static FILE *open_output_file(const char *fn
, int steps
)
630 fp
= ffopen(fn
, "w");
632 fprintf(fp
, "# Output is written every %d time steps.\n\n", steps
);
638 /* Open output file for slab COG data. Call on master only */
639 static FILE *open_slab_out(t_rot
*rot
, const char *fn
)
646 for (g
=0; g
<rot
->ngrp
; g
++)
649 if ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
650 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) )
653 fp
= open_output_file(fn
, rot
->nsttout
);
654 fprintf(fp
, "# Rotation group %d (%s), slab distance %f nm\n", g
, erotg_names
[rotg
->eType
], rotg
->slab_dist
);
660 fprintf(fp
, "# The following columns will have the syntax: (COG = center of geometry, gaussian weighted)\n");
662 print_aligned_short(fp
, "t");
663 print_aligned_short(fp
, "grp");
664 print_aligned_short(fp
, "slab");
665 print_aligned(fp
, "COG-X");
666 print_aligned(fp
, "COG-Y");
667 print_aligned(fp
, "COG-Z");
668 print_aligned_short(fp
, "slab");
669 print_aligned(fp
, "COG-X");
670 print_aligned(fp
, "COG-Y");
671 print_aligned(fp
, "COG-Z");
672 print_aligned_short(fp
, "slab");
673 fprintf(fp
, " ...\n");
681 /* Open output file and print some general information about the rotation groups.
682 * Call on master only */
683 static FILE *open_rot_out(const char *fn
, t_rot
*rot
, const output_env_t oenv
,
689 const char **setname
;
691 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
695 if (Flags
& MD_APPENDFILES
)
697 fp
= gmx_fio_fopen(fn
,"a");
701 fp
= xvgropen(fn
, "Rotation angles and energy", "Time [ps]", "angles [degree] and energies [kJ/mol]", oenv
);
702 fprintf(fp
, "# The scalar tau is the torque [kJ/mol] in the direction of the rotation vector v.\n");
703 fprintf(fp
, "# To obtain the vectorial torque, multiply tau with the group's rot_vec.\n#\n");
705 for (g
=0; g
<rot
->ngrp
; g
++)
709 bFlex
= ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
710 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) );
713 fprintf(fp
, "# Rotation group %d, potential type '%s':\n" , g
, erotg_names
[rotg
->eType
]);
714 fprintf(fp
, "# rot_massw%d %s\n" , g
, yesno_names
[rotg
->bMassW
]);
715 fprintf(fp
, "# rot_vec%d %12.5e %12.5e %12.5e\n" , g
, rotg
->vec
[XX
], rotg
->vec
[YY
], rotg
->vec
[ZZ
]);
716 fprintf(fp
, "# rot_rate%d %12.5e degree/ps\n" , g
, rotg
->rate
);
717 fprintf(fp
, "# rot_k%d %12.5e kJ/(mol*nm^2)\n" , g
, rotg
->k
);
718 if ( rotg
->eType
==erotgISO
|| rotg
->eType
==erotgPM
|| rotg
->eType
==erotgRM
|| rotg
->eType
==erotgRM2
)
719 fprintf(fp
, "# rot_pivot%d %12.5e %12.5e %12.5e nm\n", g
, rotg
->pivot
[XX
], rotg
->pivot
[YY
], rotg
->pivot
[ZZ
]);
723 fprintf(fp
, "# rot_slab_distance%d %f nm\n", g
, rotg
->slab_dist
);
724 fprintf(fp
, "# rot_min_gaussian%d %12.5e\n", g
, rotg
->min_gaussian
);
727 /* Output the centers of the rotation groups for the pivot-free potentials */
728 if ((rotg
->eType
==erotgISOPF
) || (rotg
->eType
==erotgPMPF
) || (rotg
->eType
==erotgRMPF
) || (rotg
->eType
==erotgRM2PF
729 || (rotg
->eType
==erotgFLEXT
) || (rotg
->eType
==erotgFLEX2T
)) )
731 fprintf(fp
, "# %s of ref. grp. %d %12.5e %12.5e %12.5e\n",
732 rotg
->bMassW
? "COM":"COG", g
,
733 erg
->xc_ref_center
[XX
], erg
->xc_ref_center
[YY
], erg
->xc_ref_center
[ZZ
]);
735 fprintf(fp
, "# initial %s grp. %d %12.5e %12.5e %12.5e\n",
736 rotg
->bMassW
? "COM":"COG", g
,
737 erg
->xc_center
[XX
], erg
->xc_center
[YY
], erg
->xc_center
[ZZ
]);
740 if ( (rotg
->eType
== erotgRM2
) || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) )
742 fprintf(fp
, "# rot_eps%d %12.5e nm^2\n", g
, rotg
->eps
);
746 fprintf(fp
, "#\n# Legend for the following data columns:\n");
748 print_aligned_short(fp
, "t");
750 snew(setname
, 4*rot
->ngrp
);
752 for (g
=0; g
<rot
->ngrp
; g
++)
755 sprintf(buf
, "theta_ref%d [degree]", g
);
756 print_aligned_group(fp
, "theta_ref", g
);
757 setname
[nsets
] = strdup(buf
);
760 for (g
=0; g
<rot
->ngrp
; g
++)
763 bFlex
= ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
764 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) );
766 /* For flexible axis rotation we use RMSD fitting to determine the
767 * actual angle of the rotation group */
770 sprintf(buf
, "theta-av%d [degree]", g
);
771 print_aligned_group(fp
, "theta_av", g
);
772 setname
[nsets
] = strdup(buf
);
774 sprintf(buf
, "tau%d [kJ/mol]", g
);
775 print_aligned_group(fp
, "tau", g
);
776 setname
[nsets
] = strdup(buf
);
779 sprintf(buf
, "energy%d [kJ/mol]", g
);
780 print_aligned_group(fp
, "energy", g
);
781 setname
[nsets
] = strdup(buf
);
784 fprintf(fp
, "\n#\n");
787 xvgr_legend(fp
, nsets
, setname
, oenv
);
797 /* Call on master only */
798 static FILE *open_angles_out(t_rot
*rot
, const char *fn
)
805 /* Open output file and write some information about it's structure: */
806 fp
= open_output_file(fn
, rot
->nstrout
);
807 fprintf(fp
, "# All angles given in degrees, time in ps\n");
808 for (g
=0; g
<rot
->ngrp
; g
++)
811 if ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
812 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) )
814 fprintf(fp
, "# Rotation group %d (%s), slab distance %f nm, fit type %s\n",
815 g
, erotg_names
[rotg
->eType
], rotg
->slab_dist
, erotg_fitnames
[rotg
->eFittype
]);
818 fprintf(fp
, "# The following columns will have the syntax:\n");
820 print_aligned_short(fp
, "t");
821 print_aligned_short(fp
, "grp");
822 print_aligned(fp
, "theta_ref");
823 print_aligned(fp
, "theta_fit");
824 print_aligned_short(fp
, "slab");
825 print_aligned_short(fp
, "atoms");
826 print_aligned(fp
, "theta_fit");
827 print_aligned_short(fp
, "slab");
828 print_aligned_short(fp
, "atoms");
829 print_aligned(fp
, "theta_fit");
830 fprintf(fp
, " ...\n");
836 /* Open torque output file and write some information about it's structure.
837 * Call on master only */
838 static FILE *open_torque_out(t_rot
*rot
, const char *fn
)
845 fp
= open_output_file(fn
, rot
->nsttout
);
847 for (g
=0; g
<rot
->ngrp
; g
++)
850 if ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
851 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) )
853 fprintf(fp
, "# Rotation group %d (%s), slab distance %f nm\n", g
, erotg_names
[rotg
->eType
], rotg
->slab_dist
);
854 fprintf(fp
, "# The scalar tau is the torque [kJ/mol] in the direction of the rotation vector.\n");
855 fprintf(fp
, "# To obtain the vectorial torque, multiply tau with\n");
856 fprintf(fp
, "# rot_vec%d %10.3e %10.3e %10.3e\n", g
, rotg
->vec
[XX
], rotg
->vec
[YY
], rotg
->vec
[ZZ
]);
860 fprintf(fp
, "# The following columns will have the syntax (tau=torque for that slab):\n");
862 print_aligned_short(fp
, "t");
863 print_aligned_short(fp
, "grp");
864 print_aligned_short(fp
, "slab");
865 print_aligned(fp
, "tau");
866 print_aligned_short(fp
, "slab");
867 print_aligned(fp
, "tau");
868 fprintf(fp
, " ...\n");
875 static void swap_val(double* vec
, int i
, int j
)
885 static void swap_col(double **mat
, int i
, int j
)
887 double tmp
[3] = {mat
[0][j
], mat
[1][j
], mat
[2][j
]};
900 /* Eigenvectors are stored in columns of eigen_vec */
901 static void diagonalize_symmetric(
909 jacobi(matrix
,3,eigenval
,eigen_vec
,&n_rot
);
911 /* sort in ascending order */
912 if (eigenval
[0] > eigenval
[1])
914 swap_val(eigenval
, 0, 1);
915 swap_col(eigen_vec
, 0, 1);
917 if (eigenval
[1] > eigenval
[2])
919 swap_val(eigenval
, 1, 2);
920 swap_col(eigen_vec
, 1, 2);
922 if (eigenval
[0] > eigenval
[1])
924 swap_val(eigenval
, 0, 1);
925 swap_col(eigen_vec
, 0, 1);
930 static void align_with_z(
931 rvec
* s
, /* Structure to align */
936 rvec zet
= {0.0, 0.0, 1.0};
937 rvec rot_axis
={0.0, 0.0, 0.0};
938 rvec
*rotated_str
=NULL
;
944 snew(rotated_str
, natoms
);
946 /* Normalize the axis */
947 ooanorm
= 1.0/norm(axis
);
948 svmul(ooanorm
, axis
, axis
);
950 /* Calculate the angle for the fitting procedure */
951 cprod(axis
, zet
, rot_axis
);
952 angle
= acos(axis
[2]);
956 /* Calculate the rotation matrix */
957 calc_rotmat(rot_axis
, angle
*180.0/M_PI
, rotmat
);
959 /* Apply the rotation matrix to s */
960 for (i
=0; i
<natoms
; i
++)
966 rotated_str
[i
][j
] += rotmat
[j
][k
]*s
[i
][k
];
971 /* Rewrite the rotated structure to s */
972 for(i
=0; i
<natoms
; i
++)
976 s
[i
][j
]=rotated_str
[i
][j
];
984 static void calc_correl_matrix(rvec
* Xstr
, rvec
* Ystr
, double** Rmat
, int natoms
)
995 for (k
=0; k
<natoms
; k
++)
996 Rmat
[i
][j
] += Ystr
[k
][i
] * Xstr
[k
][j
];
1000 static void weigh_coords(rvec
* str
, real
* weight
, int natoms
)
1005 for(i
=0; i
<natoms
; i
++)
1008 str
[i
][j
] *= sqrt(weight
[i
]);
1013 static double opt_angle_analytic(
1026 double **Rmat
, **RtR
, **eigvec
;
1028 double V
[3][3], WS
[3][3];
1029 double rot_matrix
[3][3];
1033 /* Do not change the original coordinates */
1034 snew(ref_s_1
, natoms
);
1035 snew(act_s_1
, natoms
);
1036 for(i
=0; i
<natoms
; i
++)
1038 copy_rvec(ref_s
[i
], ref_s_1
[i
]);
1039 copy_rvec(act_s
[i
], act_s_1
[i
]);
1042 /* Translate the structures to the origin */
1043 shift
[XX
] = -ref_com
[XX
];
1044 shift
[YY
] = -ref_com
[YY
];
1045 shift
[ZZ
] = -ref_com
[ZZ
];
1046 translate_x(ref_s_1
, natoms
, shift
);
1048 shift
[XX
] = -act_com
[XX
];
1049 shift
[YY
] = -act_com
[YY
];
1050 shift
[ZZ
] = -act_com
[ZZ
];
1051 translate_x(act_s_1
, natoms
, shift
);
1053 /* Align rotation axis with z */
1054 align_with_z(ref_s_1
, natoms
, axis
);
1055 align_with_z(act_s_1
, natoms
, axis
);
1057 /* Correlation matrix */
1058 Rmat
= allocate_square_matrix(3);
1060 for (i
=0; i
<natoms
; i
++)
1066 /* Weight positions with sqrt(weight) */
1069 weigh_coords(ref_s_1
, weight
, natoms
);
1070 weigh_coords(act_s_1
, weight
, natoms
);
1073 /* Calculate correlation matrices R=YXt (X=ref_s; Y=act_s) */
1074 calc_correl_matrix(ref_s_1
, act_s_1
, Rmat
, natoms
);
1077 RtR
= allocate_square_matrix(3);
1084 RtR
[i
][j
] += Rmat
[k
][i
] * Rmat
[k
][j
];
1088 /* Diagonalize RtR */
1093 diagonalize_symmetric(RtR
, eigvec
, eigval
);
1094 swap_col(eigvec
,0,1);
1095 swap_col(eigvec
,1,2);
1096 swap_val(eigval
,0,1);
1097 swap_val(eigval
,1,2);
1111 WS
[i
][j
] = eigvec
[i
][j
] / sqrt(eigval
[j
]);
1119 V
[i
][j
] += Rmat
[i
][k
]*WS
[k
][j
];
1123 free_square_matrix(Rmat
, 3);
1125 /* Calculate optimal rotation matrix */
1128 rot_matrix
[i
][j
] = 0.0;
1135 rot_matrix
[i
][j
] += eigvec
[i
][k
]*V
[j
][k
];
1139 rot_matrix
[2][2] = 1.0;
1141 /* Determine the optimal rotation angle: */
1142 opt_angle
= (-1.0)*acos(rot_matrix
[0][0])*180.0/M_PI
;
1143 if (rot_matrix
[0][1] < 0.0)
1144 opt_angle
= (-1.0)*opt_angle
;
1146 /* Give back some memory */
1147 free_square_matrix(RtR
, 3);
1158 /* Determine actual angle of this slab by RMSD fit to the reference */
1159 /* Not parallelized, call this routine only on the master */
1160 static void flex_fit_angle(
1167 int i
,l
,n
,islab
,ind
;
1169 rvec
*fitcoords
=NULL
;
1170 rvec act_center
; /* Center of actual positions that are passed to the fit routine */
1171 rvec ref_center
; /* Same for the reference positions */
1172 double fitangle
; /* This will be the actual angle of the rotation group derived
1173 * from an RMSD fit to the reference structure at t=0 */
1177 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1178 real OOm_av
; /* 1/average_mass of a rotation group atom */
1179 real m_rel
; /* Relative mass of a rotation group atom */
1182 erg
=rotg
->enfrotgrp
;
1184 /* Average mass of a rotation group atom: */
1185 OOm_av
= erg
->invmass
*rotg
->nat
;
1187 /**********************************/
1188 /* First collect the data we need */
1189 /**********************************/
1191 /* Collect the data for the individual slabs */
1192 for (n
= erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1194 islab
= n
- erg
->slab_first
; /* slab index */
1195 sd
= &(rotg
->enfrotgrp
->slab_data
[islab
]);
1196 sd
->nat
= erg
->lastatom
[islab
]-erg
->firstatom
[islab
]+1;
1199 /* Loop over the relevant atoms in the slab */
1200 for (l
=erg
->firstatom
[islab
]; l
<=erg
->lastatom
[islab
]; l
++)
1202 /* Current position of this atom: x[ii][XX/YY/ZZ] */
1203 copy_rvec(erg
->xc
[l
], curr_x
);
1205 /* The (unrotated) reference position of this atom is copied to ref_x.
1206 * Beware, the xc coords have been sorted in do_flexible */
1207 copy_rvec(erg
->xc_ref_sorted
[l
], ref_x
);
1209 /* Save data for doing angular RMSD fit later */
1210 /* Save the current atom position */
1211 copy_rvec(curr_x
, sd
->x
[ind
]);
1212 /* Save the corresponding reference position */
1213 copy_rvec(ref_x
, sd
->ref
[ind
]);
1215 /* Maybe also mass-weighting was requested. If yes, additionally
1216 * multiply the weights with the relative mass of the atom. If not,
1217 * multiply with unity. */
1218 m_rel
= erg
->mc_sorted
[l
]*OOm_av
;
1220 /* Save the weight for this atom in this slab */
1221 sd
->weight
[ind
] = gaussian_weight(curr_x
, rotg
, n
) * m_rel
;
1223 /* Next atom in this slab */
1228 /* Get the center of the whole rotation group. Note, again, the erg->xc have
1229 * been sorted in do_flexible */
1230 get_center(erg
->xc
, erg
->mc_sorted
, rotg
->nat
, act_center
);
1232 /******************************/
1233 /* Now do the fit calculation */
1234 /******************************/
1236 /* === Determine the optimal fit angle for the whole rotation group === */
1237 if (rotg
->eFittype
== erotgFitNORM
)
1239 /* Normalize every position to it's reference length
1240 * prior to performing the fit */
1241 for (i
=0; i
<rotg
->nat
; i
++)
1243 /* First put the center of the positions into the origin */
1244 rvec_sub(erg
->xc
[i
], act_center
, coord
);
1245 /* Determine the scaling factor for the length: */
1246 scal
= erg
->xc_ref_length
[erg
->xc_sortind
[i
]] / norm(coord
);
1247 /* Get position, multiply with the scaling factor and save in buf[i] */
1248 svmul(scal
, coord
, erg
->xc_norm
[i
]);
1250 fitcoords
= erg
->xc_norm
;
1254 fitcoords
= erg
->xc
;
1256 /* Note that from the point of view of the current positions, the reference has rotated backwards,
1257 * but we want to output the angle relative to the fixed reference, therefore the minus sign. */
1258 fitangle
= -opt_angle_analytic(erg
->xc_ref_sorted
, fitcoords
, erg
->mc_sorted
,
1259 rotg
->nat
, erg
->xc_ref_center
, act_center
, rotg
->vec
);
1260 fprintf(fp
, "%12.3e%6d%12.3f%12.3lf", t
, g
, degangle
, fitangle
);
1263 /* === Now do RMSD fitting for each slab === */
1264 /* We require at least SLAB_MIN_ATOMS in a slab, such that the fit makes sense. */
1265 #define SLAB_MIN_ATOMS 4
1267 for (n
= erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1269 islab
= n
- erg
->slab_first
; /* slab index */
1270 sd
= &(rotg
->enfrotgrp
->slab_data
[islab
]);
1271 if (sd
->nat
>= SLAB_MIN_ATOMS
)
1273 /* Get the center of the slabs reference and current positions */
1274 get_center(sd
->ref
, sd
->weight
, sd
->nat
, ref_center
);
1275 get_center(sd
->x
, sd
->weight
, sd
->nat
, act_center
);
1276 if (rotg
->eFittype
== erotgFitNORM
)
1278 /* Normalize every position to it's reference length
1279 * prior to performing the fit */
1280 for (i
=0; i
<sd
->nat
;i
++) /* Center */
1282 rvec_dec(sd
->ref
[i
], ref_center
);
1283 rvec_dec(sd
->x
[i
] , act_center
);
1284 /* Normalize x_i such that it gets the same length as ref_i */
1285 svmul( norm(sd
->ref
[i
])/norm(sd
->x
[i
]), sd
->x
[i
], sd
->x
[i
] );
1287 /* We already subtracted the centers */
1288 clear_rvec(ref_center
);
1289 clear_rvec(act_center
);
1291 fitangle
= -opt_angle_analytic(sd
->ref
, sd
->x
, sd
->weight
, sd
->nat
, ref_center
, act_center
, rotg
->vec
);
1292 fprintf(fp
, "%6d%6d%12.3f", n
, sd
->nat
, fitangle
);
1297 #undef SLAB_MIN_ATOMS
1301 /* Shift x with is */
1302 static inline void shift_single_coord(matrix box
, rvec x
, const ivec is
)
1313 x
[XX
] += tx
*box
[XX
][XX
]+ty
*box
[YY
][XX
]+tz
*box
[ZZ
][XX
];
1314 x
[YY
] += ty
*box
[YY
][YY
]+tz
*box
[ZZ
][YY
];
1315 x
[ZZ
] += tz
*box
[ZZ
][ZZ
];
1318 x
[XX
] += tx
*box
[XX
][XX
];
1319 x
[YY
] += ty
*box
[YY
][YY
];
1320 x
[ZZ
] += tz
*box
[ZZ
][ZZ
];
1325 /* Determine the 'home' slab of this atom which is the
1326 * slab with the highest Gaussian weight of all */
1327 #define round(a) (int)(a+0.5)
1328 static inline int get_homeslab(
1329 rvec curr_x
, /* The position for which the home slab shall be determined */
1330 rvec rotvec
, /* The rotation vector */
1331 real slabdist
) /* The slab distance */
1336 /* The distance of the atom to the coordinate center (where the
1337 * slab with index 0) is */
1338 dist
= iprod(rotvec
, curr_x
);
1340 return round(dist
/ slabdist
);
1344 /* For a local atom determine the relevant slabs, i.e. slabs in
1345 * which the gaussian is larger than min_gaussian
1347 static int get_single_atom_gaussians(
1355 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1358 erg
=rotg
->enfrotgrp
;
1360 /* Determine the 'home' slab of this atom: */
1361 homeslab
= get_homeslab(curr_x
, rotg
->vec
, rotg
->slab_dist
);
1363 /* First determine the weight in the atoms home slab: */
1364 g
= gaussian_weight(curr_x
, rotg
, homeslab
);
1366 erg
->gn_atom
[count
] = g
;
1367 erg
->gn_slabind
[count
] = homeslab
;
1371 /* Determine the max slab */
1373 while (g
> rotg
->min_gaussian
)
1376 g
= gaussian_weight(curr_x
, rotg
, slab
);
1377 erg
->gn_slabind
[count
]=slab
;
1378 erg
->gn_atom
[count
]=g
;
1383 /* Determine the max slab */
1388 g
= gaussian_weight(curr_x
, rotg
, slab
);
1389 erg
->gn_slabind
[count
]=slab
;
1390 erg
->gn_atom
[count
]=g
;
1393 while (g
> rotg
->min_gaussian
);
1400 static void flex2_precalc_inner_sum(t_rotgrp
*rotg
, t_commrec
*cr
)
1403 rvec xi
; /* positions in the i-sum */
1404 rvec xcn
, ycn
; /* the current and the reference slab centers */
1407 rvec rin
; /* Helper variables */
1410 real OOpsii
,OOpsiistar
;
1411 real sin_rin
; /* s_ii.r_ii */
1412 rvec s_in
,tmpvec
,tmpvec2
;
1413 real mi
,wi
; /* Mass-weighting of the positions */
1415 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1418 erg
=rotg
->enfrotgrp
;
1419 N_M
= rotg
->nat
* erg
->invmass
;
1421 /* Loop over all slabs that contain something */
1422 for (n
=erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1424 islab
= n
- erg
->slab_first
; /* slab index */
1426 /* The current center of this slab is saved in xcn: */
1427 copy_rvec(erg
->slab_center
[islab
], xcn
);
1428 /* ... and the reference center in ycn: */
1429 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
1431 /*** D. Calculate the whole inner sum used for second and third sum */
1432 /* For slab n, we need to loop over all atoms i again. Since we sorted
1433 * the atoms with respect to the rotation vector, we know that it is sufficient
1434 * to calculate from firstatom to lastatom only. All other contributions will
1436 clear_rvec(innersumvec
);
1437 for (i
= erg
->firstatom
[islab
]; i
<= erg
->lastatom
[islab
]; i
++)
1439 /* Coordinate xi of this atom */
1440 copy_rvec(erg
->xc
[i
],xi
);
1443 gaussian_xi
= gaussian_weight(xi
,rotg
,n
);
1444 mi
= erg
->mc_sorted
[i
]; /* need the sorted mass here */
1448 copy_rvec(erg
->xc_ref_sorted
[i
],yi0
); /* Reference position yi0 */
1449 rvec_sub(yi0
, ycn
, tmpvec2
); /* tmpvec2 = yi0 - ycn */
1450 mvmul(erg
->rotmat
, tmpvec2
, rin
); /* rin = Omega.(yi0 - ycn) */
1452 /* Calculate psi_i* and sin */
1453 rvec_sub(xi
, xcn
, tmpvec2
); /* tmpvec2 = xi - xcn */
1454 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x (xi - xcn) */
1455 OOpsiistar
= norm2(tmpvec
)+rotg
->eps
; /* OOpsii* = 1/psii* = |v x (xi-xcn)|^2 + eps */
1456 OOpsii
= norm(tmpvec
); /* OOpsii = 1 / psii = |v x (xi - xcn)| */
1458 /* v x (xi - xcn) */
1459 unitv(tmpvec
, s_in
); /* sin = ---------------- */
1460 /* |v x (xi - xcn)| */
1462 sin_rin
=iprod(s_in
,rin
); /* sin_rin = sin . rin */
1464 /* Now the whole sum */
1465 fac
= OOpsii
/OOpsiistar
;
1466 svmul(fac
, rin
, tmpvec
);
1467 fac2
= fac
*fac
*OOpsii
;
1468 svmul(fac2
*sin_rin
, s_in
, tmpvec2
);
1469 rvec_dec(tmpvec
, tmpvec2
);
1471 svmul(wi
*gaussian_xi
*sin_rin
, tmpvec
, tmpvec2
);
1473 rvec_inc(innersumvec
,tmpvec2
);
1474 } /* now we have the inner sum, used both for sum2 and sum3 */
1476 /* Save it to be used in do_flex2_lowlevel */
1477 copy_rvec(innersumvec
, erg
->slab_innersumvec
[islab
]);
1478 } /* END of loop over slabs */
1482 static void flex_precalc_inner_sum(t_rotgrp
*rotg
, t_commrec
*cr
)
1485 rvec xi
; /* position */
1486 rvec xcn
, ycn
; /* the current and the reference slab centers */
1487 rvec qin
,rin
; /* q_i^n and r_i^n */
1490 rvec innersumvec
; /* Inner part of sum_n2 */
1491 real gaussian_xi
; /* Gaussian weight gn(xi) */
1492 real mi
,wi
; /* Mass-weighting of the positions */
1495 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1498 erg
=rotg
->enfrotgrp
;
1499 N_M
= rotg
->nat
* erg
->invmass
;
1501 /* Loop over all slabs that contain something */
1502 for (n
=erg
->slab_first
; n
<= erg
->slab_last
; n
++)
1504 islab
= n
- erg
->slab_first
; /* slab index */
1506 /* The current center of this slab is saved in xcn: */
1507 copy_rvec(erg
->slab_center
[islab
], xcn
);
1508 /* ... and the reference center in ycn: */
1509 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
1511 /* For slab n, we need to loop over all atoms i again. Since we sorted
1512 * the atoms with respect to the rotation vector, we know that it is sufficient
1513 * to calculate from firstatom to lastatom only. All other contributions will
1515 clear_rvec(innersumvec
);
1516 for (i
=erg
->firstatom
[islab
]; i
<=erg
->lastatom
[islab
]; i
++)
1518 /* Coordinate xi of this atom */
1519 copy_rvec(erg
->xc
[i
],xi
);
1522 gaussian_xi
= gaussian_weight(xi
,rotg
,n
);
1523 mi
= erg
->mc_sorted
[i
]; /* need the sorted mass here */
1526 /* Calculate rin and qin */
1527 rvec_sub(erg
->xc_ref_sorted
[i
], ycn
, tmpvec
); /* tmpvec = yi0-ycn */
1528 mvmul(erg
->rotmat
, tmpvec
, rin
); /* rin = Omega.(yi0 - ycn) */
1529 cprod(rotg
->vec
, rin
, tmpvec
); /* tmpvec = v x Omega*(yi0-ycn) */
1531 /* v x Omega*(yi0-ycn) */
1532 unitv(tmpvec
, qin
); /* qin = --------------------- */
1533 /* |v x Omega*(yi0-ycn)| */
1536 rvec_sub(xi
, xcn
, tmpvec
); /* tmpvec = xi-xcn */
1537 bin
= iprod(qin
, tmpvec
); /* bin = qin*(xi-xcn) */
1539 svmul(wi
*gaussian_xi
*bin
, qin
, tmpvec
);
1541 /* Add this contribution to the inner sum: */
1542 rvec_add(innersumvec
, tmpvec
, innersumvec
);
1543 } /* now we have the inner sum vector S^n for this slab */
1544 /* Save it to be used in do_flex_lowlevel */
1545 copy_rvec(innersumvec
, erg
->slab_innersumvec
[islab
]);
1550 static real
do_flex2_lowlevel(
1552 real sigma
, /* The Gaussian width sigma */
1554 gmx_bool bCalcTorque
,
1558 int count
,ic
,ii
,j
,m
,n
,islab
,iigrp
;
1559 rvec xj
; /* position in the i-sum */
1560 rvec yj0
; /* the reference position in the j-sum */
1561 rvec xcn
, ycn
; /* the current and the reference slab centers */
1562 real V
; /* This node's part of the rotation pot. energy */
1563 real gaussian_xj
; /* Gaussian weight */
1567 rvec rjn
; /* Helper variables */
1570 real OOpsij
,OOpsijstar
;
1571 real OOsigma2
; /* 1/(sigma^2) */
1574 rvec sjn
,tmpvec
,tmpvec2
;
1575 rvec sum1vec_part
,sum1vec
,sum2vec_part
,sum2vec
,sum3vec
,sum4vec
,innersumvec
;
1577 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1578 real mj
,wj
; /* Mass-weighting of the positions */
1580 real Wjn
; /* g_n(x_j) m_j / Mjn */
1582 /* To calculate the torque per slab */
1583 rvec slab_force
; /* Single force from slab n on one atom */
1584 rvec slab_sum1vec_part
;
1585 real slab_sum3part
,slab_sum4part
;
1586 rvec slab_sum1vec
, slab_sum2vec
, slab_sum3vec
, slab_sum4vec
;
1589 erg
=rotg
->enfrotgrp
;
1591 /* Pre-calculate the inner sums, so that we do not have to calculate
1592 * them again for every atom */
1593 flex2_precalc_inner_sum(rotg
, cr
);
1595 /********************************************************/
1596 /* Main loop over all local atoms of the rotation group */
1597 /********************************************************/
1598 N_M
= rotg
->nat
* erg
->invmass
;
1600 OOsigma2
= 1.0 / (sigma
*sigma
);
1601 for (j
=0; j
<erg
->nat_loc
; j
++)
1603 /* Local index of a rotation group atom */
1604 ii
= erg
->ind_loc
[j
];
1605 /* Position of this atom in the collective array */
1606 iigrp
= erg
->xc_ref_ind
[j
];
1607 /* Mass-weighting */
1608 mj
= erg
->mc
[iigrp
]; /* need the unsorted mass here */
1611 /* Current position of this atom: x[ii][XX/YY/ZZ]
1612 * Note that erg->xc_center contains the center of mass in case the flex2-t
1613 * potential was chosen. For the flex2 potential erg->xc_center must be
1615 rvec_sub(x
[ii
], erg
->xc_center
, xj
);
1617 /* Shift this atom such that it is near its reference */
1618 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
1620 /* Determine the slabs to loop over, i.e. the ones with contributions
1621 * larger than min_gaussian */
1622 count
= get_single_atom_gaussians(xj
, cr
, rotg
);
1624 clear_rvec(sum1vec_part
);
1625 clear_rvec(sum2vec_part
);
1628 /* Loop over the relevant slabs for this atom */
1629 for (ic
=0; ic
< count
; ic
++)
1631 n
= erg
->gn_slabind
[ic
];
1633 /* Get the precomputed Gaussian value of curr_slab for curr_x */
1634 gaussian_xj
= erg
->gn_atom
[ic
];
1636 islab
= n
- erg
->slab_first
; /* slab index */
1638 /* The (unrotated) reference position of this atom is copied to yj0: */
1639 copy_rvec(rotg
->x_ref
[iigrp
], yj0
);
1641 beta
= calc_beta(xj
, rotg
,n
);
1643 /* The current center of this slab is saved in xcn: */
1644 copy_rvec(erg
->slab_center
[islab
], xcn
);
1645 /* ... and the reference center in ycn: */
1646 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
1648 rvec_sub(yj0
, ycn
, tmpvec2
); /* tmpvec2 = yj0 - ycn */
1651 mvmul(erg
->rotmat
, tmpvec2
, rjn
); /* rjn = Omega.(yj0 - ycn) */
1653 /* Subtract the slab center from xj */
1654 rvec_sub(xj
, xcn
, tmpvec2
); /* tmpvec2 = xj - xcn */
1657 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x (xj - xcn) */
1659 OOpsijstar
= norm2(tmpvec
)+rotg
->eps
; /* OOpsij* = 1/psij* = |v x (xj-xcn)|^2 + eps */
1661 numerator
= sqr(iprod(tmpvec
, rjn
));
1663 /*********************************/
1664 /* Add to the rotation potential */
1665 /*********************************/
1666 V
+= 0.5*rotg
->k
*wj
*gaussian_xj
*numerator
/OOpsijstar
;
1669 /*************************************/
1670 /* Now calculate the force on atom j */
1671 /*************************************/
1673 OOpsij
= norm(tmpvec
); /* OOpsij = 1 / psij = |v x (xj - xcn)| */
1675 /* v x (xj - xcn) */
1676 unitv(tmpvec
, sjn
); /* sjn = ---------------- */
1677 /* |v x (xj - xcn)| */
1679 sjn_rjn
=iprod(sjn
,rjn
); /* sjn_rjn = sjn . rjn */
1682 /*** A. Calculate the first of the four sum terms: ****************/
1683 fac
= OOpsij
/OOpsijstar
;
1684 svmul(fac
, rjn
, tmpvec
);
1685 fac2
= fac
*fac
*OOpsij
;
1686 svmul(fac2
*sjn_rjn
, sjn
, tmpvec2
);
1687 rvec_dec(tmpvec
, tmpvec2
);
1688 fac2
= wj
*gaussian_xj
; /* also needed for sum4 */
1689 svmul(fac2
*sjn_rjn
, tmpvec
, slab_sum1vec_part
);
1690 /********************/
1691 /*** Add to sum1: ***/
1692 /********************/
1693 rvec_inc(sum1vec_part
, slab_sum1vec_part
); /* sum1 still needs to vector multiplied with v */
1695 /*** B. Calculate the forth of the four sum terms: ****************/
1696 betasigpsi
= beta
*OOsigma2
*OOpsij
; /* this is also needed for sum3 */
1697 /********************/
1698 /*** Add to sum4: ***/
1699 /********************/
1700 slab_sum4part
= fac2
*betasigpsi
*fac
*sjn_rjn
*sjn_rjn
; /* Note that fac is still valid from above */
1701 sum4
+= slab_sum4part
;
1703 /*** C. Calculate Wjn for second and third sum */
1704 /* Note that we can safely divide by slab_weights since we check in
1705 * get_slab_centers that it is non-zero. */
1706 Wjn
= gaussian_xj
*mj
/erg
->slab_weights
[islab
];
1708 /* We already have precalculated the inner sum for slab n */
1709 copy_rvec(erg
->slab_innersumvec
[islab
], innersumvec
);
1711 /* Weigh the inner sum vector with Wjn */
1712 svmul(Wjn
, innersumvec
, innersumvec
);
1714 /*** E. Calculate the second of the four sum terms: */
1715 /********************/
1716 /*** Add to sum2: ***/
1717 /********************/
1718 rvec_inc(sum2vec_part
, innersumvec
); /* sum2 still needs to be vector crossproduct'ed with v */
1720 /*** F. Calculate the third of the four sum terms: */
1721 slab_sum3part
= betasigpsi
* iprod(sjn
, innersumvec
);
1722 sum3
+= slab_sum3part
; /* still needs to be multiplied with v */
1724 /*** G. Calculate the torque on the local slab's axis: */
1728 cprod(slab_sum1vec_part
, rotg
->vec
, slab_sum1vec
);
1730 cprod(innersumvec
, rotg
->vec
, slab_sum2vec
);
1732 svmul(slab_sum3part
, rotg
->vec
, slab_sum3vec
);
1734 svmul(slab_sum4part
, rotg
->vec
, slab_sum4vec
);
1736 /* The force on atom ii from slab n only: */
1737 for (m
=0; m
<DIM
; m
++)
1738 slab_force
[m
] = rotg
->k
* (-slab_sum1vec
[m
] + slab_sum2vec
[m
] - slab_sum3vec
[m
] + 0.5*slab_sum4vec
[m
]);
1740 erg
->slab_torque_v
[islab
] += torque(rotg
->vec
, slab_force
, xj
, xcn
);
1742 } /* END of loop over slabs */
1744 /* Construct the four individual parts of the vector sum: */
1745 cprod(sum1vec_part
, rotg
->vec
, sum1vec
); /* sum1vec = { } x v */
1746 cprod(sum2vec_part
, rotg
->vec
, sum2vec
); /* sum2vec = { } x v */
1747 svmul(sum3
, rotg
->vec
, sum3vec
); /* sum3vec = { } . v */
1748 svmul(sum4
, rotg
->vec
, sum4vec
); /* sum4vec = { } . v */
1750 /* Store the additional force so that it can be added to the force
1751 * array after the normal forces have been evaluated */
1752 for (m
=0; m
<DIM
; m
++)
1753 erg
->f_rot_loc
[j
][m
] = rotg
->k
* (-sum1vec
[m
] + sum2vec
[m
] - sum3vec
[m
] + 0.5*sum4vec
[m
]);
1756 fprintf(stderr
," FORCE on ATOM %d/%d = (%15.8f %15.8f %15.8f) \n",
1757 j
,ii
,erg
->f_rot_loc
[j
][XX
], erg
->f_rot_loc
[j
][YY
], erg
->f_rot_loc
[j
][ZZ
]);
1761 fprintf(stderr
, "sum1: %15.8f %15.8f %15.8f\n", -rotg
->k
*sum1vec
[XX
], -rotg
->k
*sum1vec
[YY
], -rotg
->k
*sum1vec
[ZZ
]);
1762 fprintf(stderr
, "sum2: %15.8f %15.8f %15.8f\n", rotg
->k
*sum2vec
[XX
], rotg
->k
*sum2vec
[YY
], rotg
->k
*sum2vec
[ZZ
]);
1763 fprintf(stderr
, "sum3: %15.8f %15.8f %15.8f\n", -rotg
->k
*sum3vec
[XX
], -rotg
->k
*sum3vec
[YY
], -rotg
->k
*sum3vec
[ZZ
]);
1764 fprintf(stderr
, "sum4: %15.8f %15.8f %15.8f\n", 0.5*rotg
->k
*sum4vec
[XX
], 0.5*rotg
->k
*sum4vec
[YY
], 0.5*rotg
->k
*sum4vec
[ZZ
]);
1766 } /* END of loop over local atoms */
1769 fprintf(stderr
, "THE POTENTIAL IS V=%f\n", V
);
1776 static real
do_flex_lowlevel(
1778 real sigma
, /* The Gaussian width sigma */
1780 gmx_bool bCalcTorque
,
1784 int count
,ic
,ii
,j
,m
,n
,islab
,iigrp
;
1785 rvec xj
,yj0
; /* current and reference position */
1786 rvec xcn
, ycn
; /* the current and the reference slab centers */
1787 rvec xj_xcn
; /* xj - xcn */
1788 rvec qjn
; /* q_i^n */
1789 rvec sum_n1
,sum_n2
; /* Two contributions to the rotation force */
1790 rvec innersumvec
; /* Inner part of sum_n2 */
1792 rvec force_n
; /* Single force from slab n on one atom */
1793 rvec tmpvec
,tmpvec2
,tmp_f
; /* Helper variables */
1794 real V
; /* The rotation potential energy */
1795 real OOsigma2
; /* 1/(sigma^2) */
1796 real beta
; /* beta_n(xj) */
1797 real bjn
; /* b_j^n */
1798 real gaussian_xj
; /* Gaussian weight gn(xj) */
1799 real betan_xj_sigma2
;
1800 real mj
,wj
; /* Mass-weighting of the positions */
1802 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
1805 erg
=rotg
->enfrotgrp
;
1807 /* Pre-calculate the inner sums, so that we do not have to calculate
1808 * them again for every atom */
1809 flex_precalc_inner_sum(rotg
, cr
);
1811 /********************************************************/
1812 /* Main loop over all local atoms of the rotation group */
1813 /********************************************************/
1814 OOsigma2
= 1.0/(sigma
*sigma
);
1815 N_M
= rotg
->nat
* erg
->invmass
;
1817 for (j
=0; j
<erg
->nat_loc
; j
++)
1819 /* Local index of a rotation group atom */
1820 ii
= erg
->ind_loc
[j
];
1821 /* Position of this atom in the collective array */
1822 iigrp
= erg
->xc_ref_ind
[j
];
1823 /* Mass-weighting */
1824 mj
= erg
->mc
[iigrp
]; /* need the unsorted mass here */
1827 /* Current position of this atom: x[ii][XX/YY/ZZ]
1828 * Note that erg->xc_center contains the center of mass in case the flex-t
1829 * potential was chosen. For the flex potential erg->xc_center must be
1831 rvec_sub(x
[ii
], erg
->xc_center
, xj
);
1833 /* Shift this atom such that it is near its reference */
1834 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
1836 /* Determine the slabs to loop over, i.e. the ones with contributions
1837 * larger than min_gaussian */
1838 count
= get_single_atom_gaussians(xj
, cr
, rotg
);
1843 /* Loop over the relevant slabs for this atom */
1844 for (ic
=0; ic
< count
; ic
++)
1846 n
= erg
->gn_slabind
[ic
];
1848 /* Get the precomputed Gaussian for xj in slab n */
1849 gaussian_xj
= erg
->gn_atom
[ic
];
1851 islab
= n
- erg
->slab_first
; /* slab index */
1853 /* The (unrotated) reference position of this atom is saved in yj0: */
1854 copy_rvec(rotg
->x_ref
[iigrp
], yj0
);
1856 beta
= calc_beta(xj
, rotg
, n
);
1858 /* The current center of this slab is saved in xcn: */
1859 copy_rvec(erg
->slab_center
[islab
], xcn
);
1860 /* ... and the reference center in ycn: */
1861 copy_rvec(erg
->slab_center_ref
[islab
+erg
->slab_buffer
], ycn
);
1863 rvec_sub(yj0
, ycn
, tmpvec
); /* tmpvec = yj0 - ycn */
1866 mvmul(erg
->rotmat
, tmpvec
, tmpvec2
); /* tmpvec2 = Omega.(yj0-ycn) */
1868 /* Subtract the slab center from xj */
1869 rvec_sub(xj
, xcn
, xj_xcn
); /* xj_xcn = xj - xcn */
1872 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x Omega.(xj-xcn) */
1874 /* v x Omega.(xj-xcn) */
1875 unitv(tmpvec
,qjn
); /* qjn = -------------------- */
1876 /* |v x Omega.(xj-xcn)| */
1878 bjn
= iprod(qjn
, xj_xcn
); /* bjn = qjn * (xj - xcn) */
1880 /*********************************/
1881 /* Add to the rotation potential */
1882 /*********************************/
1883 V
+= 0.5*rotg
->k
*wj
*gaussian_xj
*sqr(bjn
);
1885 /****************************************************************/
1886 /* sum_n1 will typically be the main contribution to the force: */
1887 /****************************************************************/
1888 betan_xj_sigma2
= beta
*OOsigma2
; /* beta_n(xj)/sigma^2 */
1890 /* The next lines calculate
1891 * qjn - (bjn*beta(xj)/(2sigma^2))v */
1892 svmul(bjn
*0.5*betan_xj_sigma2
, rotg
->vec
, tmpvec2
);
1893 rvec_sub(qjn
,tmpvec2
,tmpvec
);
1895 /* Multiply with gn(xj)*bjn: */
1896 svmul(gaussian_xj
*bjn
,tmpvec
,tmpvec2
);
1899 rvec_inc(sum_n1
,tmpvec2
);
1901 /* We already have precalculated the Sn term for slab n */
1902 copy_rvec(erg
->slab_innersumvec
[islab
], s_n
);
1904 svmul(betan_xj_sigma2
*iprod(s_n
, xj_xcn
), rotg
->vec
, tmpvec
); /* tmpvec = ---------- s_n (xj-xcn) */
1907 rvec_sub(s_n
, tmpvec
, innersumvec
);
1909 /* We can safely divide by slab_weights since we check in get_slab_centers
1910 * that it is non-zero. */
1911 svmul(gaussian_xj
/erg
->slab_weights
[islab
], innersumvec
, innersumvec
);
1913 rvec_add(sum_n2
, innersumvec
, sum_n2
);
1915 GMX_MPE_LOG(ev_inner_loop_finish
);
1917 /* Calculate the torque: */
1920 /* The force on atom ii from slab n only: */
1921 rvec_sub(innersumvec
, tmpvec2
, force_n
);
1922 svmul(rotg
->k
, force_n
, force_n
);
1923 erg
->slab_torque_v
[islab
] += torque(rotg
->vec
, force_n
, xj
, xcn
);
1925 } /* END of loop over slabs */
1927 /* Put both contributions together: */
1928 svmul(wj
, sum_n1
, sum_n1
);
1929 svmul(mj
, sum_n2
, sum_n2
);
1930 rvec_sub(sum_n2
,sum_n1
,tmp_f
); /* F = -grad V */
1932 /* Store the additional force so that it can be added to the force
1933 * array after the normal forces have been evaluated */
1934 for(m
=0; m
<DIM
; m
++)
1935 erg
->f_rot_loc
[j
][m
] = rotg
->k
*tmp_f
[m
];
1937 fprintf(stderr
," FORCE on atom %d = %15.8f %15.8f %15.8f 1: %15.8f %15.8f %15.8f 2: %15.8f %15.8f %15.8f\n", iigrp
,
1938 rotg
->k
*tmp_f
[XX
] , rotg
->k
*tmp_f
[YY
] , rotg
->k
*tmp_f
[ZZ
] ,
1939 -rotg
->k
*sum_n1
[XX
], -rotg
->k
*sum_n1
[YY
], -rotg
->k
*sum_n1
[ZZ
],
1940 rotg
->k
*sum_n2
[XX
], rotg
->k
*sum_n2
[YY
], rotg
->k
*sum_n2
[ZZ
]);
1942 } /* END of loop over local atoms */
1948 static void print_coordinates(t_commrec
*cr
, t_rotgrp
*rotg
, rvec x
[], matrix box
, int step
)
1952 static char buf
[STRLEN
];
1953 static gmx_bool bFirst
=1;
1958 sprintf(buf
, "coords%d.txt", cr
->nodeid
);
1959 fp
= fopen(buf
, "w");
1963 fprintf(fp
, "\nStep %d\n", step
);
1964 fprintf(fp
, "box: %f %f %f %f %f %f %f %f %f\n",
1965 box
[XX
][XX
], box
[XX
][YY
], box
[XX
][ZZ
],
1966 box
[YY
][XX
], box
[YY
][YY
], box
[YY
][ZZ
],
1967 box
[ZZ
][XX
], box
[ZZ
][ZZ
], box
[ZZ
][ZZ
]);
1968 for (i
=0; i
<rotg
->nat
; i
++)
1970 fprintf(fp
, "%4d %f %f %f\n", i
,
1971 erg
->xc
[i
][XX
], erg
->xc
[i
][YY
], erg
->xc
[i
][ZZ
]);
1979 static int projection_compare(const void *a
, const void *b
)
1981 sort_along_vec_t
*xca
, *xcb
;
1984 xca
= (sort_along_vec_t
*)a
;
1985 xcb
= (sort_along_vec_t
*)b
;
1987 if (xca
->xcproj
< xcb
->xcproj
)
1989 else if (xca
->xcproj
> xcb
->xcproj
)
1996 static void sort_collective_coordinates(
1997 t_rotgrp
*rotg
, /* Rotation group */
1998 sort_along_vec_t
*data
) /* Buffer for sorting the positions */
2001 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2004 erg
=rotg
->enfrotgrp
;
2006 /* The projection of the position vector on the rotation vector is
2007 * the relevant value for sorting. Fill the 'data' structure */
2008 for (i
=0; i
<rotg
->nat
; i
++)
2010 data
[i
].xcproj
= iprod(erg
->xc
[i
], rotg
->vec
); /* sort criterium */
2011 data
[i
].m
= erg
->mc
[i
];
2013 copy_rvec(erg
->xc
[i
] , data
[i
].x
);
2014 copy_rvec(rotg
->x_ref
[i
], data
[i
].x_ref
);
2016 /* Sort the 'data' structure */
2017 gmx_qsort(data
, rotg
->nat
, sizeof(sort_along_vec_t
), projection_compare
);
2019 /* Copy back the sorted values */
2020 for (i
=0; i
<rotg
->nat
; i
++)
2022 copy_rvec(data
[i
].x
, erg
->xc
[i
] );
2023 copy_rvec(data
[i
].x_ref
, erg
->xc_ref_sorted
[i
]);
2024 erg
->mc_sorted
[i
] = data
[i
].m
;
2025 erg
->xc_sortind
[i
] = data
[i
].ind
;
2030 /* For each slab, get the first and the last index of the sorted atom
2032 static void get_firstlast_atom_per_slab(t_rotgrp
*rotg
, t_commrec
*cr
)
2036 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2039 erg
=rotg
->enfrotgrp
;
2041 GMX_MPE_LOG(ev_get_firstlast_start
);
2043 /* Find the first atom that needs to enter the calculation for each slab */
2044 n
= erg
->slab_first
; /* slab */
2045 i
= 0; /* start with the first atom */
2048 /* Find the first atom that significantly contributes to this slab */
2049 do /* move forward in position until a large enough beta is found */
2051 beta
= calc_beta(erg
->xc
[i
], rotg
, n
);
2053 } while ((beta
< -erg
->max_beta
) && (i
< rotg
->nat
));
2055 islab
= n
- erg
->slab_first
; /* slab index */
2056 erg
->firstatom
[islab
] = i
;
2057 /* Proceed to the next slab */
2059 } while (n
<= erg
->slab_last
);
2061 /* Find the last atom for each slab */
2062 n
= erg
->slab_last
; /* start with last slab */
2063 i
= rotg
->nat
-1; /* start with the last atom */
2066 do /* move backward in position until a large enough beta is found */
2068 beta
= calc_beta(erg
->xc
[i
], rotg
, n
);
2070 } while ((beta
> erg
->max_beta
) && (i
> -1));
2072 islab
= n
- erg
->slab_first
; /* slab index */
2073 erg
->lastatom
[islab
] = i
;
2074 /* Proceed to the next slab */
2076 } while (n
>= erg
->slab_first
);
2078 GMX_MPE_LOG(ev_get_firstlast_finish
);
2082 /* Determine the very first and very last slab that needs to be considered
2083 * For the first slab that needs to be considered, we have to find the smallest
2086 * x_first * v - n*Delta_x <= beta_max
2088 * slab index n, slab distance Delta_x, rotation vector v. For the last slab we
2089 * have to find the largest n that obeys
2091 * x_last * v - n*Delta_x >= -beta_max
2094 static inline int get_first_slab(
2095 t_rotgrp
*rotg
, /* The rotation group (inputrec data) */
2096 real max_beta
, /* The max_beta value, instead of min_gaussian */
2097 rvec firstatom
) /* First atom after sorting along the rotation vector v */
2099 /* Find the first slab for the first atom */
2100 return ceil((iprod(firstatom
, rotg
->vec
) - max_beta
)/rotg
->slab_dist
);
2104 static inline int get_last_slab(
2105 t_rotgrp
*rotg
, /* The rotation group (inputrec data) */
2106 real max_beta
, /* The max_beta value, instead of min_gaussian */
2107 rvec lastatom
) /* Last atom along v */
2109 /* Find the last slab for the last atom */
2110 return floor((iprod(lastatom
, rotg
->vec
) + max_beta
)/rotg
->slab_dist
);
2114 static void get_firstlast_slab_check(
2115 t_rotgrp
*rotg
, /* The rotation group (inputrec data) */
2116 t_gmx_enfrotgrp
*erg
, /* The rotation group (data only accessible in this file) */
2117 rvec firstatom
, /* First atom after sorting along the rotation vector v */
2118 rvec lastatom
, /* Last atom along v */
2119 int g
, /* The rotation group number */
2122 erg
->slab_first
= get_first_slab(rotg
, erg
->max_beta
, firstatom
);
2123 erg
->slab_last
= get_last_slab(rotg
, erg
->max_beta
, lastatom
);
2125 /* Check whether we have reference data to compare against */
2126 if (erg
->slab_first
< erg
->slab_first_ref
)
2127 gmx_fatal(FARGS
, "%s No reference data for first slab (n=%d), unable to proceed.",
2128 RotStr
, erg
->slab_first
);
2130 /* Check whether we have reference data to compare against */
2131 if (erg
->slab_last
> erg
->slab_last_ref
)
2132 gmx_fatal(FARGS
, "%s No reference data for last slab (n=%d), unable to proceed.",
2133 RotStr
, erg
->slab_last
);
2137 /* Enforced rotation with a flexible axis */
2138 static void do_flexible(
2140 gmx_enfrot_t enfrot
, /* Other rotation data */
2141 t_rotgrp
*rotg
, /* The rotation group */
2142 int g
, /* Group number */
2143 rvec x
[], /* The local positions */
2145 double t
, /* Time in picoseconds */
2146 int step
, /* The time step */
2150 real sigma
; /* The Gaussian width sigma */
2151 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2154 erg
=rotg
->enfrotgrp
;
2156 /* Define the sigma value */
2157 sigma
= 0.7*rotg
->slab_dist
;
2159 /* Sort the collective coordinates erg->xc along the rotation vector. This is
2160 * an optimization for the inner loop.
2162 sort_collective_coordinates(rotg
, enfrot
->data
);
2164 /* Determine the first relevant slab for the first atom and the last
2165 * relevant slab for the last atom */
2166 get_firstlast_slab_check(rotg
, erg
, erg
->xc
[0], erg
->xc
[rotg
->nat
-1], g
, cr
);
2168 /* Determine for each slab depending on the min_gaussian cutoff criterium,
2169 * a first and a last atom index inbetween stuff needs to be calculated */
2170 get_firstlast_atom_per_slab(rotg
, cr
);
2172 /* Determine the gaussian-weighted center of positions for all slabs */
2173 get_slab_centers(rotg
,erg
->xc
,erg
->mc_sorted
,cr
,g
,t
,enfrot
->out_slabs
,bOutstep
,FALSE
);
2175 /* Clear the torque per slab from last time step: */
2176 nslabs
= erg
->slab_last
- erg
->slab_first
+ 1;
2177 for (l
=0; l
<nslabs
; l
++)
2178 erg
->slab_torque_v
[l
] = 0.0;
2180 /* Call the rotational forces kernel */
2181 GMX_MPE_LOG(ev_flexll_start
);
2182 if (rotg
->eType
== erotgFLEX
|| rotg
->eType
== erotgFLEXT
)
2183 erg
->V
= do_flex_lowlevel(rotg
, sigma
, x
, bOutstep
, box
, cr
);
2184 else if (rotg
->eType
== erotgFLEX2
|| rotg
->eType
== erotgFLEX2T
)
2185 erg
->V
= do_flex2_lowlevel(rotg
, sigma
, x
, bOutstep
, box
, cr
);
2187 gmx_fatal(FARGS
, "Unknown flexible rotation type");
2188 GMX_MPE_LOG(ev_flexll_finish
);
2190 /* Determine actual angle of this slab by RMSD fit and output to file - Let's hope */
2191 /* this only happens once in a while, since this is not parallelized! */
2192 if (bOutstep
&& MASTER(cr
))
2193 flex_fit_angle(g
, rotg
, t
, erg
->degangle
, enfrot
->out_angles
);
2197 /* Calculate the angle between reference and actual rotation group atom,
2198 * both projected into a plane perpendicular to the rotation vector: */
2199 static void angle(t_rotgrp
*rotg
,
2203 real
*weight
) /* atoms near the rotation axis should count less than atoms far away */
2205 rvec xp
, xrp
; /* current and reference positions projected on a plane perpendicular to pg->vec */
2209 /* Project x_ref and x into a plane through the origin perpendicular to rot_vec: */
2210 /* Project x_ref: xrp = x_ref - (vec * x_ref) * vec */
2211 svmul(iprod(rotg
->vec
, x_ref
), rotg
->vec
, dum
);
2212 rvec_sub(x_ref
, dum
, xrp
);
2213 /* Project x_act: */
2214 svmul(iprod(rotg
->vec
, x_act
), rotg
->vec
, dum
);
2215 rvec_sub(x_act
, dum
, xp
);
2217 /* Retrieve information about which vector precedes. gmx_angle always
2218 * returns a positive angle. */
2219 cprod(xp
, xrp
, dum
); /* if reference precedes, this is pointing into the same direction as vec */
2221 if (iprod(rotg
->vec
, dum
) >= 0)
2222 *alpha
= -gmx_angle(xrp
, xp
);
2224 *alpha
= +gmx_angle(xrp
, xp
);
2226 /* Also return the weight */
2231 /* Project first vector onto a plane perpendicular to the second vector
2233 * Note that v must be of unit length.
2235 static inline void project_onto_plane(rvec dr
, const rvec v
)
2240 svmul(iprod(dr
,v
),v
,tmp
); /* tmp = (dr.v)v */
2241 rvec_dec(dr
, tmp
); /* dr = dr - (dr.v)v */
2245 /* Fixed rotation: The rotation reference group rotates around an axis */
2246 /* The atoms of the actual rotation group are attached with imaginary */
2247 /* springs to the reference atoms. */
2248 static void do_fixed(
2250 t_rotgrp
*rotg
, /* The rotation group */
2251 rvec x
[], /* The positions */
2252 matrix box
, /* The simulation box */
2253 double t
, /* Time in picoseconds */
2254 int step
, /* The time step */
2259 rvec tmp_f
; /* Force */
2260 real alpha
; /* a single angle between an actual and a reference position */
2261 real weight
; /* single weight for a single angle */
2262 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2265 /* for mass weighting: */
2266 real wi
; /* Mass-weighting of the positions */
2268 real k_wi
; /* k times wi */
2273 erg
=rotg
->enfrotgrp
;
2274 bProject
= (rotg
->eType
==erotgPM
) || (rotg
->eType
==erotgPMPF
);
2276 /* Clear values from last time step */
2278 erg
->fix_torque_v
= 0.0;
2279 erg
->fix_angles_v
= 0.0;
2280 erg
->fix_weight_v
= 0.0;
2282 N_M
= rotg
->nat
* erg
->invmass
;
2284 /* Each process calculates the forces on its local atoms */
2285 for (i
=0; i
<erg
->nat_loc
; i
++)
2287 /* Calculate (x_i-x_c) resp. (x_i-u) */
2288 rvec_sub(erg
->x_loc_pbc
[i
], erg
->xc_center
, tmpvec
);
2290 /* Calculate Omega*(y_i-y_c)-(x_i-x_c) */
2291 rvec_sub(erg
->xr_loc
[i
], tmpvec
, dr
);
2294 project_onto_plane(dr
, rotg
->vec
);
2296 /* Mass-weighting */
2297 wi
= N_M
*erg
->m_loc
[i
];
2299 /* Store the additional force so that it can be added to the force
2300 * array after the normal forces have been evaluated */
2302 for (m
=0; m
<DIM
; m
++)
2304 tmp_f
[m
] = k_wi
*dr
[m
];
2305 erg
->f_rot_loc
[i
][m
] = tmp_f
[m
];
2306 erg
->V
+= 0.5*k_wi
*sqr(dr
[m
]);
2311 /* Add to the torque of this rotation group */
2312 erg
->fix_torque_v
+= torque(rotg
->vec
, tmp_f
, erg
->x_loc_pbc
[i
], erg
->xc_center
);
2314 /* Calculate the angle between reference and actual rotation group atom. */
2315 angle(rotg
, tmpvec
, erg
->xr_loc
[i
], &alpha
, &weight
); /* angle in rad, weighted */
2316 erg
->fix_angles_v
+= alpha
* weight
;
2317 erg
->fix_weight_v
+= weight
;
2319 /* If you want enforced rotation to contribute to the virial,
2320 * activate the following lines:
2323 Add the rotation contribution to the virial
2324 for(j=0; j<DIM; j++)
2326 vir[j][m] += 0.5*f[ii][j]*dr[m];
2330 fprintf(stderr
,"step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step
, cr
->nodeid
,
2331 erg
->xc_ref_ind
[i
],erg
->f_rot_loc
[i
][XX
], erg
->f_rot_loc
[i
][YY
], erg
->f_rot_loc
[i
][ZZ
],erg
->fix_torque_v
);
2333 } /* end of loop over local rotation group atoms */
2337 /* Calculate the radial motion potential and forces */
2338 static void do_radial_motion(
2340 t_rotgrp
*rotg
, /* The rotation group */
2341 rvec x
[], /* The positions */
2342 matrix box
, /* The simulation box */
2343 double t
, /* Time in picoseconds */
2344 int step
, /* The time step */
2348 rvec tmp_f
; /* Force */
2349 real alpha
; /* a single angle between an actual and a reference position */
2350 real weight
; /* single weight for a single angle */
2351 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2352 rvec xj_u
; /* xj - u */
2357 /* For mass weighting: */
2358 real wj
; /* Mass-weighting of the positions */
2362 erg
=rotg
->enfrotgrp
;
2364 /* Clear values from last time step */
2367 erg
->fix_torque_v
= 0.0;
2368 erg
->fix_angles_v
= 0.0;
2369 erg
->fix_weight_v
= 0.0;
2371 N_M
= rotg
->nat
* erg
->invmass
;
2373 /* Each process calculates the forces on its local atoms */
2374 for (j
=0; j
<erg
->nat_loc
; j
++)
2376 /* Calculate (xj-u) */
2377 rvec_sub(erg
->x_loc_pbc
[j
], erg
->xc_center
, xj_u
); /* xj_u = xj-u */
2379 /* Calculate Omega.(yj-u) */
2380 cprod(rotg
->vec
, erg
->xr_loc
[j
], tmpvec
); /* tmpvec = v x Omega.(yj-u) */
2382 /* v x Omega.(yj-u) */
2383 unitv(tmpvec
, pj
); /* pj = -------------------- */
2384 /* | v x Omega.(yj-u) | */
2386 fac
= iprod(pj
, xj_u
); /* fac = pj.(xj-u) */
2389 /* Mass-weighting */
2390 wj
= N_M
*erg
->m_loc
[j
];
2392 /* Store the additional force so that it can be added to the force
2393 * array after the normal forces have been evaluated */
2394 svmul(-rotg
->k
*wj
*fac
, pj
, tmp_f
);
2395 copy_rvec(tmp_f
, erg
->f_rot_loc
[j
]);
2399 /* Add to the torque of this rotation group */
2400 erg
->fix_torque_v
+= torque(rotg
->vec
, tmp_f
, erg
->x_loc_pbc
[j
], erg
->xc_center
);
2402 /* Calculate the angle between reference and actual rotation group atom. */
2403 angle(rotg
, xj_u
, erg
->xr_loc
[j
], &alpha
, &weight
); /* angle in rad, weighted */
2404 erg
->fix_angles_v
+= alpha
* weight
;
2405 erg
->fix_weight_v
+= weight
;
2408 fprintf(stderr
,"RM: step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step
, cr
->nodeid
,
2409 erg
->xc_ref_ind
[j
],erg
->f_rot_loc
[j
][XX
], erg
->f_rot_loc
[j
][YY
], erg
->f_rot_loc
[j
][ZZ
],erg
->fix_torque_v
);
2411 } /* end of loop over local rotation group atoms */
2412 erg
->V
= 0.5*rotg
->k
*sum
;
2416 /* Calculate the radial motion pivot-free potential and forces */
2417 static void do_radial_motion_pf(
2419 t_rotgrp
*rotg
, /* The rotation group */
2420 rvec x
[], /* The positions */
2421 matrix box
, /* The simulation box */
2422 double t
, /* Time in picoseconds */
2423 int step
, /* The time step */
2427 rvec xj
; /* Current position */
2428 rvec xj_xc
; /* xj - xc */
2429 rvec yj0_yc0
; /* yj0 - yc0 */
2430 rvec tmp_f
; /* Force */
2431 real alpha
; /* a single angle between an actual and a reference position */
2432 real weight
; /* single weight for a single angle */
2433 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2434 rvec tmpvec
, tmpvec2
;
2435 rvec innersumvec
; /* Precalculation of the inner sum */
2440 /* For mass weighting: */
2441 real mj
,wi
,wj
; /* Mass-weighting of the positions */
2445 erg
=rotg
->enfrotgrp
;
2447 /* Clear values from last time step */
2450 erg
->fix_torque_v
= 0.0;
2451 erg
->fix_angles_v
= 0.0;
2452 erg
->fix_weight_v
= 0.0;
2454 N_M
= rotg
->nat
* erg
->invmass
;
2456 /* Get the current center of the rotation group: */
2457 get_center(erg
->xc
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
2459 /* Precalculate Sum_i [ wi qi.(xi-xc) qi ] which is needed for every single j */
2460 clear_rvec(innersumvec
);
2461 for (i
=0; i
< rotg
->nat
; i
++)
2463 /* Mass-weighting */
2464 wi
= N_M
*erg
->mc
[i
];
2466 /* Calculate qi. Note that xc_ref_center has already been subtracted from
2467 * x_ref in init_rot_group.*/
2468 mvmul(erg
->rotmat
, rotg
->x_ref
[i
], tmpvec
); /* tmpvec = Omega.(yi0-yc0) */
2470 cprod(rotg
->vec
, tmpvec
, tmpvec2
); /* tmpvec2 = v x Omega.(yi0-yc0) */
2472 /* v x Omega.(yi0-yc0) */
2473 unitv(tmpvec2
, qi
); /* qi = ----------------------- */
2474 /* | v x Omega.(yi0-yc0) | */
2476 rvec_sub(erg
->xc
[i
], erg
->xc_center
, tmpvec
); /* tmpvec = xi-xc */
2478 svmul(wi
*iprod(qi
, tmpvec
), qi
, tmpvec2
);
2480 rvec_inc(innersumvec
, tmpvec2
);
2482 svmul(rotg
->k
*erg
->invmass
, innersumvec
, innersumveckM
);
2484 /* Each process calculates the forces on its local atoms */
2485 for (j
=0; j
<erg
->nat_loc
; j
++)
2487 /* Local index of a rotation group atom */
2488 ii
= erg
->ind_loc
[j
];
2489 /* Position of this atom in the collective array */
2490 iigrp
= erg
->xc_ref_ind
[j
];
2491 /* Mass-weighting */
2492 mj
= erg
->mc
[iigrp
]; /* need the unsorted mass here */
2495 /* Current position of this atom: x[ii][XX/YY/ZZ] */
2496 copy_rvec(x
[ii
], xj
);
2498 /* Shift this atom such that it is near its reference */
2499 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
2501 /* The (unrotated) reference position is yj0. yc0 has already
2502 * been subtracted in init_rot_group */
2503 copy_rvec(rotg
->x_ref
[iigrp
], yj0_yc0
); /* yj0_yc0 = yj0 - yc0 */
2505 /* Calculate Omega.(yj0-yc0) */
2506 mvmul(erg
->rotmat
, yj0_yc0
, tmpvec2
); /* tmpvec2 = Omega.(yj0 - yc0) */
2508 cprod(rotg
->vec
, tmpvec2
, tmpvec
); /* tmpvec = v x Omega.(yj0-yc0) */
2510 /* v x Omega.(yj0-yc0) */
2511 unitv(tmpvec
, qj
); /* qj = ----------------------- */
2512 /* | v x Omega.(yj0-yc0) | */
2514 /* Calculate (xj-xc) */
2515 rvec_sub(xj
, erg
->xc_center
, xj_xc
); /* xj_xc = xj-xc */
2517 fac
= iprod(qj
, xj_xc
); /* fac = qj.(xj-xc) */
2520 /* Store the additional force so that it can be added to the force
2521 * array after the normal forces have been evaluated */
2522 svmul(-rotg
->k
*wj
*fac
, qj
, tmp_f
); /* part 1 of force */
2523 svmul(mj
, innersumveckM
, tmpvec
); /* part 2 of force */
2524 rvec_inc(tmp_f
, tmpvec
);
2525 copy_rvec(tmp_f
, erg
->f_rot_loc
[j
]);
2529 /* Add to the torque of this rotation group */
2530 erg
->fix_torque_v
+= torque(rotg
->vec
, tmp_f
, xj
, erg
->xc_center
);
2532 /* Calculate the angle between reference and actual rotation group atom. */
2533 angle(rotg
, xj_xc
, yj0_yc0
, &alpha
, &weight
); /* angle in rad, weighted */
2534 erg
->fix_angles_v
+= alpha
* weight
;
2535 erg
->fix_weight_v
+= weight
;
2538 fprintf(stderr
,"RM-PF: step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step
, cr
->nodeid
,
2539 erg
->xc_ref_ind
[j
],erg
->f_rot_loc
[j
][XX
], erg
->f_rot_loc
[j
][YY
], erg
->f_rot_loc
[j
][ZZ
],erg
->fix_torque_v
);
2541 } /* end of loop over local rotation group atoms */
2542 erg
->V
= 0.5*rotg
->k
*V
;
2546 /* Precalculate the inner sum for the radial motion 2 forces */
2547 static void radial_motion2_precalc_inner_sum(t_rotgrp
*rotg
, rvec innersumvec
)
2550 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2551 rvec xi_xc
; /* xj - xc */
2552 rvec tmpvec
,tmpvec2
;
2556 rvec v_xi_xc
; /* v x (xj - u) */
2558 real wi
; /* Mass-weighting of the positions */
2562 erg
=rotg
->enfrotgrp
;
2563 N_M
= rotg
->nat
* erg
->invmass
;
2565 /* Loop over the collective set of positions */
2567 for (i
=0; i
<rotg
->nat
; i
++)
2569 /* Mass-weighting */
2570 wi
= N_M
*erg
->mc
[i
];
2572 rvec_sub(erg
->xc
[i
], erg
->xc_center
, xi_xc
); /* xi_xc = xi-xc */
2574 /* Calculate ri. Note that xc_ref_center has already been subtracted from
2575 * x_ref in init_rot_group.*/
2576 mvmul(erg
->rotmat
, rotg
->x_ref
[i
], ri
); /* ri = Omega.(yi0-yc0) */
2578 cprod(rotg
->vec
, xi_xc
, v_xi_xc
); /* v_xi_xc = v x (xi-u) */
2580 fac
= norm2(v_xi_xc
);
2582 psiistar
= 1.0/(fac
+ rotg
->eps
); /* psiistar = --------------------- */
2583 /* |v x (xi-xc)|^2 + eps */
2585 psii
= gmx_invsqrt(fac
); /* 1 */
2586 /* psii = ------------- */
2589 svmul(psii
, v_xi_xc
, si
); /* si = psii * (v x (xi-xc) ) */
2591 fac
= iprod(v_xi_xc
, ri
); /* fac = (v x (xi-xc)).ri */
2594 siri
= iprod(si
, ri
); /* siri = si.ri */
2596 svmul(psiistar
/psii
, ri
, tmpvec
);
2597 svmul(psiistar
*psiistar
/(psii
*psii
*psii
) * siri
, si
, tmpvec2
);
2598 rvec_dec(tmpvec
, tmpvec2
);
2599 cprod(tmpvec
, rotg
->vec
, tmpvec2
);
2601 svmul(wi
*siri
, tmpvec2
, tmpvec
);
2603 rvec_inc(sumvec
, tmpvec
);
2605 svmul(rotg
->k
*erg
->invmass
, sumvec
, innersumvec
);
2609 /* Calculate the radial motion 2 potential and forces */
2610 static void do_radial_motion2(
2612 t_rotgrp
*rotg
, /* The rotation group */
2613 rvec x
[], /* The positions */
2614 matrix box
, /* The simulation box */
2615 double t
, /* Time in picoseconds */
2616 int step
, /* The time step */
2620 rvec xj
; /* Position */
2621 real alpha
; /* a single angle between an actual and a reference position */
2622 real weight
; /* single weight for a single angle */
2623 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2624 rvec xj_u
; /* xj - u */
2625 rvec tmpvec
,tmpvec2
;
2626 real fac
,fac2
,Vpart
;
2629 rvec v_xj_u
; /* v x (xj - u) */
2631 real mj
,wj
; /* For mass-weighting of the positions */
2637 erg
=rotg
->enfrotgrp
;
2639 bPF
= rotg
->eType
==erotgRM2PF
;
2640 clear_rvec(innersumvec
);
2643 /* For the pivot-free variant we have to use the current center of
2644 * mass of the rotation group instead of the pivot u */
2645 get_center(erg
->xc
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
2647 /* Also, we precalculate the second term of the forces that is identical
2648 * (up to the weight factor mj) for all forces */
2649 radial_motion2_precalc_inner_sum(rotg
,innersumvec
);
2652 /* Clear values from last time step */
2655 erg
->fix_torque_v
= 0.0;
2656 erg
->fix_angles_v
= 0.0;
2657 erg
->fix_weight_v
= 0.0;
2659 N_M
= rotg
->nat
* erg
->invmass
;
2661 /* Each process calculates the forces on its local atoms */
2662 for (j
=0; j
<erg
->nat_loc
; j
++)
2666 /* Local index of a rotation group atom */
2667 ii
= erg
->ind_loc
[j
];
2668 /* Position of this atom in the collective array */
2669 iigrp
= erg
->xc_ref_ind
[j
];
2670 /* Mass-weighting */
2671 mj
= erg
->mc
[iigrp
];
2673 /* Current position of this atom: x[ii] */
2674 copy_rvec(x
[ii
], xj
);
2676 /* Shift this atom such that it is near its reference */
2677 shift_single_coord(box
, xj
, erg
->xc_shifts
[iigrp
]);
2679 /* The (unrotated) reference position is yj0. yc0 has already
2680 * been subtracted in init_rot_group */
2681 copy_rvec(rotg
->x_ref
[iigrp
], tmpvec
); /* tmpvec = yj0 - yc0 */
2683 /* Calculate Omega.(yj0-yc0) */
2684 mvmul(erg
->rotmat
, tmpvec
, rj
); /* rj = Omega.(yj0-yc0) */
2689 copy_rvec(erg
->x_loc_pbc
[j
], xj
);
2690 copy_rvec(erg
->xr_loc
[j
], rj
); /* rj = Omega.(yj0-u) */
2692 /* Mass-weighting */
2695 /* Calculate (xj-u) resp. (xj-xc) */
2696 rvec_sub(xj
, erg
->xc_center
, xj_u
); /* xj_u = xj-u */
2698 cprod(rotg
->vec
, xj_u
, v_xj_u
); /* v_xj_u = v x (xj-u) */
2700 fac
= norm2(v_xj_u
);
2702 psijstar
= 1.0/(fac
+ rotg
->eps
); /* psistar = -------------------- */
2703 /* |v x (xj-u)|^2 + eps */
2705 psij
= gmx_invsqrt(fac
); /* 1 */
2706 /* psij = ------------ */
2709 svmul(psij
, v_xj_u
, sj
); /* sj = psij * (v x (xj-u) ) */
2711 fac
= iprod(v_xj_u
, rj
); /* fac = (v x (xj-u)).rj */
2714 sjrj
= iprod(sj
, rj
); /* sjrj = sj.rj */
2716 svmul(psijstar
/psij
, rj
, tmpvec
);
2717 svmul(psijstar
*psijstar
/(psij
*psij
*psij
) * sjrj
, sj
, tmpvec2
);
2718 rvec_dec(tmpvec
, tmpvec2
);
2719 cprod(tmpvec
, rotg
->vec
, tmpvec2
);
2721 /* Store the additional force so that it can be added to the force
2722 * array after the normal forces have been evaluated */
2723 svmul(-rotg
->k
*wj
*sjrj
, tmpvec2
, tmpvec
);
2724 svmul(mj
, innersumvec
, tmpvec2
); /* This is != 0 only for the pivot-free variant */
2726 rvec_add(tmpvec2
, tmpvec
, erg
->f_rot_loc
[j
]);
2727 Vpart
+= wj
*psijstar
*fac2
;
2730 /* Add to the torque of this rotation group */
2731 erg
->fix_torque_v
+= torque(rotg
->vec
, erg
->f_rot_loc
[j
], xj
, erg
->xc_center
);
2733 /* Calculate the angle between reference and actual rotation group atom. */
2734 angle(rotg
, xj_u
, rj
, &alpha
, &weight
); /* angle in rad, weighted */
2735 erg
->fix_angles_v
+= alpha
* weight
;
2736 erg
->fix_weight_v
+= weight
;
2739 fprintf(stderr
,"RM2: step %d node%d FORCE on ATOM %d = (%15.8f %15.8f %15.8f) torque=%15.8f\n", step
, cr
->nodeid
,
2740 erg
->xc_ref_ind
[j
],erg
->f_rot_loc
[j
][XX
], erg
->f_rot_loc
[j
][YY
], erg
->f_rot_loc
[j
][ZZ
],erg
->fix_torque_v
);
2742 } /* end of loop over local rotation group atoms */
2743 erg
->V
= 0.5*rotg
->k
*Vpart
;
2747 /* Determine the smallest and largest position vector (with respect to the
2748 * rotation vector) for the reference group */
2749 static void get_firstlast_atom_ref(
2754 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2756 real xcproj
; /* The projection of a reference position on the
2758 real minproj
, maxproj
; /* Smallest and largest projection on v */
2762 erg
=rotg
->enfrotgrp
;
2764 /* Start with some value */
2765 minproj
= iprod(rotg
->x_ref
[0], rotg
->vec
);
2768 /* This is just to ensure that it still works if all the atoms of the
2769 * reference structure are situated in a plane perpendicular to the rotation
2772 *lastindex
= rotg
->nat
-1;
2774 /* Loop over all atoms of the reference group,
2775 * project them on the rotation vector to find the extremes */
2776 for (i
=0; i
<rotg
->nat
; i
++)
2778 xcproj
= iprod(rotg
->x_ref
[i
], rotg
->vec
);
2779 if (xcproj
< minproj
)
2784 if (xcproj
> maxproj
)
2793 /* Allocate memory for the slabs */
2794 static void allocate_slabs(
2801 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2805 erg
=rotg
->enfrotgrp
;
2807 /* More slabs than are defined for the reference are never needed */
2808 nslabs
= erg
->slab_last_ref
- erg
->slab_first_ref
+ 1;
2810 /* Remember how many we allocated */
2811 erg
->nslabs_alloc
= nslabs
;
2813 if (MASTER(cr
) && bVerbose
)
2814 fprintf(fplog
, "%s allocating memory to store data for %d slabs (rotation group %d).\n",
2816 snew(erg
->slab_center
, nslabs
);
2817 snew(erg
->slab_center_ref
, nslabs
);
2818 snew(erg
->slab_weights
, nslabs
);
2819 snew(erg
->slab_torque_v
, nslabs
);
2820 snew(erg
->slab_data
, nslabs
);
2821 snew(erg
->gn_atom
, nslabs
);
2822 snew(erg
->gn_slabind
, nslabs
);
2823 snew(erg
->slab_innersumvec
, nslabs
);
2824 for (i
=0; i
<nslabs
; i
++)
2826 snew(erg
->slab_data
[i
].x
, rotg
->nat
);
2827 snew(erg
->slab_data
[i
].ref
, rotg
->nat
);
2828 snew(erg
->slab_data
[i
].weight
, rotg
->nat
);
2830 snew(erg
->xc_ref_sorted
, rotg
->nat
);
2831 snew(erg
->xc_sortind
, rotg
->nat
);
2832 snew(erg
->firstatom
, nslabs
);
2833 snew(erg
->lastatom
, nslabs
);
2837 /* From the extreme coordinates of the reference group, determine the first
2838 * and last slab of the reference. We can never have more slabs in the real
2839 * simulation than calculated here for the reference.
2841 static void get_firstlast_slab_ref(t_rotgrp
*rotg
, real mc
[], int ref_firstindex
, int ref_lastindex
)
2843 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2844 int first
,last
,firststart
;
2848 erg
=rotg
->enfrotgrp
;
2849 first
= get_first_slab(rotg
, erg
->max_beta
, rotg
->x_ref
[ref_firstindex
]);
2850 last
= get_last_slab( rotg
, erg
->max_beta
, rotg
->x_ref
[ref_lastindex
]);
2853 while (get_slab_weight(first
, rotg
, rotg
->x_ref
, mc
, &dummy
) > WEIGHT_MIN
)
2857 erg
->slab_first_ref
= first
+1;
2858 while (get_slab_weight(last
, rotg
, rotg
->x_ref
, mc
, &dummy
) > WEIGHT_MIN
)
2862 erg
->slab_last_ref
= last
-1;
2864 erg
->slab_buffer
= firststart
- erg
->slab_first_ref
;
2869 static void init_rot_group(FILE *fplog
,t_commrec
*cr
,int g
,t_rotgrp
*rotg
,
2870 rvec
*x
,gmx_mtop_t
*mtop
,gmx_bool bVerbose
,FILE *out_slabs
)
2874 gmx_bool bFlex
,bColl
;
2876 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
2877 int ref_firstindex
, ref_lastindex
;
2878 real mass
,totalmass
;
2881 /* Do we have a flexible axis? */
2882 bFlex
= ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
2883 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) );
2885 /* Do we use a global set of coordinates? */
2886 bColl
= bFlex
|| (rotg
->eType
==erotgRMPF
) || (rotg
->eType
==erotgRM2PF
);
2888 erg
=rotg
->enfrotgrp
;
2890 /* Allocate space for collective coordinates if needed */
2893 snew(erg
->xc
, rotg
->nat
);
2894 snew(erg
->xc_shifts
, rotg
->nat
);
2895 snew(erg
->xc_eshifts
, rotg
->nat
);
2897 /* Save the original (whole) set of positions such that later the
2898 * molecule can always be made whole again */
2899 snew(erg
->xc_old
, rotg
->nat
);
2902 for (i
=0; i
<rotg
->nat
; i
++)
2905 copy_rvec(x
[ii
], erg
->xc_old
[i
]);
2910 gmx_bcast(rotg
->nat
*sizeof(erg
->xc_old
[0]),erg
->xc_old
, cr
);
2913 if (rotg
->eFittype
== erotgFitNORM
)
2915 snew(erg
->xc_ref_length
, rotg
->nat
); /* in case fit type NORM is chosen */
2916 snew(erg
->xc_norm
, rotg
->nat
);
2921 snew(erg
->xr_loc
, rotg
->nat
);
2922 snew(erg
->x_loc_pbc
, rotg
->nat
);
2925 snew(erg
->f_rot_loc
, rotg
->nat
);
2926 snew(erg
->xc_ref_ind
, rotg
->nat
);
2928 /* xc_ref_ind needs to be set to identity in the serial case */
2930 for (i
=0; i
<rotg
->nat
; i
++)
2931 erg
->xc_ref_ind
[i
] = i
;
2933 /* Copy the masses so that the COM can be determined. For all types of
2934 * enforced rotation, we store the masses in the erg->mc array. */
2935 snew(erg
->mc
, rotg
->nat
);
2937 snew(erg
->mc_sorted
, rotg
->nat
);
2939 snew(erg
->m_loc
, rotg
->nat
);
2941 for (i
=0; i
<rotg
->nat
; i
++)
2945 gmx_mtop_atomnr_to_atom(mtop
,rotg
->ind
[i
],&atom
);
2955 erg
->invmass
= 1.0/totalmass
;
2957 /* Set xc_ref_center for any rotation potential */
2958 if ((rotg
->eType
==erotgISO
) || (rotg
->eType
==erotgPM
) || (rotg
->eType
==erotgRM
) || (rotg
->eType
==erotgRM2
))
2960 /* Set the pivot point for the fixed, stationary-axis potentials. This
2961 * won't change during the simulation */
2962 copy_rvec(rotg
->pivot
, erg
->xc_ref_center
);
2963 copy_rvec(rotg
->pivot
, erg
->xc_center
);
2967 /* Center of the reference positions */
2968 get_center(rotg
->x_ref
, erg
->mc
, rotg
->nat
, erg
->xc_ref_center
);
2970 /* Center of the actual positions */
2973 snew(xdum
, rotg
->nat
);
2974 for (i
=0; i
<rotg
->nat
; i
++)
2977 copy_rvec(x
[ii
], xdum
[i
]);
2979 get_center(xdum
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
2984 gmx_bcast(sizeof(erg
->xc_center
), erg
->xc_center
, cr
);
2988 if ( (rotg
->eType
!= erotgFLEX
) && (rotg
->eType
!= erotgFLEX2
) )
2990 /* Put the reference positions into origin: */
2991 for (i
=0; i
<rotg
->nat
; i
++)
2992 rvec_dec(rotg
->x_ref
[i
], erg
->xc_ref_center
);
2995 /* Enforced rotation with flexible axis */
2998 /* Calculate maximum beta value from minimum gaussian (performance opt.) */
2999 erg
->max_beta
= calc_beta_max(rotg
->min_gaussian
, rotg
->slab_dist
);
3001 /* Determine the smallest and largest coordinate with respect to the rotation vector */
3002 get_firstlast_atom_ref(rotg
, &ref_firstindex
, &ref_lastindex
);
3004 /* From the extreme coordinates of the reference group, determine the first
3005 * and last slab of the reference. */
3006 get_firstlast_slab_ref(rotg
, erg
->mc
, ref_firstindex
, ref_lastindex
);
3008 /* Allocate memory for the slabs */
3009 allocate_slabs(rotg
, fplog
, g
, bVerbose
, cr
);
3011 /* Flexible rotation: determine the reference centers for the rest of the simulation */
3012 erg
->slab_first
= erg
->slab_first_ref
;
3013 erg
->slab_last
= erg
->slab_last_ref
;
3014 get_slab_centers(rotg
,rotg
->x_ref
,erg
->mc
,cr
,g
,-1,out_slabs
,TRUE
,TRUE
);
3016 /* Length of each x_rotref vector from center (needed if fit routine NORM is chosen): */
3017 if (rotg
->eFittype
== erotgFitNORM
)
3019 for (i
=0; i
<rotg
->nat
; i
++)
3021 rvec_sub(rotg
->x_ref
[i
], erg
->xc_ref_center
, coord
);
3022 erg
->xc_ref_length
[i
] = norm(coord
);
3029 extern void dd_make_local_rotation_groups(gmx_domdec_t
*dd
,t_rot
*rot
)
3034 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3038 for(g
=0; g
<rot
->ngrp
; g
++)
3040 rotg
= &rot
->grp
[g
];
3041 erg
= rotg
->enfrotgrp
;
3044 dd_make_local_group_indices(ga2la
,rotg
->nat
,rotg
->ind
,
3045 &erg
->nat_loc
,&erg
->ind_loc
,&erg
->nalloc_loc
,erg
->xc_ref_ind
);
3050 extern void init_rot(FILE *fplog
,t_inputrec
*ir
,int nfile
,const t_filenm fnm
[],
3051 t_commrec
*cr
, rvec
*x
, matrix box
, gmx_mtop_t
*mtop
, const output_env_t oenv
,
3052 gmx_bool bVerbose
, unsigned long Flags
)
3057 int nat_max
=0; /* Size of biggest rotation group */
3058 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
3059 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3060 rvec
*x_pbc
; /* Space for the pbc-correct atom positions */
3063 if ( (PAR(cr
)) && !DOMAINDECOMP(cr
) )
3064 gmx_fatal(FARGS
, "Enforced rotation is only implemented for domain decomposition!");
3066 if ( MASTER(cr
) && bVerbose
)
3067 fprintf(stdout
, "%s Initializing ...\n", RotStr
);
3071 snew(rot
->enfrot
, 1);
3074 /* Output every step for reruns */
3075 if (Flags
& MD_RERUN
)
3078 fprintf(fplog
, "%s rerun - will write rotation output every available step.\n", RotStr
);
3083 er
->out_slabs
= NULL
;
3085 er
->out_slabs
= open_slab_out(rot
, opt2fn("-rs",nfile
,fnm
));
3089 /* Remove pbc, make molecule whole.
3090 * When ir->bContinuation=TRUE this has already been done, but ok. */
3091 snew(x_pbc
,mtop
->natoms
);
3092 m_rveccopy(mtop
->natoms
,x
,x_pbc
);
3093 do_pbc_first_mtop(NULL
,ir
->ePBC
,box
,mtop
,x_pbc
);
3096 for(g
=0; g
<rot
->ngrp
; g
++)
3098 rotg
= &rot
->grp
[g
];
3101 fprintf(fplog
,"%s group %d type '%s'\n", RotStr
, g
, erotg_names
[rotg
->eType
]);
3105 /* Allocate space for the rotation group's data: */
3106 snew(rotg
->enfrotgrp
, 1);
3107 erg
= rotg
->enfrotgrp
;
3109 nat_max
=max(nat_max
, rotg
->nat
);
3114 erg
->nalloc_loc
= 0;
3115 erg
->ind_loc
= NULL
;
3119 erg
->nat_loc
= rotg
->nat
;
3120 erg
->ind_loc
= rotg
->ind
;
3122 init_rot_group(fplog
,cr
,g
,rotg
,x_pbc
,mtop
,bVerbose
,er
->out_slabs
);
3126 /* Allocate space for enforced rotation buffer variables */
3127 er
->bufsize
= nat_max
;
3128 snew(er
->data
, nat_max
);
3129 snew(er
->xbuf
, nat_max
);
3130 snew(er
->mbuf
, nat_max
);
3132 /* Buffers for MPI reducing torques, angles, weights (for each group), and V */
3133 er
->mpi_bufsize
= 4*rot
->ngrp
; /* To start with */
3134 snew(er
->mpi_inbuf
, er
->mpi_bufsize
);
3135 snew(er
->mpi_outbuf
, er
->mpi_bufsize
);
3137 /* Only do I/O on the MASTER */
3138 er
->out_angles
= NULL
;
3140 er
->out_torque
= NULL
;
3143 er
->out_rot
= open_rot_out(opt2fn("-ro",nfile
,fnm
), rot
, oenv
, Flags
);
3144 if ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
3145 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) )
3147 if (rot
->nstrout
> 0)
3148 er
->out_angles
= open_angles_out(rot
, opt2fn("-ra",nfile
,fnm
));
3149 if (rot
->nsttout
> 0)
3150 er
->out_torque
= open_torque_out(rot
, opt2fn("-rt",nfile
,fnm
));
3157 extern void finish_rot(FILE *fplog
,t_rot
*rot
)
3159 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
3164 gmx_fio_fclose(er
->out_rot
);
3166 gmx_fio_fclose(er
->out_slabs
);
3168 gmx_fio_fclose(er
->out_angles
);
3170 gmx_fio_fclose(er
->out_torque
);
3174 /* Rotate the local reference positions and store them in
3175 * erg->xr_loc[0...(nat_loc-1)]
3177 * Note that we already subtracted u or y_c from the reference positions
3178 * in init_rot_group().
3180 static void rotate_local_reference(t_rotgrp
*rotg
)
3182 gmx_enfrotgrp_t erg
;
3186 erg
=rotg
->enfrotgrp
;
3188 for (i
=0; i
<erg
->nat_loc
; i
++)
3190 /* Index of this rotation group atom with respect to the whole rotation group */
3191 ii
= erg
->xc_ref_ind
[i
];
3193 mvmul(erg
->rotmat
, rotg
->x_ref
[ii
], erg
->xr_loc
[i
]);
3198 /* Select the PBC representation for each local x position and store that
3199 * for later usage. We assume the right PBC image of an x is the one nearest to
3200 * its rotated reference */
3201 static void choose_pbc_image(rvec x
[], t_rotgrp
*rotg
, matrix box
, int npbcdim
)
3204 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3209 erg
=rotg
->enfrotgrp
;
3211 for (i
=0; i
<erg
->nat_loc
; i
++)
3215 /* Index of a rotation group atom */
3216 ii
= erg
->ind_loc
[i
];
3218 /* Get the reference position. The pivot (or COM or COG) was already
3219 * subtracted in init_rot_group() from the reference positions. Also,
3220 * the reference positions have already been rotated in
3221 * rotate_local_reference() */
3222 copy_rvec(erg
->xr_loc
[i
], xref
);
3224 /* Subtract the (old) center from the current positions
3225 * (just to determine the shifts!) */
3226 rvec_sub(x
[ii
], erg
->xc_center
, xcurr
);
3228 /* Shortest PBC distance between the atom and its reference */
3229 rvec_sub(xcurr
, xref
, dx
);
3231 /* Determine the shift for this atom */
3232 for(m
=npbcdim
-1; m
>=0; m
--)
3234 while (dx
[m
] < -0.5*box
[m
][m
])
3236 for(d
=0; d
<DIM
; d
++)
3240 while (dx
[m
] >= 0.5*box
[m
][m
])
3242 for(d
=0; d
<DIM
; d
++)
3248 /* Apply the shift to the current atom */
3249 copy_rvec(x
[ii
], erg
->x_loc_pbc
[i
]);
3250 shift_single_coord(box
, erg
->x_loc_pbc
[i
], shift
);
3255 extern void do_rotation(
3262 gmx_wallcycle_t wcycle
,
3268 gmx_bool outstep_torque
;
3269 gmx_bool bFlex
,bColl
;
3271 gmx_enfrot_t er
; /* Pointer to the enforced rotation buffer variables */
3272 gmx_enfrotgrp_t erg
; /* Pointer to enforced rotation group data */
3282 /* At which time steps do we want to output the torque */
3283 outstep_torque
= do_per_step(step
, rot
->nsttout
);
3285 /* Output time into rotation output file */
3286 if (outstep_torque
&& MASTER(cr
))
3287 fprintf(er
->out_rot
, "%12.3e",t
);
3289 /**************************************************************************/
3290 /* First do ALL the communication! */
3291 for(g
=0; g
<rot
->ngrp
; g
++)
3293 rotg
= &rot
->grp
[g
];
3294 erg
=rotg
->enfrotgrp
;
3296 /* Do we have a flexible axis? */
3297 bFlex
= ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
3298 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) );
3300 /* Do we use a collective (global) set of coordinates? */
3301 bColl
= bFlex
|| (rotg
->eType
==erotgRMPF
) || (rotg
->eType
==erotgRM2PF
);
3303 /* Calculate the rotation matrix for this angle: */
3304 erg
->degangle
= rotg
->rate
* t
;
3305 calc_rotmat(rotg
->vec
,erg
->degangle
,erg
->rotmat
);
3309 /* Transfer the rotation group's positions such that every node has
3310 * all of them. Every node contributes its local positions x and stores
3311 * it in the collective erg->xc array. */
3312 communicate_group_positions(cr
,erg
->xc
, erg
->xc_shifts
, erg
->xc_eshifts
, bNS
,
3313 x
, rotg
->nat
, erg
->nat_loc
, erg
->ind_loc
, erg
->xc_ref_ind
, erg
->xc_old
, box
);
3317 /* Fill the local masses array;
3318 * this array changes in DD/neighborsearching steps */
3321 for (i
=0; i
<erg
->nat_loc
; i
++)
3323 /* Index of local atom w.r.t. the collective rotation group */
3324 ii
= erg
->xc_ref_ind
[i
];
3325 erg
->m_loc
[i
] = erg
->mc
[ii
];
3329 /* Calculate Omega*(y_i-y_c) for the local positions */
3330 rotate_local_reference(rotg
);
3332 /* Choose the nearest PBC images of the group atoms with respect
3333 * to the rotated reference positions */
3334 choose_pbc_image(x
, rotg
, box
, 3);
3336 /* Get the center of the rotation group */
3337 if ( (rotg
->eType
==erotgISOPF
) || (rotg
->eType
==erotgPMPF
) )
3338 get_center_comm(cr
, erg
->x_loc_pbc
, erg
->m_loc
, erg
->nat_loc
, rotg
->nat
, erg
->xc_center
);
3341 } /* End of loop over rotation groups */
3343 /**************************************************************************/
3344 /* Done communicating, we can start to count cycles now ... */
3345 wallcycle_start(wcycle
, ewcROT
);
3346 GMX_MPE_LOG(ev_rotcycles_start
);
3352 for(g
=0; g
<rot
->ngrp
; g
++)
3354 rotg
= &rot
->grp
[g
];
3355 erg
=rotg
->enfrotgrp
;
3357 bFlex
= ( (rotg
->eType
==erotgFLEX
) || (rotg
->eType
==erotgFLEXT
)
3358 || (rotg
->eType
==erotgFLEX2
) || (rotg
->eType
==erotgFLEX2T
) );
3360 bColl
= bFlex
|| (rotg
->eType
==erotgRMPF
) || (rotg
->eType
==erotgRM2PF
);
3362 if (outstep_torque
&& MASTER(cr
))
3363 fprintf(er
->out_rot
, "%12.4f", erg
->degangle
);
3371 do_fixed(cr
,rotg
,x
,box
,t
,step
,outstep_torque
);
3374 do_radial_motion(cr
,rotg
,x
,box
,t
,step
,outstep_torque
);
3377 do_radial_motion_pf(cr
,rotg
,x
,box
,t
,step
,outstep_torque
);
3381 do_radial_motion2(cr
,rotg
,x
,box
,t
,step
,outstep_torque
);
3385 /* Subtract the center of the rotation group from the collective positions array
3386 * Also store the center in erg->xc_center since it needs to be subtracted
3387 * in the low level routines from the local coordinates as well */
3388 get_center(erg
->xc
, erg
->mc
, rotg
->nat
, erg
->xc_center
);
3389 svmul(-1.0, erg
->xc_center
, transvec
);
3390 translate_x(erg
->xc
, rotg
->nat
, transvec
);
3391 do_flexible(cr
,er
,rotg
,g
,x
,box
,t
,step
,outstep_torque
);
3395 /* Do NOT subtract the center of mass in the low level routines! */
3396 clear_rvec(erg
->xc_center
);
3397 do_flexible(cr
,er
,rotg
,g
,x
,box
,t
,step
,outstep_torque
);
3400 gmx_fatal(FARGS
, "No such rotation potential.");
3407 fprintf(stderr
, "%s calculation (step %d) took %g seconds.\n", RotStr
, step
, MPI_Wtime()-t0
);
3410 /* Stop the cycle counter and add to the force cycles for load balancing */
3411 cycles_rot
= wallcycle_stop(wcycle
,ewcROT
);
3412 if (DOMAINDECOMP(cr
) && wcycle
)
3413 dd_cycles_add(cr
->dd
,cycles_rot
,ddCyclF
);
3414 GMX_MPE_LOG(ev_rotcycles_finish
);