File: | gromacs/essentialdynamics/edsam.c |
Location: | line 2202, column 9 |
Description: | Value stored to 'rad' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2004, The GROMACS development team. |
6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
8 | * and including many others, as listed in the AUTHORS file in the |
9 | * top-level source directory and at http://www.gromacs.org. |
10 | * |
11 | * GROMACS is free software; you can redistribute it and/or |
12 | * modify it under the terms of the GNU Lesser General Public License |
13 | * as published by the Free Software Foundation; either version 2.1 |
14 | * of the License, or (at your option) any later version. |
15 | * |
16 | * GROMACS is distributed in the hope that it will be useful, |
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
19 | * Lesser General Public License for more details. |
20 | * |
21 | * You should have received a copy of the GNU Lesser General Public |
22 | * License along with GROMACS; if not, see |
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
25 | * |
26 | * If you want to redistribute modifications to GROMACS, please |
27 | * consider that scientific software is very special. Version |
28 | * control is crucial - bugs must be traceable. We will be happy to |
29 | * consider code for inclusion in the official distribution, but |
30 | * derived work must not be called official GROMACS. Details are found |
31 | * in the README & COPYING files - if they are missing, get the |
32 | * official version at http://www.gromacs.org. |
33 | * |
34 | * To help us fund GROMACS development, we humbly ask that you cite |
35 | * the research papers on the package. Check out http://www.gromacs.org. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <stdio.h> |
42 | #include <string.h> |
43 | #include <time.h> |
44 | |
45 | #include "typedefs.h" |
46 | #include "gromacs/utility/cstringutil.h" |
47 | #include "gromacs/utility/smalloc.h" |
48 | #include "names.h" |
49 | #include "gromacs/fileio/confio.h" |
50 | #include "txtdump.h" |
51 | #include "gromacs/math/vec.h" |
52 | #include "nrnb.h" |
53 | #include "mshift.h" |
54 | #include "mdrun.h" |
55 | #include "update.h" |
56 | #include "physics.h" |
57 | #include "mtop_util.h" |
58 | #include "gromacs/essentialdynamics/edsam.h" |
59 | #include "gromacs/fileio/gmxfio.h" |
60 | #include "gromacs/fileio/xvgr.h" |
61 | #include "gromacs/mdlib/groupcoord.h" |
62 | |
63 | #include "gromacs/linearalgebra/nrjac.h" |
64 | #include "gromacs/utility/fatalerror.h" |
65 | |
66 | /* We use the same defines as in mvdata.c here */ |
67 | #define block_bc(cr, d)gmx_bcast( sizeof(d), &(d), (cr)) gmx_bcast( sizeof(d), &(d), (cr)) |
68 | #define nblock_bc(cr, nr, d)gmx_bcast((nr)*sizeof((d)[0]), (d), (cr)) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr)) |
69 | #define snew_bc(cr, d, nr){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((d)) = save_calloc("(d)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 69, ((nr)), sizeof(*((d)))); }} { if (!MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {snew((d), (nr))((d)) = save_calloc("(d)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 69, ((nr)), sizeof(*((d)))); }} |
70 | |
71 | /* These macros determine the column width in the output file */ |
72 | #define EDcol_sfmt"%17s" "%17s" |
73 | #define EDcol_efmt"%17.5e" "%17.5e" |
74 | #define EDcol_ffmt"%17f" "%17f" |
75 | |
76 | /* enum to identify the type of ED: none, normal ED, flooding */ |
77 | enum { |
78 | eEDnone, eEDedsam, eEDflood, eEDnr |
79 | }; |
80 | |
81 | /* enum to identify operations on reference, average, origin, target structures */ |
82 | enum { |
83 | eedREF, eedAV, eedORI, eedTAR, eedNR |
84 | }; |
85 | |
86 | |
87 | typedef struct |
88 | { |
89 | int neig; /* nr of eigenvectors */ |
90 | int *ieig; /* index nrs of eigenvectors */ |
91 | real *stpsz; /* stepsizes (per eigenvector) */ |
92 | rvec **vec; /* eigenvector components */ |
93 | real *xproj; /* instantaneous x projections */ |
94 | real *fproj; /* instantaneous f projections */ |
95 | real radius; /* instantaneous radius */ |
96 | real *refproj; /* starting or target projections */ |
97 | /* When using flooding as harmonic restraint: The current reference projection |
98 | * is at each step calculated from the initial refproj0 and the slope. */ |
99 | real *refproj0, *refprojslope; |
100 | } t_eigvec; |
101 | |
102 | |
103 | typedef struct |
104 | { |
105 | t_eigvec mon; /* only monitored, no constraints */ |
106 | t_eigvec linfix; /* fixed linear constraints */ |
107 | t_eigvec linacc; /* acceptance linear constraints */ |
108 | t_eigvec radfix; /* fixed radial constraints (exp) */ |
109 | t_eigvec radacc; /* acceptance radial constraints (exp) */ |
110 | t_eigvec radcon; /* acceptance rad. contraction constr. */ |
111 | } t_edvecs; |
112 | |
113 | |
114 | typedef struct |
115 | { |
116 | real deltaF0; |
117 | gmx_bool bHarmonic; /* Use flooding for harmonic restraint on |
118 | the eigenvector */ |
119 | gmx_bool bConstForce; /* Do not calculate a flooding potential, |
120 | instead flood with a constant force */ |
121 | real tau; |
122 | real deltaF; |
123 | real Efl; |
124 | real kT; |
125 | real Vfl; |
126 | real dt; |
127 | real constEfl; |
128 | real alpha2; |
129 | rvec *forces_cartesian; |
130 | t_eigvec vecs; /* use flooding for these */ |
131 | } t_edflood; |
132 | |
133 | |
134 | /* This type is for the average, reference, target, and origin structure */ |
135 | typedef struct gmx_edx |
136 | { |
137 | int nr; /* number of atoms this structure contains */ |
138 | int nr_loc; /* number of atoms on local node */ |
139 | int *anrs; /* atom index numbers */ |
140 | int *anrs_loc; /* local atom index numbers */ |
141 | int nalloc_loc; /* allocation size of anrs_loc */ |
142 | int *c_ind; /* at which position of the whole anrs |
143 | * array is a local atom?, i.e. |
144 | * c_ind[0...nr_loc-1] gives the atom index |
145 | * with respect to the collective |
146 | * anrs[0...nr-1] array */ |
147 | rvec *x; /* positions for this structure */ |
148 | rvec *x_old; /* Last positions which have the correct PBC |
149 | representation of the ED group. In |
150 | combination with keeping track of the |
151 | shift vectors, the ED group can always |
152 | be made whole */ |
153 | real *m; /* masses */ |
154 | real mtot; /* total mass (only used in sref) */ |
155 | real *sqrtm; /* sqrt of the masses used for mass- |
156 | * weighting of analysis (only used in sav) */ |
157 | } t_gmx_edx; |
158 | |
159 | |
160 | typedef struct edpar |
161 | { |
162 | int nini; /* total Nr of atoms */ |
163 | gmx_bool fitmas; /* true if trans fit with cm */ |
164 | gmx_bool pcamas; /* true if mass-weighted PCA */ |
165 | int presteps; /* number of steps to run without any |
166 | * perturbations ... just monitoring */ |
167 | int outfrq; /* freq (in steps) of writing to edo */ |
168 | int maxedsteps; /* max nr of steps per cycle */ |
169 | |
170 | /* all gmx_edx datasets are copied to all nodes in the parallel case */ |
171 | struct gmx_edx sref; /* reference positions, to these fitting |
172 | * will be done */ |
173 | gmx_bool bRefEqAv; /* If true, reference & average indices |
174 | * are the same. Used for optimization */ |
175 | struct gmx_edx sav; /* average positions */ |
176 | struct gmx_edx star; /* target positions */ |
177 | struct gmx_edx sori; /* origin positions */ |
178 | |
179 | t_edvecs vecs; /* eigenvectors */ |
180 | real slope; /* minimal slope in acceptance radexp */ |
181 | |
182 | t_edflood flood; /* parameters especially for flooding */ |
183 | struct t_ed_buffer *buf; /* handle to local buffers */ |
184 | struct edpar *next_edi; /* Pointer to another ED group */ |
185 | } t_edpar; |
186 | |
187 | |
188 | typedef struct gmx_edsam |
189 | { |
190 | int eEDtype; /* Type of ED: see enums above */ |
191 | FILE *edo; /* output file pointer */ |
192 | t_edpar *edpar; |
193 | gmx_bool bFirst; |
194 | } t_gmx_edsam; |
195 | |
196 | |
197 | struct t_do_edsam |
198 | { |
199 | matrix old_rotmat; |
200 | real oldrad; |
201 | rvec old_transvec, older_transvec, transvec_compact; |
202 | rvec *xcoll; /* Positions from all nodes, this is the |
203 | collective set we work on. |
204 | These are the positions of atoms with |
205 | average structure indices */ |
206 | rvec *xc_ref; /* same but with reference structure indices */ |
207 | ivec *shifts_xcoll; /* Shifts for xcoll */ |
208 | ivec *extra_shifts_xcoll; /* xcoll shift changes since last NS step */ |
209 | ivec *shifts_xc_ref; /* Shifts for xc_ref */ |
210 | ivec *extra_shifts_xc_ref; /* xc_ref shift changes since last NS step */ |
211 | gmx_bool bUpdateShifts; /* TRUE in NS steps to indicate that the |
212 | ED shifts for this ED group need to |
213 | be updated */ |
214 | }; |
215 | |
216 | |
217 | /* definition of ED buffer structure */ |
218 | struct t_ed_buffer |
219 | { |
220 | struct t_fit_to_ref * fit_to_ref; |
221 | struct t_do_edfit * do_edfit; |
222 | struct t_do_edsam * do_edsam; |
223 | struct t_do_radcon * do_radcon; |
224 | }; |
225 | |
226 | |
227 | /* Function declarations */ |
228 | static void fit_to_reference(rvec *xcoll, rvec transvec, matrix rotmat, t_edpar *edi); |
229 | static void translate_and_rotate(rvec *x, int nat, rvec transvec, matrix rotmat); |
230 | static real rmsd_from_structure(rvec *x, struct gmx_edx *s); |
231 | static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms); |
232 | static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate); |
233 | static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate); |
234 | static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv); |
235 | /* End function declarations */ |
236 | |
237 | |
238 | /* Do we have to perform essential dynamics constraints or possibly only flooding |
239 | * for any of the ED groups? */ |
240 | static gmx_bool bNeedDoEdsam(t_edpar *edi) |
241 | { |
242 | return edi->vecs.mon.neig |
243 | || edi->vecs.linfix.neig |
244 | || edi->vecs.linacc.neig |
245 | || edi->vecs.radfix.neig |
246 | || edi->vecs.radacc.neig |
247 | || edi->vecs.radcon.neig; |
248 | } |
249 | |
250 | |
251 | /* Multiple ED groups will be labeled with letters instead of numbers |
252 | * to avoid confusion with eigenvector indices */ |
253 | static char get_EDgroupChar(int nr_edi, int nED) |
254 | { |
255 | if (nED == 1) |
256 | { |
257 | return ' '; |
258 | } |
259 | |
260 | /* nr_edi = 1 -> A |
261 | * nr_edi = 2 -> B ... |
262 | */ |
263 | return 'A' + nr_edi - 1; |
264 | } |
265 | |
266 | |
267 | /* Does not subtract average positions, projection on single eigenvector is returned |
268 | * used by: do_linfix, do_linacc, do_radfix, do_radacc, do_radcon |
269 | * Average position is subtracted in ed_apply_constraints prior to calling projectx |
270 | */ |
271 | static real projectx(t_edpar *edi, rvec *xcoll, rvec *vec) |
272 | { |
273 | int i; |
274 | real proj = 0.0; |
275 | |
276 | |
277 | for (i = 0; i < edi->sav.nr; i++) |
278 | { |
279 | proj += edi->sav.sqrtm[i]*iprod(vec[i], xcoll[i]); |
280 | } |
281 | |
282 | return proj; |
283 | } |
284 | |
285 | |
286 | /* Specialized: projection is stored in vec->refproj |
287 | * -> used for radacc, radfix, radcon and center of flooding potential |
288 | * subtracts average positions, projects vector x */ |
289 | static void rad_project(t_edpar *edi, rvec *x, t_eigvec *vec) |
290 | { |
291 | int i; |
292 | real rad = 0.0; |
293 | |
294 | /* Subtract average positions */ |
295 | for (i = 0; i < edi->sav.nr; i++) |
296 | { |
297 | rvec_dec(x[i], edi->sav.x[i]); |
298 | } |
299 | |
300 | for (i = 0; i < vec->neig; i++) |
301 | { |
302 | vec->refproj[i] = projectx(edi, x, vec->vec[i]); |
303 | rad += pow((vec->refproj[i]-vec->xproj[i]), 2); |
304 | } |
305 | vec->radius = sqrt(rad); |
306 | |
307 | /* Add average positions */ |
308 | for (i = 0; i < edi->sav.nr; i++) |
309 | { |
310 | rvec_inc(x[i], edi->sav.x[i]); |
311 | } |
312 | } |
313 | |
314 | |
315 | /* Project vector x, subtract average positions prior to projection and add |
316 | * them afterwards to retain the unchanged vector. Store in xproj. Mass-weighting |
317 | * is applied. */ |
318 | static void project_to_eigvectors(rvec *x, /* The positions to project to an eigenvector */ |
319 | t_eigvec *vec, /* The eigenvectors */ |
320 | t_edpar *edi) |
321 | { |
322 | int i; |
323 | |
324 | |
325 | if (!vec->neig) |
326 | { |
327 | return; |
328 | } |
329 | |
330 | /* Subtract average positions */ |
331 | for (i = 0; i < edi->sav.nr; i++) |
332 | { |
333 | rvec_dec(x[i], edi->sav.x[i]); |
334 | } |
335 | |
336 | for (i = 0; i < vec->neig; i++) |
337 | { |
338 | vec->xproj[i] = projectx(edi, x, vec->vec[i]); |
339 | } |
340 | |
341 | /* Add average positions */ |
342 | for (i = 0; i < edi->sav.nr; i++) |
343 | { |
344 | rvec_inc(x[i], edi->sav.x[i]); |
345 | } |
346 | } |
347 | |
348 | |
349 | /* Project vector x onto all edi->vecs (mon, linfix,...) */ |
350 | static void project(rvec *x, /* positions to project */ |
351 | t_edpar *edi) /* edi data set */ |
352 | { |
353 | /* It is not more work to subtract the average position in every |
354 | * subroutine again, because these routines are rarely used simultaneously */ |
355 | project_to_eigvectors(x, &edi->vecs.mon, edi); |
356 | project_to_eigvectors(x, &edi->vecs.linfix, edi); |
357 | project_to_eigvectors(x, &edi->vecs.linacc, edi); |
358 | project_to_eigvectors(x, &edi->vecs.radfix, edi); |
359 | project_to_eigvectors(x, &edi->vecs.radacc, edi); |
360 | project_to_eigvectors(x, &edi->vecs.radcon, edi); |
361 | } |
362 | |
363 | |
364 | static real calc_radius(t_eigvec *vec) |
365 | { |
366 | int i; |
367 | real rad = 0.0; |
368 | |
369 | |
370 | for (i = 0; i < vec->neig; i++) |
371 | { |
372 | rad += pow((vec->refproj[i]-vec->xproj[i]), 2); |
373 | } |
374 | |
375 | return rad = sqrt(rad); |
376 | } |
377 | |
378 | |
379 | /* Debug helper */ |
380 | #ifdef DEBUGHELPERS |
381 | static void dump_xcoll(t_edpar *edi, struct t_do_edsam *buf, t_commrec *cr, |
382 | int step) |
383 | { |
384 | int i; |
385 | FILE *fp; |
386 | char fn[STRLEN4096]; |
387 | rvec *xcoll; |
388 | ivec *shifts, *eshifts; |
389 | |
390 | |
391 | if (!MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
392 | { |
393 | return; |
394 | } |
395 | |
396 | xcoll = buf->xcoll; |
397 | shifts = buf->shifts_xcoll; |
398 | eshifts = buf->extra_shifts_xcoll; |
399 | |
400 | sprintf(fn, "xcolldump_step%d.txt", step); |
401 | fp = fopen(fn, "w"); |
402 | |
403 | for (i = 0; i < edi->sav.nr; i++) |
404 | { |
405 | fprintf(fp, "%d %9.5f %9.5f %9.5f %d %d %d %d %d %d\n", |
406 | edi->sav.anrs[i]+1, |
407 | xcoll[i][XX0], xcoll[i][YY1], xcoll[i][ZZ2], |
408 | shifts[i][XX0], shifts[i][YY1], shifts[i][ZZ2], |
409 | eshifts[i][XX0], eshifts[i][YY1], eshifts[i][ZZ2]); |
410 | } |
411 | |
412 | fclose(fp); |
413 | } |
414 | |
415 | |
416 | /* Debug helper */ |
417 | static void dump_edi_positions(FILE *out, struct gmx_edx *s, const char name[]) |
418 | { |
419 | int i; |
420 | |
421 | |
422 | fprintf(out, "#%s positions:\n%d\n", name, s->nr); |
423 | if (s->nr == 0) |
424 | { |
425 | return; |
426 | } |
427 | |
428 | fprintf(out, "#index, x, y, z"); |
429 | if (s->sqrtm) |
430 | { |
431 | fprintf(out, ", sqrt(m)"); |
432 | } |
433 | for (i = 0; i < s->nr; i++) |
434 | { |
435 | fprintf(out, "\n%6d %11.6f %11.6f %11.6f", s->anrs[i], s->x[i][XX0], s->x[i][YY1], s->x[i][ZZ2]); |
436 | if (s->sqrtm) |
437 | { |
438 | fprintf(out, "%9.3f", s->sqrtm[i]); |
439 | } |
440 | } |
441 | fprintf(out, "\n"); |
442 | } |
443 | |
444 | |
445 | /* Debug helper */ |
446 | static void dump_edi_eigenvecs(FILE *out, t_eigvec *ev, |
447 | const char name[], int length) |
448 | { |
449 | int i, j; |
450 | |
451 | |
452 | fprintf(out, "#%s eigenvectors:\n%d\n", name, ev->neig); |
453 | /* Dump the data for every eigenvector: */ |
454 | for (i = 0; i < ev->neig; i++) |
455 | { |
456 | fprintf(out, "EV %4d\ncomponents %d\nstepsize %f\nxproj %f\nfproj %f\nrefproj %f\nradius %f\nComponents:\n", |
457 | ev->ieig[i], length, ev->stpsz[i], ev->xproj[i], ev->fproj[i], ev->refproj[i], ev->radius); |
458 | for (j = 0; j < length; j++) |
459 | { |
460 | fprintf(out, "%11.6f %11.6f %11.6f\n", ev->vec[i][j][XX0], ev->vec[i][j][YY1], ev->vec[i][j][ZZ2]); |
461 | } |
462 | } |
463 | } |
464 | |
465 | |
466 | /* Debug helper */ |
467 | static void dump_edi(t_edpar *edpars, t_commrec *cr, int nr_edi) |
468 | { |
469 | FILE *out; |
470 | char fn[STRLEN4096]; |
471 | |
472 | |
473 | sprintf(fn, "EDdump_node%d_edi%d", cr->nodeid, nr_edi); |
474 | out = gmx_ffopen(fn, "w"); |
475 | |
476 | fprintf(out, "#NINI\n %d\n#FITMAS\n %d\n#ANALYSIS_MAS\n %d\n", |
477 | edpars->nini, edpars->fitmas, edpars->pcamas); |
478 | fprintf(out, "#OUTFRQ\n %d\n#MAXLEN\n %d\n#SLOPECRIT\n %f\n", |
479 | edpars->outfrq, edpars->maxedsteps, edpars->slope); |
480 | fprintf(out, "#PRESTEPS\n %d\n#DELTA_F0\n %f\n#TAU\n %f\n#EFL_NULL\n %f\n#ALPHA2\n %f\n", |
481 | edpars->presteps, edpars->flood.deltaF0, edpars->flood.tau, |
482 | edpars->flood.constEfl, edpars->flood.alpha2); |
483 | |
484 | /* Dump reference, average, target, origin positions */ |
485 | dump_edi_positions(out, &edpars->sref, "REFERENCE"); |
486 | dump_edi_positions(out, &edpars->sav, "AVERAGE" ); |
487 | dump_edi_positions(out, &edpars->star, "TARGET" ); |
488 | dump_edi_positions(out, &edpars->sori, "ORIGIN" ); |
489 | |
490 | /* Dump eigenvectors */ |
491 | dump_edi_eigenvecs(out, &edpars->vecs.mon, "MONITORED", edpars->sav.nr); |
492 | dump_edi_eigenvecs(out, &edpars->vecs.linfix, "LINFIX", edpars->sav.nr); |
493 | dump_edi_eigenvecs(out, &edpars->vecs.linacc, "LINACC", edpars->sav.nr); |
494 | dump_edi_eigenvecs(out, &edpars->vecs.radfix, "RADFIX", edpars->sav.nr); |
495 | dump_edi_eigenvecs(out, &edpars->vecs.radacc, "RADACC", edpars->sav.nr); |
496 | dump_edi_eigenvecs(out, &edpars->vecs.radcon, "RADCON", edpars->sav.nr); |
497 | |
498 | /* Dump flooding eigenvectors */ |
499 | dump_edi_eigenvecs(out, &edpars->flood.vecs, "FLOODING", edpars->sav.nr); |
500 | |
501 | /* Dump ed local buffer */ |
502 | fprintf(out, "buf->do_edfit =%p\n", (void*)edpars->buf->do_edfit ); |
503 | fprintf(out, "buf->do_edsam =%p\n", (void*)edpars->buf->do_edsam ); |
504 | fprintf(out, "buf->do_radcon =%p\n", (void*)edpars->buf->do_radcon ); |
505 | |
506 | gmx_ffclose(out); |
507 | } |
508 | |
509 | |
510 | /* Debug helper */ |
511 | static void dump_rotmat(FILE* out, matrix rotmat) |
512 | { |
513 | fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[XX0][XX0], rotmat[XX0][YY1], rotmat[XX0][ZZ2]); |
514 | fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[YY1][XX0], rotmat[YY1][YY1], rotmat[YY1][ZZ2]); |
515 | fprintf(out, "ROTMAT: %12.8f %12.8f %12.8f\n", rotmat[ZZ2][XX0], rotmat[ZZ2][YY1], rotmat[ZZ2][ZZ2]); |
516 | } |
517 | |
518 | |
519 | /* Debug helper */ |
520 | static void dump_rvec(FILE *out, int dim, rvec *x) |
521 | { |
522 | int i; |
523 | |
524 | |
525 | for (i = 0; i < dim; i++) |
526 | { |
527 | fprintf(out, "%4d %f %f %f\n", i, x[i][XX0], x[i][YY1], x[i][ZZ2]); |
528 | } |
529 | } |
530 | |
531 | |
532 | /* Debug helper */ |
533 | static void dump_mat(FILE* out, int dim, double** mat) |
534 | { |
535 | int i, j; |
536 | |
537 | |
538 | fprintf(out, "MATRIX:\n"); |
539 | for (i = 0; i < dim; i++) |
540 | { |
541 | for (j = 0; j < dim; j++) |
542 | { |
543 | fprintf(out, "%f ", mat[i][j]); |
544 | } |
545 | fprintf(out, "\n"); |
546 | } |
547 | } |
548 | #endif |
549 | |
550 | |
551 | struct t_do_edfit { |
552 | double **omega; |
553 | double **om; |
554 | }; |
555 | |
556 | static void do_edfit(int natoms, rvec *xp, rvec *x, matrix R, t_edpar *edi) |
557 | { |
558 | /* this is a copy of do_fit with some modifications */ |
559 | int c, r, n, j, i, irot; |
560 | double d[6], xnr, xpc; |
561 | matrix vh, vk, u; |
562 | int index; |
563 | real max_d; |
564 | |
565 | struct t_do_edfit *loc; |
566 | gmx_bool bFirst; |
567 | |
568 | if (edi->buf->do_edfit != NULL((void*)0)) |
569 | { |
570 | bFirst = FALSE0; |
571 | } |
572 | else |
573 | { |
574 | bFirst = TRUE1; |
575 | snew(edi->buf->do_edfit, 1)(edi->buf->do_edfit) = save_calloc("edi->buf->do_edfit" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 575, (1), sizeof(*(edi->buf->do_edfit))); |
576 | } |
577 | loc = edi->buf->do_edfit; |
578 | |
579 | if (bFirst) |
580 | { |
581 | snew(loc->omega, 2*DIM)(loc->omega) = save_calloc("loc->omega", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 581, (2*3), sizeof(*(loc->omega))); |
582 | snew(loc->om, 2*DIM)(loc->om) = save_calloc("loc->om", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 582, (2*3), sizeof(*(loc->om))); |
583 | for (i = 0; i < 2*DIM3; i++) |
584 | { |
585 | snew(loc->omega[i], 2*DIM)(loc->omega[i]) = save_calloc("loc->omega[i]", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 585, (2*3), sizeof(*(loc->omega[i]))); |
586 | snew(loc->om[i], 2*DIM)(loc->om[i]) = save_calloc("loc->om[i]", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 586, (2*3), sizeof(*(loc->om[i]))); |
587 | } |
588 | } |
589 | |
590 | for (i = 0; (i < 6); i++) |
591 | { |
592 | d[i] = 0; |
593 | for (j = 0; (j < 6); j++) |
594 | { |
595 | loc->omega[i][j] = 0; |
596 | loc->om[i][j] = 0; |
597 | } |
598 | } |
599 | |
600 | /* calculate the matrix U */ |
601 | clear_mat(u); |
602 | for (n = 0; (n < natoms); n++) |
603 | { |
604 | for (c = 0; (c < DIM3); c++) |
605 | { |
606 | xpc = xp[n][c]; |
607 | for (r = 0; (r < DIM3); r++) |
608 | { |
609 | xnr = x[n][r]; |
610 | u[c][r] += xnr*xpc; |
611 | } |
612 | } |
613 | } |
614 | |
615 | /* construct loc->omega */ |
616 | /* loc->omega is symmetric -> loc->omega==loc->omega' */ |
617 | for (r = 0; (r < 6); r++) |
618 | { |
619 | for (c = 0; (c <= r); c++) |
620 | { |
621 | if ((r >= 3) && (c < 3)) |
622 | { |
623 | loc->omega[r][c] = u[r-3][c]; |
624 | loc->omega[c][r] = u[r-3][c]; |
625 | } |
626 | else |
627 | { |
628 | loc->omega[r][c] = 0; |
629 | loc->omega[c][r] = 0; |
630 | } |
631 | } |
632 | } |
633 | |
634 | /* determine h and k */ |
635 | #ifdef DEBUG |
636 | { |
637 | int i; |
638 | dump_mat(stderrstderr, 2*DIM3, loc->omega); |
639 | for (i = 0; i < 6; i++) |
640 | { |
641 | fprintf(stderrstderr, "d[%d] = %f\n", i, d[i]); |
642 | } |
643 | } |
644 | #endif |
645 | jacobi(loc->omega, 6, d, loc->om, &irot); |
646 | |
647 | if (irot == 0) |
648 | { |
649 | fprintf(stderrstderr, "IROT=0\n"); |
650 | } |
651 | |
652 | index = 0; /* For the compiler only */ |
653 | |
654 | for (j = 0; (j < 3); j++) |
655 | { |
656 | max_d = -1000; |
657 | for (i = 0; (i < 6); i++) |
658 | { |
659 | if (d[i] > max_d) |
660 | { |
661 | max_d = d[i]; |
662 | index = i; |
663 | } |
664 | } |
665 | d[index] = -10000; |
666 | for (i = 0; (i < 3); i++) |
667 | { |
668 | vh[j][i] = M_SQRT21.41421356237309504880*loc->om[i][index]; |
669 | vk[j][i] = M_SQRT21.41421356237309504880*loc->om[i+DIM3][index]; |
670 | } |
671 | } |
672 | |
673 | /* determine R */ |
674 | for (c = 0; (c < 3); c++) |
675 | { |
676 | for (r = 0; (r < 3); r++) |
677 | { |
678 | R[c][r] = vk[0][r]*vh[0][c]+ |
679 | vk[1][r]*vh[1][c]+ |
680 | vk[2][r]*vh[2][c]; |
681 | } |
682 | } |
683 | if (det(R) < 0) |
684 | { |
685 | for (c = 0; (c < 3); c++) |
686 | { |
687 | for (r = 0; (r < 3); r++) |
688 | { |
689 | R[c][r] = vk[0][r]*vh[0][c]+ |
690 | vk[1][r]*vh[1][c]- |
691 | vk[2][r]*vh[2][c]; |
692 | } |
693 | } |
694 | } |
695 | } |
696 | |
697 | |
698 | static void rmfit(int nat, rvec *xcoll, rvec transvec, matrix rotmat) |
699 | { |
700 | rvec vec; |
701 | matrix tmat; |
702 | |
703 | |
704 | /* Remove rotation. |
705 | * The inverse rotation is described by the transposed rotation matrix */ |
706 | transpose(rotmat, tmat); |
707 | rotate_x(xcoll, nat, tmat); |
708 | |
709 | /* Remove translation */ |
710 | vec[XX0] = -transvec[XX0]; |
711 | vec[YY1] = -transvec[YY1]; |
712 | vec[ZZ2] = -transvec[ZZ2]; |
713 | translate_x(xcoll, nat, vec); |
714 | } |
715 | |
716 | |
717 | /********************************************************************************** |
718 | ******************** FLOODING **************************************************** |
719 | ********************************************************************************** |
720 | |
721 | The flooding ability was added later to edsam. Many of the edsam functionality could be reused for that purpose. |
722 | The flooding covariance matrix, i.e. the selected eigenvectors and their corresponding eigenvalues are |
723 | read as 7th Component Group. The eigenvalues are coded into the stepsize parameter (as used by -linfix or -linacc). |
724 | |
725 | do_md clls right in the beginning the function init_edsam, which reads the edi file, saves all the necessary information in |
726 | the edi structure and calls init_flood, to initialise some extra fields in the edi->flood structure. |
727 | |
728 | since the flooding acts on forces do_flood is called from the function force() (force.c), while the other |
729 | edsam functionality is hooked into md via the update() (update.c) function acting as constraint on positions. |
730 | |
731 | do_flood makes a copy of the positions, |
732 | fits them, projects them computes flooding_energy, and flooding forces. The forces are computed in the |
733 | space of the eigenvectors and are then blown up to the full cartesian space and rotated back to remove the |
734 | fit. Then do_flood adds these forces to the forcefield-forces |
735 | (given as parameter) and updates the adaptive flooding parameters Efl and deltaF. |
736 | |
737 | To center the flooding potential at a different location one can use the -ori option in make_edi. The ori |
738 | structure is projected to the system of eigenvectors and then this position in the subspace is used as |
739 | center of the flooding potential. If the option is not used, the center will be zero in the subspace, |
740 | i.e. the average structure as given in the make_edi file. |
741 | |
742 | To use the flooding potential as restraint, make_edi has the option -restrain, which leads to inverted |
743 | signs of alpha2 and Efl, such that the sign in the exponential of Vfl is not inverted but the sign of |
744 | Vfl is inverted. Vfl = Efl * exp (- .../Efl/alpha2*x^2...) With tau>0 the negative Efl will grow slowly |
745 | so that the restraint is switched off slowly. When Efl==0 and inverted flooding is ON is reached no |
746 | further adaption is applied, Efl will stay constant at zero. |
747 | |
748 | To use restraints with harmonic potentials switch -restrain and -harmonic. Then the eigenvalues are |
749 | used as spring constants for the harmonic potential. |
750 | Note that eq3 in the flooding paper (J. Comp. Chem. 2006, 27, 1693-1702) defines the parameter lambda \ |
751 | as the inverse of the spring constant, whereas the implementation uses lambda as the spring constant. |
752 | |
753 | To use more than one flooding matrix just concatenate several .edi files (cat flood1.edi flood2.edi > flood_all.edi) |
754 | the routine read_edi_file reads all of theses flooding files. |
755 | The structure t_edi is now organized as a list of t_edis and the function do_flood cycles through the list |
756 | calling the do_single_flood() routine for every single entry. Since every state variables have been kept in one |
757 | edi there is no interdependence whatsoever. The forces are added together. |
758 | |
759 | To write energies into the .edr file, call the function |
760 | get_flood_enx_names(char**, int *nnames) to get the Header (Vfl1 Vfl2... Vfln) |
761 | and call |
762 | get_flood_energies(real Vfl[],int nnames); |
763 | |
764 | TODO: |
765 | - one could program the whole thing such that Efl, Vfl and deltaF is written to the .edr file. -- i dont know how to do that, yet. |
766 | |
767 | Maybe one should give a range of atoms for which to remove motion, so that motion is removed with |
768 | two edsam files from two peptide chains |
769 | */ |
770 | |
771 | static void write_edo_flood(t_edpar *edi, FILE *fp, real rmsd) |
772 | { |
773 | int i; |
774 | |
775 | |
776 | /* Output how well we fit to the reference structure */ |
777 | fprintf(fp, EDcol_ffmt"%17f", rmsd); |
778 | |
779 | for (i = 0; i < edi->flood.vecs.neig; i++) |
780 | { |
781 | fprintf(fp, EDcol_efmt"%17.5e", edi->flood.vecs.xproj[i]); |
782 | |
783 | /* Check whether the reference projection changes with time (this can happen |
784 | * in case flooding is used as harmonic restraint). If so, output the |
785 | * current reference projection */ |
786 | if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0) |
787 | { |
788 | fprintf(fp, EDcol_efmt"%17.5e", edi->flood.vecs.refproj[i]); |
789 | } |
790 | |
791 | /* Output Efl if we are doing adaptive flooding */ |
792 | if (0 != edi->flood.tau) |
793 | { |
794 | fprintf(fp, EDcol_efmt"%17.5e", edi->flood.Efl); |
795 | } |
796 | fprintf(fp, EDcol_efmt"%17.5e", edi->flood.Vfl); |
797 | |
798 | /* Output deltaF if we are doing adaptive flooding */ |
799 | if (0 != edi->flood.tau) |
800 | { |
801 | fprintf(fp, EDcol_efmt"%17.5e", edi->flood.deltaF); |
802 | } |
803 | fprintf(fp, EDcol_efmt"%17.5e", edi->flood.vecs.fproj[i]); |
804 | } |
805 | } |
806 | |
807 | |
808 | /* From flood.xproj compute the Vfl(x) at this point */ |
809 | static real flood_energy(t_edpar *edi, gmx_int64_t step) |
810 | { |
811 | /* compute flooding energy Vfl |
812 | Vfl = Efl * exp( - \frac {kT} {2Efl alpha^2} * sum_i { \lambda_i c_i^2 } ) |
813 | \lambda_i is the reciprocal eigenvalue 1/\sigma_i |
814 | it is already computed by make_edi and stored in stpsz[i] |
815 | bHarmonic: |
816 | Vfl = - Efl * 1/2(sum _i {\frac 1{\lambda_i} c_i^2}) |
817 | */ |
818 | real sum; |
819 | real Vfl; |
820 | int i; |
821 | |
822 | |
823 | /* Each time this routine is called (i.e. each time step), we add a small |
824 | * value to the reference projection. This way a harmonic restraint towards |
825 | * a moving reference is realized. If no value for the additive constant |
826 | * is provided in the edi file, the reference will not change. */ |
827 | if (edi->flood.bHarmonic) |
828 | { |
829 | for (i = 0; i < edi->flood.vecs.neig; i++) |
830 | { |
831 | edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i] + step * edi->flood.vecs.refprojslope[i]; |
832 | } |
833 | } |
834 | |
835 | sum = 0.0; |
836 | /* Compute sum which will be the exponent of the exponential */ |
837 | for (i = 0; i < edi->flood.vecs.neig; i++) |
838 | { |
839 | /* stpsz stores the reciprocal eigenvalue 1/sigma_i */ |
840 | sum += edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i])*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]); |
841 | } |
842 | |
843 | /* Compute the Gauss function*/ |
844 | if (edi->flood.bHarmonic) |
845 | { |
846 | Vfl = -0.5*edi->flood.Efl*sum; /* minus sign because Efl is negative, if restrain is on. */ |
847 | } |
848 | else |
849 | { |
850 | Vfl = edi->flood.Efl != 0 ? edi->flood.Efl*exp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2*sum) : 0; |
851 | } |
852 | |
853 | return Vfl; |
854 | } |
855 | |
856 | |
857 | /* From the position and from Vfl compute forces in subspace -> store in edi->vec.flood.fproj */ |
858 | static void flood_forces(t_edpar *edi) |
859 | { |
860 | /* compute the forces in the subspace of the flooding eigenvectors |
861 | * by the formula F_i= V_{fl}(c) * ( \frac {kT} {E_{fl}} \lambda_i c_i */ |
862 | |
863 | int i; |
864 | real energy = edi->flood.Vfl; |
865 | |
866 | |
867 | if (edi->flood.bHarmonic) |
868 | { |
869 | for (i = 0; i < edi->flood.vecs.neig; i++) |
870 | { |
871 | edi->flood.vecs.fproj[i] = edi->flood.Efl* edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]); |
872 | } |
873 | } |
874 | else |
875 | { |
876 | for (i = 0; i < edi->flood.vecs.neig; i++) |
877 | { |
878 | /* if Efl is zero the forces are zero if not use the formula */ |
879 | edi->flood.vecs.fproj[i] = edi->flood.Efl != 0 ? edi->flood.kT/edi->flood.Efl/edi->flood.alpha2*energy*edi->flood.vecs.stpsz[i]*(edi->flood.vecs.xproj[i]-edi->flood.vecs.refproj[i]) : 0; |
880 | } |
881 | } |
882 | } |
883 | |
884 | |
885 | /* Raise forces from subspace into cartesian space */ |
886 | static void flood_blowup(t_edpar *edi, rvec *forces_cart) |
887 | { |
888 | /* this function lifts the forces from the subspace to the cartesian space |
889 | all the values not contained in the subspace are assumed to be zero and then |
890 | a coordinate transformation from eigenvector to cartesian vectors is performed |
891 | The nonexistent values don't have to be set to zero explicitly, they would occur |
892 | as zero valued summands, hence we just stop to compute this part of the sum. |
893 | |
894 | for every atom we add all the contributions to this atom from all the different eigenvectors. |
895 | |
896 | NOTE: one could add directly to the forcefield forces, would mean we wouldn't have to clear the |
897 | field forces_cart prior the computation, but we compute the forces separately |
898 | to have them accessible for diagnostics |
899 | */ |
900 | int j, eig; |
901 | rvec dum; |
902 | real *forces_sub; |
903 | |
904 | |
905 | forces_sub = edi->flood.vecs.fproj; |
906 | |
907 | |
908 | /* Calculate the cartesian forces for the local atoms */ |
909 | |
910 | /* Clear forces first */ |
911 | for (j = 0; j < edi->sav.nr_loc; j++) |
912 | { |
913 | clear_rvec(forces_cart[j]); |
914 | } |
915 | |
916 | /* Now compute atomwise */ |
917 | for (j = 0; j < edi->sav.nr_loc; j++) |
918 | { |
919 | /* Compute forces_cart[edi->sav.anrs[j]] */ |
920 | for (eig = 0; eig < edi->flood.vecs.neig; eig++) |
921 | { |
922 | /* Force vector is force * eigenvector (compute only atom j) */ |
923 | svmul(forces_sub[eig], edi->flood.vecs.vec[eig][edi->sav.c_ind[j]], dum); |
924 | /* Add this vector to the cartesian forces */ |
925 | rvec_inc(forces_cart[j], dum); |
926 | } |
927 | } |
928 | } |
929 | |
930 | |
931 | /* Update the values of Efl, deltaF depending on tau and Vfl */ |
932 | static void update_adaption(t_edpar *edi) |
933 | { |
934 | /* this function updates the parameter Efl and deltaF according to the rules given in |
935 | * 'predicting unimolecular chemical reactions: chemical flooding' M Mueller et al, |
936 | * J. chem Phys. */ |
937 | |
938 | if ((edi->flood.tau < 0 ? -edi->flood.tau : edi->flood.tau ) > 0.00000001) |
939 | { |
940 | edi->flood.Efl = edi->flood.Efl+edi->flood.dt/edi->flood.tau*(edi->flood.deltaF0-edi->flood.deltaF); |
941 | /* check if restrain (inverted flooding) -> don't let EFL become positive */ |
942 | if (edi->flood.alpha2 < 0 && edi->flood.Efl > -0.00000001) |
943 | { |
944 | edi->flood.Efl = 0; |
945 | } |
946 | |
947 | edi->flood.deltaF = (1-edi->flood.dt/edi->flood.tau)*edi->flood.deltaF+edi->flood.dt/edi->flood.tau*edi->flood.Vfl; |
948 | } |
949 | } |
950 | |
951 | |
952 | static void do_single_flood( |
953 | FILE *edo, |
954 | rvec x[], |
955 | rvec force[], |
956 | t_edpar *edi, |
957 | gmx_int64_t step, |
958 | matrix box, |
959 | t_commrec *cr, |
960 | gmx_bool bNS) /* Are we in a neighbor searching step? */ |
961 | { |
962 | int i; |
963 | matrix rotmat; /* rotation matrix */ |
964 | matrix tmat; /* inverse rotation */ |
965 | rvec transvec; /* translation vector */ |
966 | real rmsdev; |
967 | struct t_do_edsam *buf; |
968 | |
969 | |
970 | buf = edi->buf->do_edsam; |
971 | |
972 | |
973 | /* Broadcast the positions of the AVERAGE structure such that they are known on |
974 | * every processor. Each node contributes its local positions x and stores them in |
975 | * the collective ED array buf->xcoll */ |
976 | communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x, |
977 | edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box); |
978 | |
979 | /* Only assembly REFERENCE positions if their indices differ from the average ones */ |
980 | if (!edi->bRefEqAv) |
981 | { |
982 | communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x, |
983 | edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box); |
984 | } |
985 | |
986 | /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions. |
987 | * We do not need to update the shifts until the next NS step */ |
988 | buf->bUpdateShifts = FALSE0; |
989 | |
990 | /* Now all nodes have all of the ED/flooding positions in edi->sav->xcoll, |
991 | * as well as the indices in edi->sav.anrs */ |
992 | |
993 | /* Fit the reference indices to the reference structure */ |
994 | if (edi->bRefEqAv) |
995 | { |
996 | fit_to_reference(buf->xcoll, transvec, rotmat, edi); |
997 | } |
998 | else |
999 | { |
1000 | fit_to_reference(buf->xc_ref, transvec, rotmat, edi); |
1001 | } |
1002 | |
1003 | /* Now apply the translation and rotation to the ED structure */ |
1004 | translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat); |
1005 | |
1006 | /* Project fitted structure onto supbspace -> store in edi->flood.vecs.xproj */ |
1007 | project_to_eigvectors(buf->xcoll, &edi->flood.vecs, edi); |
1008 | |
1009 | if (FALSE0 == edi->flood.bConstForce) |
1010 | { |
1011 | /* Compute Vfl(x) from flood.xproj */ |
1012 | edi->flood.Vfl = flood_energy(edi, step); |
1013 | |
1014 | update_adaption(edi); |
1015 | |
1016 | /* Compute the flooding forces */ |
1017 | flood_forces(edi); |
1018 | } |
1019 | |
1020 | /* Translate them into cartesian positions */ |
1021 | flood_blowup(edi, edi->flood.forces_cartesian); |
1022 | |
1023 | /* Rotate forces back so that they correspond to the given structure and not to the fitted one */ |
1024 | /* Each node rotates back its local forces */ |
1025 | transpose(rotmat, tmat); |
1026 | rotate_x(edi->flood.forces_cartesian, edi->sav.nr_loc, tmat); |
1027 | |
1028 | /* Finally add forces to the main force variable */ |
1029 | for (i = 0; i < edi->sav.nr_loc; i++) |
1030 | { |
1031 | rvec_inc(force[edi->sav.anrs_loc[i]], edi->flood.forces_cartesian[i]); |
1032 | } |
1033 | |
1034 | /* Output is written by the master process */ |
1035 | if (do_per_step(step, edi->outfrq) && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
1036 | { |
1037 | /* Output how well we fit to the reference */ |
1038 | if (edi->bRefEqAv) |
1039 | { |
1040 | /* Indices of reference and average structures are identical, |
1041 | * thus we can calculate the rmsd to SREF using xcoll */ |
1042 | rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref); |
1043 | } |
1044 | else |
1045 | { |
1046 | /* We have to translate & rotate the reference atoms first */ |
1047 | translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat); |
1048 | rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref); |
1049 | } |
1050 | |
1051 | write_edo_flood(edi, edo, rmsdev); |
1052 | } |
1053 | } |
1054 | |
1055 | |
1056 | /* Main flooding routine, called from do_force */ |
1057 | extern void do_flood( |
1058 | t_commrec *cr, |
1059 | t_inputrec *ir, |
1060 | rvec x[], |
1061 | rvec force[], |
1062 | gmx_edsam_t ed, |
1063 | matrix box, |
1064 | gmx_int64_t step, |
1065 | gmx_bool bNS) |
1066 | { |
1067 | t_edpar *edi; |
1068 | |
1069 | |
1070 | edi = ed->edpar; |
1071 | |
1072 | /* Write time to edo, when required. Output the time anyhow since we need |
1073 | * it in the output file for ED constraints. */ |
1074 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && do_per_step(step, edi->outfrq)) |
1075 | { |
1076 | fprintf(ed->edo, "\n%12f", ir->init_t + step*ir->delta_t); |
1077 | } |
1078 | |
1079 | if (ed->eEDtype != eEDflood) |
1080 | { |
1081 | return; |
1082 | } |
1083 | |
1084 | while (edi) |
1085 | { |
1086 | /* Call flooding for one matrix */ |
1087 | if (edi->flood.vecs.neig) |
1088 | { |
1089 | do_single_flood(ed->edo, x, force, edi, step, box, cr, bNS); |
1090 | } |
1091 | edi = edi->next_edi; |
1092 | } |
1093 | } |
1094 | |
1095 | |
1096 | /* Called by init_edi, configure some flooding related variables and structures, |
1097 | * print headers to output files */ |
1098 | static void init_flood(t_edpar *edi, gmx_edsam_t ed, real dt) |
1099 | { |
1100 | int i; |
1101 | |
1102 | |
1103 | edi->flood.Efl = edi->flood.constEfl; |
1104 | edi->flood.Vfl = 0; |
1105 | edi->flood.dt = dt; |
1106 | |
1107 | if (edi->flood.vecs.neig) |
1108 | { |
1109 | /* If in any of the ED groups we find a flooding vector, flooding is turned on */ |
1110 | ed->eEDtype = eEDflood; |
1111 | |
1112 | fprintf(stderrstderr, "ED: Flooding %d eigenvector%s.\n", edi->flood.vecs.neig, edi->flood.vecs.neig > 1 ? "s" : ""); |
1113 | |
1114 | if (edi->flood.bConstForce) |
1115 | { |
1116 | /* We have used stpsz as a vehicle to carry the fproj values for constant |
1117 | * force flooding. Now we copy that to flood.vecs.fproj. Note that |
1118 | * in const force flooding, fproj is never changed. */ |
1119 | for (i = 0; i < edi->flood.vecs.neig; i++) |
1120 | { |
1121 | edi->flood.vecs.fproj[i] = edi->flood.vecs.stpsz[i]; |
1122 | |
1123 | fprintf(stderrstderr, "ED: applying on eigenvector %d a constant force of %g\n", |
1124 | edi->flood.vecs.ieig[i], edi->flood.vecs.fproj[i]); |
1125 | } |
1126 | } |
1127 | } |
1128 | } |
1129 | |
1130 | |
1131 | #ifdef DEBUGHELPERS |
1132 | /*********** Energy book keeping ******/ |
1133 | static void get_flood_enx_names(t_edpar *edi, char** names, int *nnames) /* get header of energies */ |
1134 | { |
1135 | t_edpar *actual; |
1136 | int count; |
1137 | char buf[STRLEN4096]; |
1138 | actual = edi; |
1139 | count = 1; |
1140 | while (actual) |
1141 | { |
1142 | srenew(names, count)(names) = save_realloc("names", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1142, (names), (count), sizeof(*(names))); |
1143 | sprintf(buf, "Vfl_%d", count); |
1144 | names[count-1] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t )(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1 ) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, buf, __len); __retval ; })) : __strdup (buf))); |
1145 | actual = actual->next_edi; |
1146 | count++; |
1147 | } |
1148 | *nnames = count-1; |
1149 | } |
1150 | |
1151 | |
1152 | static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames) |
1153 | { |
1154 | /*fl has to be big enough to capture nnames-many entries*/ |
1155 | t_edpar *actual; |
1156 | int count; |
1157 | |
1158 | |
1159 | actual = edi; |
1160 | count = 1; |
1161 | while (actual) |
1162 | { |
1163 | Vfl[count-1] = actual->flood.Vfl; |
1164 | actual = actual->next_edi; |
1165 | count++; |
1166 | } |
1167 | if (nnames != count-1) |
1168 | { |
1169 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1169, "Number of energies is not consistent with t_edi structure"); |
1170 | } |
1171 | } |
1172 | /************* END of FLOODING IMPLEMENTATION ****************************/ |
1173 | #endif |
1174 | |
1175 | |
1176 | gmx_edsam_t ed_open(int natoms, edsamstate_t *EDstate, int nfile, const t_filenm fnm[], unsigned long Flags, const output_env_t oenv, t_commrec *cr) |
1177 | { |
1178 | gmx_edsam_t ed; |
1179 | int nED; |
1180 | |
1181 | |
1182 | /* Allocate space for the ED data structure */ |
1183 | snew(ed, 1)(ed) = save_calloc("ed", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1183, (1), sizeof(*(ed))); |
1184 | |
1185 | /* We want to perform ED (this switch might later be upgraded to eEDflood) */ |
1186 | ed->eEDtype = eEDedsam; |
1187 | |
1188 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
1189 | { |
1190 | fprintf(stderrstderr, "ED sampling will be performed!\n"); |
1191 | snew(ed->edpar, 1)(ed->edpar) = save_calloc("ed->edpar", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1191, (1), sizeof(*(ed->edpar))); |
1192 | |
1193 | /* Read the edi input file: */ |
1194 | nED = read_edi_file(ftp2fn(efEDI, nfile, fnm), ed->edpar, natoms); |
1195 | |
1196 | /* Make sure the checkpoint was produced in a run using this .edi file */ |
1197 | if (EDstate->bFromCpt) |
1198 | { |
1199 | crosscheck_edi_file_vs_checkpoint(ed, EDstate); |
1200 | } |
1201 | else |
1202 | { |
1203 | EDstate->nED = nED; |
1204 | } |
1205 | init_edsamstate(ed, EDstate); |
1206 | |
1207 | /* The master opens the ED output file */ |
1208 | if (Flags & MD_APPENDFILES(1<<15)) |
1209 | { |
1210 | ed->edo = gmx_fio_fopen(opt2fn("-eo", nfile, fnm), "a+"); |
1211 | } |
1212 | else |
1213 | { |
1214 | ed->edo = xvgropen(opt2fn("-eo", nfile, fnm), |
1215 | "Essential dynamics / flooding output", |
1216 | "Time (ps)", |
1217 | "RMSDs (nm), projections on EVs (nm), ...", oenv); |
1218 | |
1219 | /* Make a descriptive legend */ |
1220 | write_edo_legend(ed, EDstate->nED, oenv); |
1221 | } |
1222 | } |
1223 | return ed; |
1224 | } |
1225 | |
1226 | |
1227 | /* Broadcasts the structure data */ |
1228 | static void bc_ed_positions(t_commrec *cr, struct gmx_edx *s, int stype) |
1229 | { |
1230 | snew_bc(cr, s->anrs, s->nr ){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((s->anrs)) = save_calloc("(s->anrs)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1230, ((s->nr)), sizeof(*((s->anrs)))); }}; /* Index numbers */ |
1231 | snew_bc(cr, s->x, s->nr ){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((s->x)) = save_calloc("(s->x)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1231, ((s->nr)), sizeof(*((s->x)))); }}; /* Positions */ |
1232 | nblock_bc(cr, s->nr, s->anrs )gmx_bcast((s->nr)*sizeof((s->anrs)[0]), (s->anrs), ( cr)); |
1233 | nblock_bc(cr, s->nr, s->x )gmx_bcast((s->nr)*sizeof((s->x)[0]), (s->x), (cr)); |
1234 | |
1235 | /* For the average & reference structures we need an array for the collective indices, |
1236 | * and we need to broadcast the masses as well */ |
1237 | if (stype == eedAV || stype == eedREF) |
1238 | { |
1239 | /* We need these additional variables in the parallel case: */ |
1240 | snew(s->c_ind, s->nr )(s->c_ind) = save_calloc("s->c_ind", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1240, (s->nr), sizeof(*(s->c_ind))); /* Collective indices */ |
1241 | /* Local atom indices get assigned in dd_make_local_group_indices. |
1242 | * There, also memory is allocated */ |
1243 | s->nalloc_loc = 0; /* allocation size of s->anrs_loc */ |
1244 | snew_bc(cr, s->x_old, s->nr){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((s->x_old)) = save_calloc("(s->x_old)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1244, ((s->nr)), sizeof(*((s->x_old)))); }}; /* To be able to always make the ED molecule whole, ... */ |
1245 | nblock_bc(cr, s->nr, s->x_old)gmx_bcast((s->nr)*sizeof((s->x_old)[0]), (s->x_old), (cr)); /* ... keep track of shift changes with the help of old coords */ |
1246 | } |
1247 | |
1248 | /* broadcast masses for the reference structure (for mass-weighted fitting) */ |
1249 | if (stype == eedREF) |
1250 | { |
1251 | snew_bc(cr, s->m, s->nr){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((s->m)) = save_calloc("(s->m)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1251, ((s->nr)), sizeof(*((s->m)))); }}; |
1252 | nblock_bc(cr, s->nr, s->m)gmx_bcast((s->nr)*sizeof((s->m)[0]), (s->m), (cr)); |
1253 | } |
1254 | |
1255 | /* For the average structure we might need the masses for mass-weighting */ |
1256 | if (stype == eedAV) |
1257 | { |
1258 | snew_bc(cr, s->sqrtm, s->nr){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((s->sqrtm)) = save_calloc("(s->sqrtm)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1258, ((s->nr)), sizeof(*((s->sqrtm)))); }}; |
1259 | nblock_bc(cr, s->nr, s->sqrtm)gmx_bcast((s->nr)*sizeof((s->sqrtm)[0]), (s->sqrtm), (cr)); |
1260 | snew_bc(cr, s->m, s->nr){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((s->m)) = save_calloc("(s->m)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1260, ((s->nr)), sizeof(*((s->m)))); }}; |
1261 | nblock_bc(cr, s->nr, s->m)gmx_bcast((s->nr)*sizeof((s->m)[0]), (s->m), (cr)); |
1262 | } |
1263 | } |
1264 | |
1265 | |
1266 | /* Broadcasts the eigenvector data */ |
1267 | static void bc_ed_vecs(t_commrec *cr, t_eigvec *ev, int length, gmx_bool bHarmonic) |
1268 | { |
1269 | int i; |
1270 | |
1271 | snew_bc(cr, ev->ieig, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->ieig)) = save_calloc("(ev->ieig)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1271, ((ev->neig)), sizeof(*((ev->ieig)))); }}; /* index numbers of eigenvector */ |
1272 | snew_bc(cr, ev->stpsz, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->stpsz)) = save_calloc("(ev->stpsz)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1272, ((ev->neig)), sizeof(*((ev->stpsz)))); }}; /* stepsizes per eigenvector */ |
1273 | snew_bc(cr, ev->xproj, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->xproj)) = save_calloc("(ev->xproj)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1273, ((ev->neig)), sizeof(*((ev->xproj)))); }}; /* instantaneous x projection */ |
1274 | snew_bc(cr, ev->fproj, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->fproj)) = save_calloc("(ev->fproj)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1274, ((ev->neig)), sizeof(*((ev->fproj)))); }}; /* instantaneous f projection */ |
1275 | snew_bc(cr, ev->refproj, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->refproj)) = save_calloc("(ev->refproj)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1275, ((ev->neig)), sizeof(*((ev->refproj)))); }}; /* starting or target projection */ |
1276 | |
1277 | nblock_bc(cr, ev->neig, ev->ieig )gmx_bcast((ev->neig)*sizeof((ev->ieig)[0]), (ev->ieig ), (cr)); |
1278 | nblock_bc(cr, ev->neig, ev->stpsz )gmx_bcast((ev->neig)*sizeof((ev->stpsz)[0]), (ev->stpsz ), (cr)); |
1279 | nblock_bc(cr, ev->neig, ev->xproj )gmx_bcast((ev->neig)*sizeof((ev->xproj)[0]), (ev->xproj ), (cr)); |
1280 | nblock_bc(cr, ev->neig, ev->fproj )gmx_bcast((ev->neig)*sizeof((ev->fproj)[0]), (ev->fproj ), (cr)); |
1281 | nblock_bc(cr, ev->neig, ev->refproj)gmx_bcast((ev->neig)*sizeof((ev->refproj)[0]), (ev-> refproj), (cr)); |
1282 | |
1283 | snew_bc(cr, ev->vec, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->vec)) = save_calloc("(ev->vec)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1283, ((ev->neig)), sizeof(*((ev->vec)))); }}; /* Eigenvector components */ |
1284 | for (i = 0; i < ev->neig; i++) |
1285 | { |
1286 | snew_bc(cr, ev->vec[i], length){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->vec[i])) = save_calloc("(ev->vec[i])", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1286, ((length)), sizeof(*((ev->vec[i])))); }}; |
1287 | nblock_bc(cr, length, ev->vec[i])gmx_bcast((length)*sizeof((ev->vec[i])[0]), (ev->vec[i] ), (cr)); |
1288 | } |
1289 | |
1290 | /* For harmonic restraints the reference projections can change with time */ |
1291 | if (bHarmonic) |
1292 | { |
1293 | snew_bc(cr, ev->refproj0, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->refproj0)) = save_calloc("(ev->refproj0)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1293, ((ev->neig)), sizeof(*((ev->refproj0)))); }}; |
1294 | snew_bc(cr, ev->refprojslope, ev->neig){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ev->refprojslope)) = save_calloc("(ev->refprojslope)" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1294, ((ev->neig)), sizeof(*((ev->refprojslope)))); } }; |
1295 | nblock_bc(cr, ev->neig, ev->refproj0 )gmx_bcast((ev->neig)*sizeof((ev->refproj0)[0]), (ev-> refproj0), (cr)); |
1296 | nblock_bc(cr, ev->neig, ev->refprojslope)gmx_bcast((ev->neig)*sizeof((ev->refprojslope)[0]), (ev ->refprojslope), (cr)); |
1297 | } |
1298 | } |
1299 | |
1300 | |
1301 | /* Broadcasts the ED / flooding data to other nodes |
1302 | * and allocates memory where needed */ |
1303 | static void broadcast_ed_data(t_commrec *cr, gmx_edsam_t ed, int numedis) |
1304 | { |
1305 | int nr; |
1306 | t_edpar *edi; |
1307 | |
1308 | |
1309 | /* Master lets the other nodes know if its ED only or also flooding */ |
1310 | gmx_bcast(sizeof(ed->eEDtype), &(ed->eEDtype), cr); |
1311 | |
1312 | snew_bc(cr, ed->edpar, 1){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((ed->edpar)) = save_calloc("(ed->edpar)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1312, ((1)), sizeof(*((ed->edpar)))); }}; |
1313 | /* Now transfer the ED data set(s) */ |
1314 | edi = ed->edpar; |
1315 | for (nr = 0; nr < numedis; nr++) |
1316 | { |
1317 | /* Broadcast a single ED data set */ |
1318 | block_bc(cr, *edi)gmx_bcast( sizeof(*edi), &(*edi), (cr)); |
1319 | |
1320 | /* Broadcast positions */ |
1321 | bc_ed_positions(cr, &(edi->sref), eedREF); /* reference positions (don't broadcast masses) */ |
1322 | bc_ed_positions(cr, &(edi->sav ), eedAV ); /* average positions (do broadcast masses as well) */ |
1323 | bc_ed_positions(cr, &(edi->star), eedTAR); /* target positions */ |
1324 | bc_ed_positions(cr, &(edi->sori), eedORI); /* origin positions */ |
1325 | |
1326 | /* Broadcast eigenvectors */ |
1327 | bc_ed_vecs(cr, &edi->vecs.mon, edi->sav.nr, FALSE0); |
1328 | bc_ed_vecs(cr, &edi->vecs.linfix, edi->sav.nr, FALSE0); |
1329 | bc_ed_vecs(cr, &edi->vecs.linacc, edi->sav.nr, FALSE0); |
1330 | bc_ed_vecs(cr, &edi->vecs.radfix, edi->sav.nr, FALSE0); |
1331 | bc_ed_vecs(cr, &edi->vecs.radacc, edi->sav.nr, FALSE0); |
1332 | bc_ed_vecs(cr, &edi->vecs.radcon, edi->sav.nr, FALSE0); |
1333 | /* Broadcast flooding eigenvectors and, if needed, values for the moving reference */ |
1334 | bc_ed_vecs(cr, &edi->flood.vecs, edi->sav.nr, edi->flood.bHarmonic); |
1335 | |
1336 | /* Set the pointer to the next ED group */ |
1337 | if (edi->next_edi) |
1338 | { |
1339 | snew_bc(cr, edi->next_edi, 1){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((edi->next_edi)) = save_calloc("(edi->next_edi)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1339, ((1)), sizeof(*((edi->next_edi)))); }}; |
1340 | edi = edi->next_edi; |
1341 | } |
1342 | } |
1343 | } |
1344 | |
1345 | |
1346 | /* init-routine called for every *.edi-cycle, initialises t_edpar structure */ |
1347 | static void init_edi(gmx_mtop_t *mtop, t_edpar *edi) |
1348 | { |
1349 | int i; |
1350 | real totalmass = 0.0; |
1351 | rvec com; |
1352 | gmx_mtop_atomlookup_t alook = NULL((void*)0); |
1353 | t_atom *atom; |
1354 | |
1355 | /* NOTE Init_edi is executed on the master process only |
1356 | * The initialized data sets are then transmitted to the |
1357 | * other nodes in broadcast_ed_data */ |
1358 | |
1359 | alook = gmx_mtop_atomlookup_init(mtop); |
1360 | |
1361 | /* evaluate masses (reference structure) */ |
1362 | snew(edi->sref.m, edi->sref.nr)(edi->sref.m) = save_calloc("edi->sref.m", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1362, (edi->sref.nr), sizeof(*(edi->sref.m))); |
1363 | for (i = 0; i < edi->sref.nr; i++) |
1364 | { |
1365 | if (edi->fitmas) |
1366 | { |
1367 | gmx_mtop_atomnr_to_atom(alook, edi->sref.anrs[i], &atom); |
1368 | edi->sref.m[i] = atom->m; |
1369 | } |
1370 | else |
1371 | { |
1372 | edi->sref.m[i] = 1.0; |
1373 | } |
1374 | |
1375 | /* Check that every m > 0. Bad things will happen otherwise. */ |
1376 | if (edi->sref.m[i] <= 0.0) |
1377 | { |
1378 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1378, "Reference structure atom %d (sam.edi index %d) has a mass of %g.\n" |
1379 | "For a mass-weighted fit, all reference structure atoms need to have a mass >0.\n" |
1380 | "Either make the covariance analysis non-mass-weighted, or exclude massless\n" |
1381 | "atoms from the reference structure by creating a proper index group.\n", |
1382 | i, edi->sref.anrs[i]+1, edi->sref.m[i]); |
1383 | } |
1384 | |
1385 | totalmass += edi->sref.m[i]; |
1386 | } |
1387 | edi->sref.mtot = totalmass; |
1388 | |
1389 | /* Masses m and sqrt(m) for the average structure. Note that m |
1390 | * is needed if forces have to be evaluated in do_edsam */ |
1391 | snew(edi->sav.sqrtm, edi->sav.nr )(edi->sav.sqrtm) = save_calloc("edi->sav.sqrtm", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1391, (edi->sav.nr), sizeof(*(edi->sav.sqrtm))); |
1392 | snew(edi->sav.m, edi->sav.nr )(edi->sav.m) = save_calloc("edi->sav.m", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1392, (edi->sav.nr), sizeof(*(edi->sav.m))); |
1393 | for (i = 0; i < edi->sav.nr; i++) |
1394 | { |
1395 | gmx_mtop_atomnr_to_atom(alook, edi->sav.anrs[i], &atom); |
1396 | edi->sav.m[i] = atom->m; |
1397 | if (edi->pcamas) |
1398 | { |
1399 | edi->sav.sqrtm[i] = sqrt(atom->m); |
1400 | } |
1401 | else |
1402 | { |
1403 | edi->sav.sqrtm[i] = 1.0; |
1404 | } |
1405 | |
1406 | /* Check that every m > 0. Bad things will happen otherwise. */ |
1407 | if (edi->sav.sqrtm[i] <= 0.0) |
1408 | { |
1409 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1409, "Average structure atom %d (sam.edi index %d) has a mass of %g.\n" |
1410 | "For ED with mass-weighting, all average structure atoms need to have a mass >0.\n" |
1411 | "Either make the covariance analysis non-mass-weighted, or exclude massless\n" |
1412 | "atoms from the average structure by creating a proper index group.\n", |
1413 | i, edi->sav.anrs[i]+1, atom->m); |
1414 | } |
1415 | } |
1416 | |
1417 | gmx_mtop_atomlookup_destroy(alook); |
1418 | |
1419 | /* put reference structure in origin */ |
1420 | get_center(edi->sref.x, edi->sref.m, edi->sref.nr, com); |
1421 | com[XX0] = -com[XX0]; |
1422 | com[YY1] = -com[YY1]; |
1423 | com[ZZ2] = -com[ZZ2]; |
1424 | translate_x(edi->sref.x, edi->sref.nr, com); |
1425 | |
1426 | /* Init ED buffer */ |
1427 | snew(edi->buf, 1)(edi->buf) = save_calloc("edi->buf", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1427, (1), sizeof(*(edi->buf))); |
1428 | } |
1429 | |
1430 | |
1431 | static void check(const char *line, const char *label) |
1432 | { |
1433 | if (!strstr(line, label)) |
1434 | { |
1435 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1435, "Could not find input parameter %s at expected position in edsam input-file (.edi)\nline read instead is %s", label, line); |
1436 | } |
1437 | } |
1438 | |
1439 | |
1440 | static int read_checked_edint(FILE *file, const char *label) |
1441 | { |
1442 | char line[STRLEN4096+1]; |
1443 | int idum; |
1444 | |
1445 | |
1446 | fgets2 (line, STRLEN4096, file); |
1447 | check(line, label); |
1448 | fgets2 (line, STRLEN4096, file); |
1449 | sscanf (line, "%d", &idum); |
1450 | return idum; |
1451 | } |
1452 | |
1453 | |
1454 | static int read_edint(FILE *file, gmx_bool *bEOF) |
1455 | { |
1456 | char line[STRLEN4096+1]; |
1457 | int idum; |
1458 | char *eof; |
1459 | |
1460 | |
1461 | eof = fgets2 (line, STRLEN4096, file); |
1462 | if (eof == NULL((void*)0)) |
1463 | { |
1464 | *bEOF = TRUE1; |
1465 | return -1; |
1466 | } |
1467 | eof = fgets2 (line, STRLEN4096, file); |
1468 | if (eof == NULL((void*)0)) |
1469 | { |
1470 | *bEOF = TRUE1; |
1471 | return -1; |
1472 | } |
1473 | sscanf (line, "%d", &idum); |
1474 | *bEOF = FALSE0; |
1475 | return idum; |
1476 | } |
1477 | |
1478 | |
1479 | static real read_checked_edreal(FILE *file, const char *label) |
1480 | { |
1481 | char line[STRLEN4096+1]; |
1482 | double rdum; |
1483 | |
1484 | |
1485 | fgets2 (line, STRLEN4096, file); |
1486 | check(line, label); |
1487 | fgets2 (line, STRLEN4096, file); |
1488 | sscanf (line, "%lf", &rdum); |
1489 | return (real) rdum; /* always read as double and convert to single */ |
1490 | } |
1491 | |
1492 | |
1493 | static void read_edx(FILE *file, int number, int *anrs, rvec *x) |
1494 | { |
1495 | int i, j; |
1496 | char line[STRLEN4096+1]; |
1497 | double d[3]; |
1498 | |
1499 | |
1500 | for (i = 0; i < number; i++) |
1501 | { |
1502 | fgets2 (line, STRLEN4096, file); |
1503 | sscanf (line, "%d%lf%lf%lf", &anrs[i], &d[0], &d[1], &d[2]); |
1504 | anrs[i]--; /* we are reading FORTRAN indices */ |
1505 | for (j = 0; j < 3; j++) |
1506 | { |
1507 | x[i][j] = d[j]; /* always read as double and convert to single */ |
1508 | } |
1509 | } |
1510 | } |
1511 | |
1512 | |
1513 | static void scan_edvec(FILE *in, int nr, rvec *vec) |
1514 | { |
1515 | char line[STRLEN4096+1]; |
1516 | int i; |
1517 | double x, y, z; |
1518 | |
1519 | |
1520 | for (i = 0; (i < nr); i++) |
1521 | { |
1522 | fgets2 (line, STRLEN4096, in); |
1523 | sscanf (line, "%le%le%le", &x, &y, &z); |
1524 | vec[i][XX0] = x; |
1525 | vec[i][YY1] = y; |
1526 | vec[i][ZZ2] = z; |
1527 | } |
1528 | } |
1529 | |
1530 | |
1531 | static void read_edvec(FILE *in, int nr, t_eigvec *tvec, gmx_bool bReadRefproj, gmx_bool *bHaveReference) |
1532 | { |
1533 | int i, idum, nscan; |
1534 | double rdum, refproj_dum = 0.0, refprojslope_dum = 0.0; |
1535 | char line[STRLEN4096+1]; |
1536 | |
1537 | |
1538 | tvec->neig = read_checked_edint(in, "NUMBER OF EIGENVECTORS"); |
1539 | if (tvec->neig > 0) |
1540 | { |
1541 | snew(tvec->ieig, tvec->neig)(tvec->ieig) = save_calloc("tvec->ieig", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1541, (tvec->neig), sizeof(*(tvec->ieig))); |
1542 | snew(tvec->stpsz, tvec->neig)(tvec->stpsz) = save_calloc("tvec->stpsz", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1542, (tvec->neig), sizeof(*(tvec->stpsz))); |
1543 | snew(tvec->vec, tvec->neig)(tvec->vec) = save_calloc("tvec->vec", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1543, (tvec->neig), sizeof(*(tvec->vec))); |
1544 | snew(tvec->xproj, tvec->neig)(tvec->xproj) = save_calloc("tvec->xproj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1544, (tvec->neig), sizeof(*(tvec->xproj))); |
1545 | snew(tvec->fproj, tvec->neig)(tvec->fproj) = save_calloc("tvec->fproj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1545, (tvec->neig), sizeof(*(tvec->fproj))); |
1546 | snew(tvec->refproj, tvec->neig)(tvec->refproj) = save_calloc("tvec->refproj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1546, (tvec->neig), sizeof(*(tvec->refproj))); |
1547 | if (bReadRefproj) |
1548 | { |
1549 | snew(tvec->refproj0, tvec->neig)(tvec->refproj0) = save_calloc("tvec->refproj0", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1549, (tvec->neig), sizeof(*(tvec->refproj0))); |
1550 | snew(tvec->refprojslope, tvec->neig)(tvec->refprojslope) = save_calloc("tvec->refprojslope" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1550, (tvec->neig), sizeof(*(tvec->refprojslope))); |
1551 | } |
1552 | |
1553 | for (i = 0; (i < tvec->neig); i++) |
1554 | { |
1555 | fgets2 (line, STRLEN4096, in); |
1556 | if (bReadRefproj) /* ONLY when using flooding as harmonic restraint */ |
1557 | { |
1558 | nscan = sscanf(line, "%d%lf%lf%lf", &idum, &rdum, &refproj_dum, &refprojslope_dum); |
1559 | /* Zero out values which were not scanned */ |
1560 | switch (nscan) |
1561 | { |
1562 | case 4: |
1563 | /* Every 4 values read, including reference position */ |
1564 | *bHaveReference = TRUE1; |
1565 | break; |
1566 | case 3: |
1567 | /* A reference position is provided */ |
1568 | *bHaveReference = TRUE1; |
1569 | /* No value for slope, set to 0 */ |
1570 | refprojslope_dum = 0.0; |
1571 | break; |
1572 | case 2: |
1573 | /* No values for reference projection and slope, set to 0 */ |
1574 | refproj_dum = 0.0; |
1575 | refprojslope_dum = 0.0; |
1576 | break; |
1577 | default: |
1578 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1578, "Expected 2 - 4 (not %d) values for flooding vec: <nr> <spring const> <refproj> <refproj-slope>\n", nscan); |
1579 | break; |
1580 | } |
1581 | tvec->refproj[i] = refproj_dum; |
1582 | tvec->refproj0[i] = refproj_dum; |
1583 | tvec->refprojslope[i] = refprojslope_dum; |
1584 | } |
1585 | else /* Normal flooding */ |
1586 | { |
1587 | nscan = sscanf(line, "%d%lf", &idum, &rdum); |
1588 | if (nscan != 2) |
1589 | { |
1590 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1590, "Expected 2 values for flooding vec: <nr> <stpsz>\n"); |
1591 | } |
1592 | } |
1593 | tvec->ieig[i] = idum; |
1594 | tvec->stpsz[i] = rdum; |
1595 | } /* end of loop over eigenvectors */ |
1596 | |
1597 | for (i = 0; (i < tvec->neig); i++) |
1598 | { |
1599 | snew(tvec->vec[i], nr)(tvec->vec[i]) = save_calloc("tvec->vec[i]", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1599, (nr), sizeof(*(tvec->vec[i]))); |
1600 | scan_edvec(in, nr, tvec->vec[i]); |
1601 | } |
1602 | } |
1603 | } |
1604 | |
1605 | |
1606 | /* calls read_edvec for the vector groups, only for flooding there is an extra call */ |
1607 | static void read_edvecs(FILE *in, int nr, t_edvecs *vecs) |
1608 | { |
1609 | gmx_bool bHaveReference = FALSE0; |
1610 | |
1611 | |
1612 | read_edvec(in, nr, &vecs->mon, FALSE0, &bHaveReference); |
1613 | read_edvec(in, nr, &vecs->linfix, FALSE0, &bHaveReference); |
1614 | read_edvec(in, nr, &vecs->linacc, FALSE0, &bHaveReference); |
1615 | read_edvec(in, nr, &vecs->radfix, FALSE0, &bHaveReference); |
1616 | read_edvec(in, nr, &vecs->radacc, FALSE0, &bHaveReference); |
1617 | read_edvec(in, nr, &vecs->radcon, FALSE0, &bHaveReference); |
1618 | } |
1619 | |
1620 | |
1621 | /* Check if the same atom indices are used for reference and average positions */ |
1622 | static gmx_bool check_if_same(struct gmx_edx sref, struct gmx_edx sav) |
1623 | { |
1624 | int i; |
1625 | |
1626 | |
1627 | /* If the number of atoms differs between the two structures, |
1628 | * they cannot be identical */ |
1629 | if (sref.nr != sav.nr) |
1630 | { |
1631 | return FALSE0; |
1632 | } |
1633 | |
1634 | /* Now that we know that both stuctures have the same number of atoms, |
1635 | * check if also the indices are identical */ |
1636 | for (i = 0; i < sav.nr; i++) |
1637 | { |
1638 | if (sref.anrs[i] != sav.anrs[i]) |
1639 | { |
1640 | return FALSE0; |
1641 | } |
1642 | } |
1643 | fprintf(stderrstderr, "ED: Note: Reference and average structure are composed of the same atom indices.\n"); |
1644 | |
1645 | return TRUE1; |
1646 | } |
1647 | |
1648 | |
1649 | static int read_edi(FILE* in, t_edpar *edi, int nr_mdatoms, const char *fn) |
1650 | { |
1651 | int readmagic; |
1652 | const int magic = 670; |
1653 | gmx_bool bEOF; |
1654 | |
1655 | /* Was a specific reference point for the flooding/umbrella potential provided in the edi file? */ |
1656 | gmx_bool bHaveReference = FALSE0; |
1657 | |
1658 | |
1659 | /* the edi file is not free format, so expect problems if the input is corrupt. */ |
1660 | |
1661 | /* check the magic number */ |
1662 | readmagic = read_edint(in, &bEOF); |
1663 | /* Check whether we have reached the end of the input file */ |
1664 | if (bEOF) |
1665 | { |
1666 | return 0; |
1667 | } |
1668 | |
1669 | if (readmagic != magic) |
1670 | { |
1671 | if (readmagic == 666 || readmagic == 667 || readmagic == 668) |
1672 | { |
1673 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1673, "Wrong magic number: Use newest version of make_edi to produce edi file"); |
1674 | } |
1675 | else if (readmagic != 669) |
1676 | { |
1677 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1677, "Wrong magic number %d in %s", readmagic, fn); |
1678 | } |
1679 | } |
1680 | |
1681 | /* check the number of atoms */ |
1682 | edi->nini = read_edint(in, &bEOF); |
1683 | if (edi->nini != nr_mdatoms) |
1684 | { |
1685 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1685, "Nr of atoms in %s (%d) does not match nr of md atoms (%d)", fn, edi->nini, nr_mdatoms); |
1686 | } |
1687 | |
1688 | /* Done checking. For the rest we blindly trust the input */ |
1689 | edi->fitmas = read_checked_edint(in, "FITMAS"); |
1690 | edi->pcamas = read_checked_edint(in, "ANALYSIS_MAS"); |
1691 | edi->outfrq = read_checked_edint(in, "OUTFRQ"); |
1692 | edi->maxedsteps = read_checked_edint(in, "MAXLEN"); |
1693 | edi->slope = read_checked_edreal(in, "SLOPECRIT"); |
1694 | |
1695 | edi->presteps = read_checked_edint(in, "PRESTEPS"); |
1696 | edi->flood.deltaF0 = read_checked_edreal(in, "DELTA_F0"); |
1697 | edi->flood.deltaF = read_checked_edreal(in, "INIT_DELTA_F"); |
1698 | edi->flood.tau = read_checked_edreal(in, "TAU"); |
1699 | edi->flood.constEfl = read_checked_edreal(in, "EFL_NULL"); |
1700 | edi->flood.alpha2 = read_checked_edreal(in, "ALPHA2"); |
1701 | edi->flood.kT = read_checked_edreal(in, "KT"); |
1702 | edi->flood.bHarmonic = read_checked_edint(in, "HARMONIC"); |
1703 | if (readmagic > 669) |
1704 | { |
1705 | edi->flood.bConstForce = read_checked_edint(in, "CONST_FORCE_FLOODING"); |
1706 | } |
1707 | else |
1708 | { |
1709 | edi->flood.bConstForce = FALSE0; |
1710 | } |
1711 | edi->sref.nr = read_checked_edint(in, "NREF"); |
1712 | |
1713 | /* allocate space for reference positions and read them */ |
1714 | snew(edi->sref.anrs, edi->sref.nr)(edi->sref.anrs) = save_calloc("edi->sref.anrs", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1714, (edi->sref.nr), sizeof(*(edi->sref.anrs))); |
1715 | snew(edi->sref.x, edi->sref.nr)(edi->sref.x) = save_calloc("edi->sref.x", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1715, (edi->sref.nr), sizeof(*(edi->sref.x))); |
1716 | snew(edi->sref.x_old, edi->sref.nr)(edi->sref.x_old) = save_calloc("edi->sref.x_old", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1716, (edi->sref.nr), sizeof(*(edi->sref.x_old))); |
1717 | edi->sref.sqrtm = NULL((void*)0); |
1718 | read_edx(in, edi->sref.nr, edi->sref.anrs, edi->sref.x); |
1719 | |
1720 | /* average positions. they define which atoms will be used for ED sampling */ |
1721 | edi->sav.nr = read_checked_edint(in, "NAV"); |
1722 | snew(edi->sav.anrs, edi->sav.nr)(edi->sav.anrs) = save_calloc("edi->sav.anrs", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1722, (edi->sav.nr), sizeof(*(edi->sav.anrs))); |
1723 | snew(edi->sav.x, edi->sav.nr)(edi->sav.x) = save_calloc("edi->sav.x", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1723, (edi->sav.nr), sizeof(*(edi->sav.x))); |
1724 | snew(edi->sav.x_old, edi->sav.nr)(edi->sav.x_old) = save_calloc("edi->sav.x_old", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1724, (edi->sav.nr), sizeof(*(edi->sav.x_old))); |
1725 | read_edx(in, edi->sav.nr, edi->sav.anrs, edi->sav.x); |
1726 | |
1727 | /* Check if the same atom indices are used for reference and average positions */ |
1728 | edi->bRefEqAv = check_if_same(edi->sref, edi->sav); |
1729 | |
1730 | /* eigenvectors */ |
1731 | read_edvecs(in, edi->sav.nr, &edi->vecs); |
1732 | read_edvec(in, edi->sav.nr, &edi->flood.vecs, edi->flood.bHarmonic, &bHaveReference); |
1733 | |
1734 | /* target positions */ |
1735 | edi->star.nr = read_edint(in, &bEOF); |
1736 | if (edi->star.nr > 0) |
1737 | { |
1738 | snew(edi->star.anrs, edi->star.nr)(edi->star.anrs) = save_calloc("edi->star.anrs", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1738, (edi->star.nr), sizeof(*(edi->star.anrs))); |
1739 | snew(edi->star.x, edi->star.nr)(edi->star.x) = save_calloc("edi->star.x", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1739, (edi->star.nr), sizeof(*(edi->star.x))); |
1740 | edi->star.sqrtm = NULL((void*)0); |
1741 | read_edx(in, edi->star.nr, edi->star.anrs, edi->star.x); |
1742 | } |
1743 | |
1744 | /* positions defining origin of expansion circle */ |
1745 | edi->sori.nr = read_edint(in, &bEOF); |
1746 | if (edi->sori.nr > 0) |
1747 | { |
1748 | if (bHaveReference) |
1749 | { |
1750 | /* Both an -ori structure and a at least one manual reference point have been |
1751 | * specified. That's ambiguous and probably not intentional. */ |
1752 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1752, "ED: An origin structure has been provided and a at least one (moving) reference\n" |
1753 | " point was manually specified in the edi file. That is ambiguous. Aborting.\n"); |
1754 | } |
1755 | snew(edi->sori.anrs, edi->sori.nr)(edi->sori.anrs) = save_calloc("edi->sori.anrs", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1755, (edi->sori.nr), sizeof(*(edi->sori.anrs))); |
1756 | snew(edi->sori.x, edi->sori.nr)(edi->sori.x) = save_calloc("edi->sori.x", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1756, (edi->sori.nr), sizeof(*(edi->sori.x))); |
1757 | edi->sori.sqrtm = NULL((void*)0); |
1758 | read_edx(in, edi->sori.nr, edi->sori.anrs, edi->sori.x); |
1759 | } |
1760 | |
1761 | /* all done */ |
1762 | return 1; |
1763 | } |
1764 | |
1765 | |
1766 | |
1767 | /* Read in the edi input file. Note that it may contain several ED data sets which were |
1768 | * achieved by concatenating multiple edi files. The standard case would be a single ED |
1769 | * data set, though. */ |
1770 | static int read_edi_file(const char *fn, t_edpar *edi, int nr_mdatoms) |
1771 | { |
1772 | FILE *in; |
1773 | t_edpar *curr_edi, *last_edi; |
1774 | t_edpar *edi_read; |
1775 | int edi_nr = 0; |
1776 | |
1777 | |
1778 | /* This routine is executed on the master only */ |
1779 | |
1780 | /* Open the .edi parameter input file */ |
1781 | in = gmx_fio_fopen(fn, "r"); |
1782 | fprintf(stderrstderr, "ED: Reading edi file %s\n", fn); |
1783 | |
1784 | /* Now read a sequence of ED input parameter sets from the edi file */ |
1785 | curr_edi = edi; |
1786 | last_edi = edi; |
1787 | while (read_edi(in, curr_edi, nr_mdatoms, fn) ) |
1788 | { |
1789 | edi_nr++; |
1790 | |
1791 | /* Since we arrived within this while loop we know that there is still another data set to be read in */ |
1792 | /* We need to allocate space for the data: */ |
1793 | snew(edi_read, 1)(edi_read) = save_calloc("edi_read", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1793, (1), sizeof(*(edi_read))); |
1794 | /* Point the 'next_edi' entry to the next edi: */ |
1795 | curr_edi->next_edi = edi_read; |
1796 | /* Keep the curr_edi pointer for the case that the next group is empty: */ |
1797 | last_edi = curr_edi; |
1798 | /* Let's prepare to read in the next edi data set: */ |
1799 | curr_edi = edi_read; |
1800 | } |
1801 | if (edi_nr == 0) |
1802 | { |
1803 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1803, "No complete ED data set found in edi file %s.", fn); |
1804 | } |
1805 | |
1806 | /* Terminate the edi group list with a NULL pointer: */ |
1807 | last_edi->next_edi = NULL((void*)0); |
1808 | |
1809 | fprintf(stderrstderr, "ED: Found %d ED group%s.\n", edi_nr, edi_nr > 1 ? "s" : ""); |
1810 | |
1811 | /* Close the .edi file again */ |
1812 | gmx_fio_fclose(in); |
1813 | |
1814 | return edi_nr; |
1815 | } |
1816 | |
1817 | |
1818 | struct t_fit_to_ref { |
1819 | rvec *xcopy; /* Working copy of the positions in fit_to_reference */ |
1820 | }; |
1821 | |
1822 | /* Fit the current positions to the reference positions |
1823 | * Do not actually do the fit, just return rotation and translation. |
1824 | * Note that the COM of the reference structure was already put into |
1825 | * the origin by init_edi. */ |
1826 | static void fit_to_reference(rvec *xcoll, /* The positions to be fitted */ |
1827 | rvec transvec, /* The translation vector */ |
1828 | matrix rotmat, /* The rotation matrix */ |
1829 | t_edpar *edi) /* Just needed for do_edfit */ |
1830 | { |
1831 | rvec com; /* center of mass */ |
1832 | int i; |
1833 | struct t_fit_to_ref *loc; |
1834 | |
1835 | |
1836 | /* Allocate memory the first time this routine is called for each edi group */ |
1837 | if (NULL((void*)0) == edi->buf->fit_to_ref) |
1838 | { |
1839 | snew(edi->buf->fit_to_ref, 1)(edi->buf->fit_to_ref) = save_calloc("edi->buf->fit_to_ref" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1839, (1), sizeof(*(edi->buf->fit_to_ref))); |
1840 | snew(edi->buf->fit_to_ref->xcopy, edi->sref.nr)(edi->buf->fit_to_ref->xcopy) = save_calloc("edi->buf->fit_to_ref->xcopy" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 1840, (edi->sref.nr), sizeof(*(edi->buf->fit_to_ref ->xcopy))); |
1841 | } |
1842 | loc = edi->buf->fit_to_ref; |
1843 | |
1844 | /* We do not touch the original positions but work on a copy. */ |
1845 | for (i = 0; i < edi->sref.nr; i++) |
1846 | { |
1847 | copy_rvec(xcoll[i], loc->xcopy[i]); |
1848 | } |
1849 | |
1850 | /* Calculate the center of mass */ |
1851 | get_center(loc->xcopy, edi->sref.m, edi->sref.nr, com); |
1852 | |
1853 | transvec[XX0] = -com[XX0]; |
1854 | transvec[YY1] = -com[YY1]; |
1855 | transvec[ZZ2] = -com[ZZ2]; |
1856 | |
1857 | /* Subtract the center of mass from the copy */ |
1858 | translate_x(loc->xcopy, edi->sref.nr, transvec); |
1859 | |
1860 | /* Determine the rotation matrix */ |
1861 | do_edfit(edi->sref.nr, edi->sref.x, loc->xcopy, rotmat, edi); |
1862 | } |
1863 | |
1864 | |
1865 | static void translate_and_rotate(rvec *x, /* The positions to be translated and rotated */ |
1866 | int nat, /* How many positions are there? */ |
1867 | rvec transvec, /* The translation vector */ |
1868 | matrix rotmat) /* The rotation matrix */ |
1869 | { |
1870 | /* Translation */ |
1871 | translate_x(x, nat, transvec); |
1872 | |
1873 | /* Rotation */ |
1874 | rotate_x(x, nat, rotmat); |
1875 | } |
1876 | |
1877 | |
1878 | /* Gets the rms deviation of the positions to the structure s */ |
1879 | /* fit_to_structure has to be called before calling this routine! */ |
1880 | static real rmsd_from_structure(rvec *x, /* The positions under consideration */ |
1881 | struct gmx_edx *s) /* The structure from which the rmsd shall be computed */ |
1882 | { |
1883 | real rmsd = 0.0; |
1884 | int i; |
1885 | |
1886 | |
1887 | for (i = 0; i < s->nr; i++) |
1888 | { |
1889 | rmsd += distance2(s->x[i], x[i]); |
1890 | } |
1891 | |
1892 | rmsd /= (real) s->nr; |
1893 | rmsd = sqrt(rmsd); |
1894 | |
1895 | return rmsd; |
1896 | } |
1897 | |
1898 | |
1899 | void dd_make_local_ed_indices(gmx_domdec_t *dd, struct gmx_edsam *ed) |
1900 | { |
1901 | t_edpar *edi; |
1902 | |
1903 | |
1904 | if (ed->eEDtype != eEDnone) |
1905 | { |
1906 | /* Loop over ED groups */ |
1907 | edi = ed->edpar; |
1908 | while (edi) |
1909 | { |
1910 | /* Local atoms of the reference structure (for fitting), need only be assembled |
1911 | * if their indices differ from the average ones */ |
1912 | if (!edi->bRefEqAv) |
1913 | { |
1914 | dd_make_local_group_indices(dd->ga2la, edi->sref.nr, edi->sref.anrs, |
1915 | &edi->sref.nr_loc, &edi->sref.anrs_loc, &edi->sref.nalloc_loc, edi->sref.c_ind); |
1916 | } |
1917 | |
1918 | /* Local atoms of the average structure (on these ED will be performed) */ |
1919 | dd_make_local_group_indices(dd->ga2la, edi->sav.nr, edi->sav.anrs, |
1920 | &edi->sav.nr_loc, &edi->sav.anrs_loc, &edi->sav.nalloc_loc, edi->sav.c_ind); |
1921 | |
1922 | /* Indicate that the ED shift vectors for this structure need to be updated |
1923 | * at the next call to communicate_group_positions, since obviously we are in a NS step */ |
1924 | edi->buf->do_edsam->bUpdateShifts = TRUE1; |
1925 | |
1926 | /* Set the pointer to the next ED group (if any) */ |
1927 | edi = edi->next_edi; |
1928 | } |
1929 | } |
1930 | } |
1931 | |
1932 | |
1933 | static gmx_inlineinline void ed_unshift_single_coord(matrix box, const rvec x, const ivec is, rvec xu) |
1934 | { |
1935 | int tx, ty, tz; |
1936 | |
1937 | |
1938 | tx = is[XX0]; |
1939 | ty = is[YY1]; |
1940 | tz = is[ZZ2]; |
1941 | |
1942 | if (TRICLINIC(box)(box[1][0] != 0 || box[2][0] != 0 || box[2][1] != 0)) |
1943 | { |
1944 | xu[XX0] = x[XX0]-tx*box[XX0][XX0]-ty*box[YY1][XX0]-tz*box[ZZ2][XX0]; |
1945 | xu[YY1] = x[YY1]-ty*box[YY1][YY1]-tz*box[ZZ2][YY1]; |
1946 | xu[ZZ2] = x[ZZ2]-tz*box[ZZ2][ZZ2]; |
1947 | } |
1948 | else |
1949 | { |
1950 | xu[XX0] = x[XX0]-tx*box[XX0][XX0]; |
1951 | xu[YY1] = x[YY1]-ty*box[YY1][YY1]; |
1952 | xu[ZZ2] = x[ZZ2]-tz*box[ZZ2][ZZ2]; |
1953 | } |
1954 | } |
1955 | |
1956 | |
1957 | static void do_linfix(rvec *xcoll, t_edpar *edi, gmx_int64_t step) |
1958 | { |
1959 | int i, j; |
1960 | real proj, add; |
1961 | rvec vec_dum; |
1962 | |
1963 | |
1964 | /* loop over linfix vectors */ |
1965 | for (i = 0; i < edi->vecs.linfix.neig; i++) |
1966 | { |
1967 | /* calculate the projection */ |
1968 | proj = projectx(edi, xcoll, edi->vecs.linfix.vec[i]); |
1969 | |
1970 | /* calculate the correction */ |
1971 | add = edi->vecs.linfix.refproj[i] + step*edi->vecs.linfix.stpsz[i] - proj; |
1972 | |
1973 | /* apply the correction */ |
1974 | add /= edi->sav.sqrtm[i]; |
1975 | for (j = 0; j < edi->sav.nr; j++) |
1976 | { |
1977 | svmul(add, edi->vecs.linfix.vec[i][j], vec_dum); |
1978 | rvec_inc(xcoll[j], vec_dum); |
1979 | } |
1980 | } |
1981 | } |
1982 | |
1983 | |
1984 | static void do_linacc(rvec *xcoll, t_edpar *edi) |
1985 | { |
1986 | int i, j; |
1987 | real proj, add; |
1988 | rvec vec_dum; |
1989 | |
1990 | |
1991 | /* loop over linacc vectors */ |
1992 | for (i = 0; i < edi->vecs.linacc.neig; i++) |
1993 | { |
1994 | /* calculate the projection */ |
1995 | proj = projectx(edi, xcoll, edi->vecs.linacc.vec[i]); |
1996 | |
1997 | /* calculate the correction */ |
1998 | add = 0.0; |
1999 | if (edi->vecs.linacc.stpsz[i] > 0.0) |
2000 | { |
2001 | if ((proj-edi->vecs.linacc.refproj[i]) < 0.0) |
2002 | { |
2003 | add = edi->vecs.linacc.refproj[i] - proj; |
2004 | } |
2005 | } |
2006 | if (edi->vecs.linacc.stpsz[i] < 0.0) |
2007 | { |
2008 | if ((proj-edi->vecs.linacc.refproj[i]) > 0.0) |
2009 | { |
2010 | add = edi->vecs.linacc.refproj[i] - proj; |
2011 | } |
2012 | } |
2013 | |
2014 | /* apply the correction */ |
2015 | add /= edi->sav.sqrtm[i]; |
2016 | for (j = 0; j < edi->sav.nr; j++) |
2017 | { |
2018 | svmul(add, edi->vecs.linacc.vec[i][j], vec_dum); |
2019 | rvec_inc(xcoll[j], vec_dum); |
2020 | } |
2021 | |
2022 | /* new positions will act as reference */ |
2023 | edi->vecs.linacc.refproj[i] = proj + add; |
2024 | } |
2025 | } |
2026 | |
2027 | |
2028 | static void do_radfix(rvec *xcoll, t_edpar *edi) |
2029 | { |
2030 | int i, j; |
2031 | real *proj, rad = 0.0, ratio; |
2032 | rvec vec_dum; |
2033 | |
2034 | |
2035 | if (edi->vecs.radfix.neig == 0) |
2036 | { |
2037 | return; |
2038 | } |
2039 | |
2040 | snew(proj, edi->vecs.radfix.neig)(proj) = save_calloc("proj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2040, (edi->vecs.radfix.neig), sizeof(*(proj))); |
2041 | |
2042 | /* loop over radfix vectors */ |
2043 | for (i = 0; i < edi->vecs.radfix.neig; i++) |
2044 | { |
2045 | /* calculate the projections, radius */ |
2046 | proj[i] = projectx(edi, xcoll, edi->vecs.radfix.vec[i]); |
2047 | rad += pow(proj[i] - edi->vecs.radfix.refproj[i], 2); |
2048 | } |
2049 | |
2050 | rad = sqrt(rad); |
2051 | ratio = (edi->vecs.radfix.stpsz[0]+edi->vecs.radfix.radius)/rad - 1.0; |
2052 | edi->vecs.radfix.radius += edi->vecs.radfix.stpsz[0]; |
2053 | |
2054 | /* loop over radfix vectors */ |
2055 | for (i = 0; i < edi->vecs.radfix.neig; i++) |
2056 | { |
2057 | proj[i] -= edi->vecs.radfix.refproj[i]; |
2058 | |
2059 | /* apply the correction */ |
2060 | proj[i] /= edi->sav.sqrtm[i]; |
2061 | proj[i] *= ratio; |
2062 | for (j = 0; j < edi->sav.nr; j++) |
2063 | { |
2064 | svmul(proj[i], edi->vecs.radfix.vec[i][j], vec_dum); |
2065 | rvec_inc(xcoll[j], vec_dum); |
2066 | } |
2067 | } |
2068 | |
2069 | sfree(proj)save_free("proj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2069, (proj)); |
2070 | } |
2071 | |
2072 | |
2073 | static void do_radacc(rvec *xcoll, t_edpar *edi) |
2074 | { |
2075 | int i, j; |
2076 | real *proj, rad = 0.0, ratio = 0.0; |
2077 | rvec vec_dum; |
2078 | |
2079 | |
2080 | if (edi->vecs.radacc.neig == 0) |
2081 | { |
2082 | return; |
2083 | } |
2084 | |
2085 | snew(proj, edi->vecs.radacc.neig)(proj) = save_calloc("proj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2085, (edi->vecs.radacc.neig), sizeof(*(proj))); |
2086 | |
2087 | /* loop over radacc vectors */ |
2088 | for (i = 0; i < edi->vecs.radacc.neig; i++) |
2089 | { |
2090 | /* calculate the projections, radius */ |
2091 | proj[i] = projectx(edi, xcoll, edi->vecs.radacc.vec[i]); |
2092 | rad += pow(proj[i] - edi->vecs.radacc.refproj[i], 2); |
2093 | } |
2094 | rad = sqrt(rad); |
2095 | |
2096 | /* only correct when radius decreased */ |
2097 | if (rad < edi->vecs.radacc.radius) |
2098 | { |
2099 | ratio = edi->vecs.radacc.radius/rad - 1.0; |
2100 | rad = edi->vecs.radacc.radius; |
2101 | } |
2102 | else |
2103 | { |
2104 | edi->vecs.radacc.radius = rad; |
2105 | } |
2106 | |
2107 | /* loop over radacc vectors */ |
2108 | for (i = 0; i < edi->vecs.radacc.neig; i++) |
2109 | { |
2110 | proj[i] -= edi->vecs.radacc.refproj[i]; |
2111 | |
2112 | /* apply the correction */ |
2113 | proj[i] /= edi->sav.sqrtm[i]; |
2114 | proj[i] *= ratio; |
2115 | for (j = 0; j < edi->sav.nr; j++) |
2116 | { |
2117 | svmul(proj[i], edi->vecs.radacc.vec[i][j], vec_dum); |
2118 | rvec_inc(xcoll[j], vec_dum); |
2119 | } |
2120 | } |
2121 | sfree(proj)save_free("proj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2121, (proj)); |
2122 | } |
2123 | |
2124 | |
2125 | struct t_do_radcon { |
2126 | real *proj; |
2127 | }; |
2128 | |
2129 | static void do_radcon(rvec *xcoll, t_edpar *edi) |
2130 | { |
2131 | int i, j; |
2132 | real rad = 0.0, ratio = 0.0; |
2133 | struct t_do_radcon *loc; |
2134 | gmx_bool bFirst; |
2135 | rvec vec_dum; |
2136 | |
2137 | |
2138 | if (edi->buf->do_radcon != NULL((void*)0)) |
2139 | { |
2140 | bFirst = FALSE0; |
2141 | loc = edi->buf->do_radcon; |
2142 | } |
2143 | else |
2144 | { |
2145 | bFirst = TRUE1; |
2146 | snew(edi->buf->do_radcon, 1)(edi->buf->do_radcon) = save_calloc("edi->buf->do_radcon" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2146, (1), sizeof(*(edi->buf->do_radcon))); |
2147 | } |
2148 | loc = edi->buf->do_radcon; |
2149 | |
2150 | if (edi->vecs.radcon.neig == 0) |
2151 | { |
2152 | return; |
2153 | } |
2154 | |
2155 | if (bFirst) |
2156 | { |
2157 | snew(loc->proj, edi->vecs.radcon.neig)(loc->proj) = save_calloc("loc->proj", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2157, (edi->vecs.radcon.neig), sizeof(*(loc->proj))); |
2158 | } |
2159 | |
2160 | /* loop over radcon vectors */ |
2161 | for (i = 0; i < edi->vecs.radcon.neig; i++) |
2162 | { |
2163 | /* calculate the projections, radius */ |
2164 | loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]); |
2165 | rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2); |
2166 | } |
2167 | rad = sqrt(rad); |
2168 | /* only correct when radius increased */ |
2169 | if (rad > edi->vecs.radcon.radius) |
2170 | { |
2171 | ratio = edi->vecs.radcon.radius/rad - 1.0; |
2172 | |
2173 | /* loop over radcon vectors */ |
2174 | for (i = 0; i < edi->vecs.radcon.neig; i++) |
2175 | { |
2176 | /* apply the correction */ |
2177 | loc->proj[i] -= edi->vecs.radcon.refproj[i]; |
2178 | loc->proj[i] /= edi->sav.sqrtm[i]; |
2179 | loc->proj[i] *= ratio; |
2180 | |
2181 | for (j = 0; j < edi->sav.nr; j++) |
2182 | { |
2183 | svmul(loc->proj[i], edi->vecs.radcon.vec[i][j], vec_dum); |
2184 | rvec_inc(xcoll[j], vec_dum); |
2185 | } |
2186 | } |
2187 | } |
2188 | else |
2189 | { |
2190 | edi->vecs.radcon.radius = rad; |
2191 | } |
2192 | |
2193 | if (rad != edi->vecs.radcon.radius) |
2194 | { |
2195 | rad = 0.0; |
2196 | for (i = 0; i < edi->vecs.radcon.neig; i++) |
2197 | { |
2198 | /* calculate the projections, radius */ |
2199 | loc->proj[i] = projectx(edi, xcoll, edi->vecs.radcon.vec[i]); |
2200 | rad += pow(loc->proj[i] - edi->vecs.radcon.refproj[i], 2); |
2201 | } |
2202 | rad = sqrt(rad); |
Value stored to 'rad' is never read | |
2203 | } |
2204 | } |
2205 | |
2206 | |
2207 | static void ed_apply_constraints(rvec *xcoll, t_edpar *edi, gmx_int64_t step) |
2208 | { |
2209 | int i; |
2210 | |
2211 | |
2212 | /* subtract the average positions */ |
2213 | for (i = 0; i < edi->sav.nr; i++) |
2214 | { |
2215 | rvec_dec(xcoll[i], edi->sav.x[i]); |
2216 | } |
2217 | |
2218 | /* apply the constraints */ |
2219 | if (step >= 0) |
2220 | { |
2221 | do_linfix(xcoll, edi, step); |
2222 | } |
2223 | do_linacc(xcoll, edi); |
2224 | if (step >= 0) |
2225 | { |
2226 | do_radfix(xcoll, edi); |
2227 | } |
2228 | do_radacc(xcoll, edi); |
2229 | do_radcon(xcoll, edi); |
2230 | |
2231 | /* add back the average positions */ |
2232 | for (i = 0; i < edi->sav.nr; i++) |
2233 | { |
2234 | rvec_inc(xcoll[i], edi->sav.x[i]); |
2235 | } |
2236 | } |
2237 | |
2238 | |
2239 | /* Write out the projections onto the eigenvectors. The order of output |
2240 | * corresponds to ed_output_legend() */ |
2241 | static void write_edo(t_edpar *edi, FILE *fp, real rmsd) |
2242 | { |
2243 | int i; |
2244 | |
2245 | |
2246 | /* Output how well we fit to the reference structure */ |
2247 | fprintf(fp, EDcol_ffmt"%17f", rmsd); |
2248 | |
2249 | for (i = 0; i < edi->vecs.mon.neig; i++) |
2250 | { |
2251 | fprintf(fp, EDcol_efmt"%17.5e", edi->vecs.mon.xproj[i]); |
2252 | } |
2253 | |
2254 | for (i = 0; i < edi->vecs.linfix.neig; i++) |
2255 | { |
2256 | fprintf(fp, EDcol_efmt"%17.5e", edi->vecs.linfix.xproj[i]); |
2257 | } |
2258 | |
2259 | for (i = 0; i < edi->vecs.linacc.neig; i++) |
2260 | { |
2261 | fprintf(fp, EDcol_efmt"%17.5e", edi->vecs.linacc.xproj[i]); |
2262 | } |
2263 | |
2264 | for (i = 0; i < edi->vecs.radfix.neig; i++) |
2265 | { |
2266 | fprintf(fp, EDcol_efmt"%17.5e", edi->vecs.radfix.xproj[i]); |
2267 | } |
2268 | if (edi->vecs.radfix.neig) |
2269 | { |
2270 | fprintf(fp, EDcol_ffmt"%17f", calc_radius(&edi->vecs.radfix)); /* fixed increment radius */ |
2271 | } |
2272 | |
2273 | for (i = 0; i < edi->vecs.radacc.neig; i++) |
2274 | { |
2275 | fprintf(fp, EDcol_efmt"%17.5e", edi->vecs.radacc.xproj[i]); |
2276 | } |
2277 | if (edi->vecs.radacc.neig) |
2278 | { |
2279 | fprintf(fp, EDcol_ffmt"%17f", calc_radius(&edi->vecs.radacc)); /* acceptance radius */ |
2280 | } |
2281 | |
2282 | for (i = 0; i < edi->vecs.radcon.neig; i++) |
2283 | { |
2284 | fprintf(fp, EDcol_efmt"%17.5e", edi->vecs.radcon.xproj[i]); |
2285 | } |
2286 | if (edi->vecs.radcon.neig) |
2287 | { |
2288 | fprintf(fp, EDcol_ffmt"%17f", calc_radius(&edi->vecs.radcon)); /* contracting radius */ |
2289 | } |
2290 | } |
2291 | |
2292 | /* Returns if any constraints are switched on */ |
2293 | static int ed_constraints(gmx_bool edtype, t_edpar *edi) |
2294 | { |
2295 | if (edtype == eEDedsam || edtype == eEDflood) |
2296 | { |
2297 | return (edi->vecs.linfix.neig || edi->vecs.linacc.neig || |
2298 | edi->vecs.radfix.neig || edi->vecs.radacc.neig || |
2299 | edi->vecs.radcon.neig); |
2300 | } |
2301 | return 0; |
2302 | } |
2303 | |
2304 | |
2305 | /* Copies reference projection 'refproj' to fixed 'refproj0' variable for flooding/ |
2306 | * umbrella sampling simulations. */ |
2307 | static void copyEvecReference(t_eigvec* floodvecs) |
2308 | { |
2309 | int i; |
2310 | |
2311 | |
2312 | if (NULL((void*)0) == floodvecs->refproj0) |
2313 | { |
2314 | snew(floodvecs->refproj0, floodvecs->neig)(floodvecs->refproj0) = save_calloc("floodvecs->refproj0" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2314, (floodvecs->neig), sizeof(*(floodvecs->refproj0 ))); |
2315 | } |
2316 | |
2317 | for (i = 0; i < floodvecs->neig; i++) |
2318 | { |
2319 | floodvecs->refproj0[i] = floodvecs->refproj[i]; |
2320 | } |
2321 | } |
2322 | |
2323 | |
2324 | /* Call on MASTER only. Check whether the essential dynamics / flooding |
2325 | * groups of the checkpoint file are consistent with the provided .edi file. */ |
2326 | static void crosscheck_edi_file_vs_checkpoint(gmx_edsam_t ed, edsamstate_t *EDstate) |
2327 | { |
2328 | t_edpar *edi = NULL((void*)0); /* points to a single edi data set */ |
2329 | int edinum; |
2330 | |
2331 | |
2332 | if (NULL((void*)0) == EDstate->nref || NULL((void*)0) == EDstate->nav) |
2333 | { |
2334 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2334, "Essential dynamics and flooding can only be switched on (or off) at the\n" |
2335 | "start of a new simulation. If a simulation runs with/without ED constraints,\n" |
2336 | "it must also continue with/without ED constraints when checkpointing.\n" |
2337 | "To switch on (or off) ED constraints, please prepare a new .tpr to start\n" |
2338 | "from without a checkpoint.\n"); |
2339 | } |
2340 | |
2341 | edi = ed->edpar; |
2342 | edinum = 0; |
2343 | while (edi != NULL((void*)0)) |
2344 | { |
2345 | /* Check number of atoms in the reference and average structures */ |
2346 | if (EDstate->nref[edinum] != edi->sref.nr) |
2347 | { |
2348 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2348, "The number of reference structure atoms in ED group %c is\n" |
2349 | "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n", |
2350 | get_EDgroupChar(edinum+1, 0), EDstate->nref[edinum], edi->sref.nr); |
2351 | } |
2352 | if (EDstate->nav[edinum] != edi->sav.nr) |
2353 | { |
2354 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2354, "The number of average structure atoms in ED group %c is\n" |
2355 | "not the same in .cpt (NREF=%d) and .edi (NREF=%d) files!\n", |
2356 | get_EDgroupChar(edinum+1, 0), EDstate->nav[edinum], edi->sav.nr); |
2357 | } |
2358 | edi = edi->next_edi; |
2359 | edinum++; |
2360 | } |
2361 | |
2362 | if (edinum != EDstate->nED) |
2363 | { |
2364 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2364, "The number of essential dynamics / flooding groups is not consistent.\n" |
2365 | "There are %d ED groups in the .cpt file, but %d in the .edi file!\n" |
2366 | "Are you sure this is the correct .edi file?\n", EDstate->nED, edinum); |
2367 | } |
2368 | } |
2369 | |
2370 | |
2371 | /* The edsamstate struct stores the information we need to make the ED group |
2372 | * whole again after restarts from a checkpoint file. Here we do the following: |
2373 | * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing, |
2374 | * b) if we did start from .cpt, we copy over the last whole structures from .cpt, |
2375 | * c) in any case, for subsequent checkpoint writing, we set the pointers in |
2376 | * edsamstate to the x_old arrays, which contain the correct PBC representation of |
2377 | * all ED structures at the last time step. */ |
2378 | static void init_edsamstate(gmx_edsam_t ed, edsamstate_t *EDstate) |
2379 | { |
2380 | int i, nr_edi; |
2381 | t_edpar *edi; |
2382 | |
2383 | |
2384 | snew(EDstate->old_sref_p, EDstate->nED)(EDstate->old_sref_p) = save_calloc("EDstate->old_sref_p" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2384, (EDstate->nED), sizeof(*(EDstate->old_sref_p))); |
2385 | snew(EDstate->old_sav_p, EDstate->nED)(EDstate->old_sav_p) = save_calloc("EDstate->old_sav_p" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2385, (EDstate->nED), sizeof(*(EDstate->old_sav_p))); |
2386 | |
2387 | /* If we did not read in a .cpt file, these arrays are not yet allocated */ |
2388 | if (!EDstate->bFromCpt) |
2389 | { |
2390 | snew(EDstate->nref, EDstate->nED)(EDstate->nref) = save_calloc("EDstate->nref", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2390, (EDstate->nED), sizeof(*(EDstate->nref))); |
2391 | snew(EDstate->nav, EDstate->nED)(EDstate->nav) = save_calloc("EDstate->nav", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2391, (EDstate->nED), sizeof(*(EDstate->nav))); |
2392 | } |
2393 | |
2394 | /* Loop over all ED/flooding data sets (usually only one, though) */ |
2395 | edi = ed->edpar; |
2396 | for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++) |
2397 | { |
2398 | /* We always need the last reference and average positions such that |
2399 | * in the next time step we can make the ED group whole again |
2400 | * if the atoms do not have the correct PBC representation */ |
2401 | if (EDstate->bFromCpt) |
2402 | { |
2403 | /* Copy the last whole positions of reference and average group from .cpt */ |
2404 | for (i = 0; i < edi->sref.nr; i++) |
2405 | { |
2406 | copy_rvec(EDstate->old_sref[nr_edi-1][i], edi->sref.x_old[i]); |
2407 | } |
2408 | for (i = 0; i < edi->sav.nr; i++) |
2409 | { |
2410 | copy_rvec(EDstate->old_sav [nr_edi-1][i], edi->sav.x_old [i]); |
2411 | } |
2412 | } |
2413 | else |
2414 | { |
2415 | EDstate->nref[nr_edi-1] = edi->sref.nr; |
2416 | EDstate->nav [nr_edi-1] = edi->sav.nr; |
2417 | } |
2418 | |
2419 | /* For subsequent checkpoint writing, set the edsamstate pointers to the edi arrays: */ |
2420 | EDstate->old_sref_p[nr_edi-1] = edi->sref.x_old; |
2421 | EDstate->old_sav_p [nr_edi-1] = edi->sav.x_old; |
2422 | |
2423 | edi = edi->next_edi; |
2424 | } |
2425 | } |
2426 | |
2427 | |
2428 | /* Adds 'buf' to 'str' */ |
2429 | static void add_to_string(char **str, char *buf) |
2430 | { |
2431 | int len; |
2432 | |
2433 | |
2434 | len = strlen(*str) + strlen(buf) + 1; |
2435 | srenew(*str, len)(*str) = save_realloc("*str", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2435, (*str), (len), sizeof(*(*str))); |
2436 | strcat(*str, buf); |
2437 | } |
2438 | |
2439 | |
2440 | static void add_to_string_aligned(char **str, char *buf) |
2441 | { |
2442 | char buf_aligned[STRLEN4096]; |
2443 | |
2444 | sprintf(buf_aligned, EDcol_sfmt"%17s", buf); |
2445 | add_to_string(str, buf_aligned); |
2446 | } |
2447 | |
2448 | |
2449 | static void nice_legend(const char ***setname, int *nsets, char **LegendStr, char *value, char *unit, char EDgroupchar) |
2450 | { |
2451 | char tmp[STRLEN4096], tmp2[STRLEN4096]; |
2452 | |
2453 | |
2454 | sprintf(tmp, "%c %s", EDgroupchar, value); |
2455 | add_to_string_aligned(LegendStr, tmp); |
2456 | sprintf(tmp2, "%s (%s)", tmp, unit); |
2457 | (*setname)[*nsets] = strdup(tmp2)(__extension__ (__builtin_constant_p (tmp2) && ((size_t )(const void *)((tmp2) + 1) - (size_t)(const void *)(tmp2) == 1) ? (((const char *) (tmp2))[0] == '\0' ? (char *) calloc ( (size_t) 1, (size_t) 1) : ({ size_t __len = strlen (tmp2) + 1 ; char *__retval = (char *) malloc (__len); if (__retval != ( (void*)0)) __retval = (char *) memcpy (__retval, tmp2, __len) ; __retval; })) : __strdup (tmp2))); |
2458 | (*nsets)++; |
2459 | } |
2460 | |
2461 | |
2462 | static void nice_legend_evec(const char ***setname, int *nsets, char **LegendStr, t_eigvec *evec, char EDgroupChar, const char *EDtype) |
2463 | { |
2464 | int i; |
2465 | char tmp[STRLEN4096]; |
2466 | |
2467 | |
2468 | for (i = 0; i < evec->neig; i++) |
2469 | { |
2470 | sprintf(tmp, "EV%dprj%s", evec->ieig[i], EDtype); |
2471 | nice_legend(setname, nsets, LegendStr, tmp, "nm", EDgroupChar); |
2472 | } |
2473 | } |
2474 | |
2475 | |
2476 | /* Makes a legend for the xvg output file. Call on MASTER only! */ |
2477 | static void write_edo_legend(gmx_edsam_t ed, int nED, const output_env_t oenv) |
2478 | { |
2479 | t_edpar *edi = NULL((void*)0); |
2480 | int i; |
2481 | int nr_edi, nsets, n_flood, n_edsam; |
2482 | const char **setname; |
2483 | char buf[STRLEN4096]; |
2484 | char *LegendStr = NULL((void*)0); |
2485 | |
2486 | |
2487 | edi = ed->edpar; |
2488 | |
2489 | fprintf(ed->edo, "# Output will be written every %d step%s\n", ed->edpar->outfrq, ed->edpar->outfrq != 1 ? "s" : ""); |
2490 | |
2491 | for (nr_edi = 1; nr_edi <= nED; nr_edi++) |
2492 | { |
2493 | fprintf(ed->edo, "#\n"); |
2494 | fprintf(ed->edo, "# Summary of applied con/restraints for the ED group %c\n", get_EDgroupChar(nr_edi, nED)); |
2495 | fprintf(ed->edo, "# Atoms in average structure: %d\n", edi->sav.nr); |
2496 | fprintf(ed->edo, "# monitor : %d vec%s\n", edi->vecs.mon.neig, edi->vecs.mon.neig != 1 ? "s" : ""); |
2497 | fprintf(ed->edo, "# LINFIX : %d vec%s\n", edi->vecs.linfix.neig, edi->vecs.linfix.neig != 1 ? "s" : ""); |
2498 | fprintf(ed->edo, "# LINACC : %d vec%s\n", edi->vecs.linacc.neig, edi->vecs.linacc.neig != 1 ? "s" : ""); |
2499 | fprintf(ed->edo, "# RADFIX : %d vec%s\n", edi->vecs.radfix.neig, edi->vecs.radfix.neig != 1 ? "s" : ""); |
2500 | fprintf(ed->edo, "# RADACC : %d vec%s\n", edi->vecs.radacc.neig, edi->vecs.radacc.neig != 1 ? "s" : ""); |
2501 | fprintf(ed->edo, "# RADCON : %d vec%s\n", edi->vecs.radcon.neig, edi->vecs.radcon.neig != 1 ? "s" : ""); |
2502 | fprintf(ed->edo, "# FLOODING : %d vec%s ", edi->flood.vecs.neig, edi->flood.vecs.neig != 1 ? "s" : ""); |
2503 | |
2504 | if (edi->flood.vecs.neig) |
2505 | { |
2506 | /* If in any of the groups we find a flooding vector, flooding is turned on */ |
2507 | ed->eEDtype = eEDflood; |
2508 | |
2509 | /* Print what flavor of flooding we will do */ |
2510 | if (0 == edi->flood.tau) /* constant flooding strength */ |
2511 | { |
2512 | fprintf(ed->edo, "Efl_null = %g", edi->flood.constEfl); |
2513 | if (edi->flood.bHarmonic) |
2514 | { |
2515 | fprintf(ed->edo, ", harmonic"); |
2516 | } |
2517 | } |
2518 | else /* adaptive flooding */ |
2519 | { |
2520 | fprintf(ed->edo, ", adaptive"); |
2521 | } |
2522 | } |
2523 | fprintf(ed->edo, "\n"); |
2524 | |
2525 | edi = edi->next_edi; |
2526 | } |
2527 | |
2528 | /* Print a nice legend */ |
2529 | snew(LegendStr, 1)(LegendStr) = save_calloc("LegendStr", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2529, (1), sizeof(*(LegendStr))); |
2530 | LegendStr[0] = '\0'; |
2531 | sprintf(buf, "# %6s", "time"); |
2532 | add_to_string(&LegendStr, buf); |
2533 | |
2534 | /* Calculate the maximum number of columns we could end up with */ |
2535 | edi = ed->edpar; |
2536 | nsets = 0; |
2537 | for (nr_edi = 1; nr_edi <= nED; nr_edi++) |
2538 | { |
2539 | nsets += 5 +edi->vecs.mon.neig |
2540 | +edi->vecs.linfix.neig |
2541 | +edi->vecs.linacc.neig |
2542 | +edi->vecs.radfix.neig |
2543 | +edi->vecs.radacc.neig |
2544 | +edi->vecs.radcon.neig |
2545 | + 6*edi->flood.vecs.neig; |
2546 | edi = edi->next_edi; |
2547 | } |
2548 | snew(setname, nsets)(setname) = save_calloc("setname", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2548, (nsets), sizeof(*(setname))); |
2549 | |
2550 | /* In the mdrun time step in a first function call (do_flood()) the flooding |
2551 | * forces are calculated and in a second function call (do_edsam()) the |
2552 | * ED constraints. To get a corresponding legend, we need to loop twice |
2553 | * over the edi groups and output first the flooding, then the ED part */ |
2554 | |
2555 | /* The flooding-related legend entries, if flooding is done */ |
2556 | nsets = 0; |
2557 | if (eEDflood == ed->eEDtype) |
2558 | { |
2559 | edi = ed->edpar; |
2560 | for (nr_edi = 1; nr_edi <= nED; nr_edi++) |
2561 | { |
2562 | /* Always write out the projection on the flooding EVs. Of course, this can also |
2563 | * be achieved with the monitoring option in do_edsam() (if switched on by the |
2564 | * user), but in that case the positions need to be communicated in do_edsam(), |
2565 | * which is not necessary when doing flooding only. */ |
2566 | nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) ); |
2567 | |
2568 | for (i = 0; i < edi->flood.vecs.neig; i++) |
2569 | { |
2570 | sprintf(buf, "EV%dprjFLOOD", edi->flood.vecs.ieig[i]); |
2571 | nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED)); |
2572 | |
2573 | /* Output the current reference projection if it changes with time; |
2574 | * this can happen when flooding is used as harmonic restraint */ |
2575 | if (edi->flood.bHarmonic && edi->flood.vecs.refprojslope[i] != 0.0) |
2576 | { |
2577 | sprintf(buf, "EV%d ref.prj.", edi->flood.vecs.ieig[i]); |
2578 | nice_legend(&setname, &nsets, &LegendStr, buf, "nm", get_EDgroupChar(nr_edi, nED)); |
2579 | } |
2580 | |
2581 | /* For flooding we also output Efl, Vfl, deltaF, and the flooding forces */ |
2582 | if (0 != edi->flood.tau) /* only output Efl for adaptive flooding (constant otherwise) */ |
2583 | { |
2584 | sprintf(buf, "EV%d-Efl", edi->flood.vecs.ieig[i]); |
2585 | nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED)); |
2586 | } |
2587 | |
2588 | sprintf(buf, "EV%d-Vfl", edi->flood.vecs.ieig[i]); |
2589 | nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED)); |
2590 | |
2591 | if (0 != edi->flood.tau) /* only output deltaF for adaptive flooding (zero otherwise) */ |
2592 | { |
2593 | sprintf(buf, "EV%d-deltaF", edi->flood.vecs.ieig[i]); |
2594 | nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol", get_EDgroupChar(nr_edi, nED)); |
2595 | } |
2596 | |
2597 | sprintf(buf, "EV%d-FLforces", edi->flood.vecs.ieig[i]); |
2598 | nice_legend(&setname, &nsets, &LegendStr, buf, "kJ/mol/nm", get_EDgroupChar(nr_edi, nED)); |
2599 | } |
2600 | |
2601 | edi = edi->next_edi; |
2602 | } /* End of flooding-related legend entries */ |
2603 | } |
2604 | n_flood = nsets; |
2605 | |
2606 | /* Now the ED-related entries, if essential dynamics is done */ |
2607 | edi = ed->edpar; |
2608 | for (nr_edi = 1; nr_edi <= nED; nr_edi++) |
2609 | { |
2610 | if (bNeedDoEdsam(edi)) /* Only print ED legend if at least one ED option is on */ |
2611 | { |
2612 | nice_legend(&setname, &nsets, &LegendStr, "RMSD to ref", "nm", get_EDgroupChar(nr_edi, nED) ); |
2613 | |
2614 | /* Essential dynamics, projections on eigenvectors */ |
2615 | nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.mon, get_EDgroupChar(nr_edi, nED), "MON" ); |
2616 | nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linfix, get_EDgroupChar(nr_edi, nED), "LINFIX"); |
2617 | nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.linacc, get_EDgroupChar(nr_edi, nED), "LINACC"); |
2618 | nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radfix, get_EDgroupChar(nr_edi, nED), "RADFIX"); |
2619 | if (edi->vecs.radfix.neig) |
2620 | { |
2621 | nice_legend(&setname, &nsets, &LegendStr, "RADFIX radius", "nm", get_EDgroupChar(nr_edi, nED)); |
2622 | } |
2623 | nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radacc, get_EDgroupChar(nr_edi, nED), "RADACC"); |
2624 | if (edi->vecs.radacc.neig) |
2625 | { |
2626 | nice_legend(&setname, &nsets, &LegendStr, "RADACC radius", "nm", get_EDgroupChar(nr_edi, nED)); |
2627 | } |
2628 | nice_legend_evec(&setname, &nsets, &LegendStr, &edi->vecs.radcon, get_EDgroupChar(nr_edi, nED), "RADCON"); |
2629 | if (edi->vecs.radcon.neig) |
2630 | { |
2631 | nice_legend(&setname, &nsets, &LegendStr, "RADCON radius", "nm", get_EDgroupChar(nr_edi, nED)); |
2632 | } |
2633 | } |
2634 | edi = edi->next_edi; |
2635 | } /* end of 'pure' essential dynamics legend entries */ |
2636 | n_edsam = nsets - n_flood; |
2637 | |
2638 | xvgr_legend(ed->edo, nsets, setname, oenv); |
2639 | sfree(setname)save_free("setname", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2639, (setname)); |
2640 | |
2641 | fprintf(ed->edo, "#\n" |
2642 | "# Legend for %d column%s of flooding plus %d column%s of essential dynamics data:\n", |
2643 | n_flood, 1 == n_flood ? "" : "s", |
2644 | n_edsam, 1 == n_edsam ? "" : "s"); |
2645 | fprintf(ed->edo, "%s", LegendStr); |
2646 | sfree(LegendStr)save_free("LegendStr", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2646, (LegendStr)); |
2647 | |
2648 | fflush(ed->edo); |
2649 | } |
2650 | |
2651 | |
2652 | /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle |
2653 | * contained in the input file, creates a NULL terminated list of t_edpar structures */ |
2654 | void init_edsam(gmx_mtop_t *mtop, |
2655 | t_inputrec *ir, |
2656 | t_commrec *cr, |
2657 | gmx_edsam_t ed, |
2658 | rvec x[], |
2659 | matrix box, |
2660 | edsamstate_t *EDstate) |
2661 | { |
2662 | t_edpar *edi = NULL((void*)0); /* points to a single edi data set */ |
2663 | int i, nr_edi, avindex; |
2664 | rvec *x_pbc = NULL((void*)0); /* positions of the whole MD system with pbc removed */ |
2665 | rvec *xfit = NULL((void*)0), *xstart = NULL((void*)0); /* dummy arrays to determine initial RMSDs */ |
2666 | rvec fit_transvec; /* translation ... */ |
2667 | matrix fit_rotmat; /* ... and rotation from fit to reference structure */ |
2668 | rvec *ref_x_old = NULL((void*)0); /* helper pointer */ |
2669 | |
2670 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
2671 | { |
2672 | fprintf(stderrstderr, "ED: Initializing essential dynamics constraints.\n"); |
2673 | |
2674 | if (NULL((void*)0) == ed) |
2675 | { |
2676 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2676, "The checkpoint file you provided is from an essential dynamics or\n" |
2677 | "flooding simulation. Please also provide the correct .edi file with -ei.\n"); |
2678 | } |
2679 | } |
2680 | |
2681 | /* Needed for initializing radacc radius in do_edsam */ |
2682 | ed->bFirst = TRUE1; |
2683 | |
2684 | /* The input file is read by the master and the edi structures are |
2685 | * initialized here. Input is stored in ed->edpar. Then the edi |
2686 | * structures are transferred to the other nodes */ |
2687 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
2688 | { |
2689 | /* Initialization for every ED/flooding group. Flooding uses one edi group per |
2690 | * flooding vector, Essential dynamics can be applied to more than one structure |
2691 | * as well, but will be done in the order given in the edi file, so |
2692 | * expect different results for different order of edi file concatenation! */ |
2693 | edi = ed->edpar; |
2694 | while (edi != NULL((void*)0)) |
2695 | { |
2696 | init_edi(mtop, edi); |
2697 | init_flood(edi, ed, ir->delta_t); |
2698 | edi = edi->next_edi; |
2699 | } |
2700 | } |
2701 | |
2702 | /* The master does the work here. The other nodes get the positions |
2703 | * not before dd_partition_system which is called after init_edsam */ |
2704 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
2705 | { |
2706 | if (!EDstate->bFromCpt) |
2707 | { |
2708 | /* Remove PBC, make molecule(s) subject to ED whole. */ |
2709 | snew(x_pbc, mtop->natoms)(x_pbc) = save_calloc("x_pbc", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2709, (mtop->natoms), sizeof(*(x_pbc))); |
2710 | m_rveccopy(mtop->natoms, x, x_pbc); |
2711 | do_pbc_first_mtop(NULL((void*)0), ir->ePBC, box, mtop, x_pbc); |
2712 | } |
2713 | /* Reset pointer to first ED data set which contains the actual ED data */ |
2714 | edi = ed->edpar; |
2715 | /* Loop over all ED/flooding data sets (usually only one, though) */ |
2716 | for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++) |
2717 | { |
2718 | /* For multiple ED groups we use the output frequency that was specified |
2719 | * in the first set */ |
2720 | if (nr_edi > 1) |
2721 | { |
2722 | edi->outfrq = ed->edpar->outfrq; |
2723 | } |
2724 | |
2725 | /* Extract the initial reference and average positions. When starting |
2726 | * from .cpt, these have already been read into sref.x_old |
2727 | * in init_edsamstate() */ |
2728 | if (!EDstate->bFromCpt) |
2729 | { |
2730 | /* If this is the first run (i.e. no checkpoint present) we assume |
2731 | * that the starting positions give us the correct PBC representation */ |
2732 | for (i = 0; i < edi->sref.nr; i++) |
2733 | { |
2734 | copy_rvec(x_pbc[edi->sref.anrs[i]], edi->sref.x_old[i]); |
2735 | } |
2736 | |
2737 | for (i = 0; i < edi->sav.nr; i++) |
2738 | { |
2739 | copy_rvec(x_pbc[edi->sav.anrs[i]], edi->sav.x_old[i]); |
2740 | } |
2741 | } |
2742 | |
2743 | /* Now we have the PBC-correct start positions of the reference and |
2744 | average structure. We copy that over to dummy arrays on which we |
2745 | can apply fitting to print out the RMSD. We srenew the memory since |
2746 | the size of the buffers is likely different for every ED group */ |
2747 | srenew(xfit, edi->sref.nr )(xfit) = save_realloc("xfit", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2747, (xfit), (edi->sref.nr), sizeof(*(xfit))); |
2748 | srenew(xstart, edi->sav.nr )(xstart) = save_realloc("xstart", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2748, (xstart), (edi->sav.nr), sizeof(*(xstart))); |
2749 | if (edi->bRefEqAv) |
2750 | { |
2751 | /* Reference indices are the same as average indices */ |
2752 | ref_x_old = edi->sav.x_old; |
2753 | } |
2754 | else |
2755 | { |
2756 | ref_x_old = edi->sref.x_old; |
2757 | } |
2758 | copy_rvecn(ref_x_old, xfit, 0, edi->sref.nr); |
2759 | copy_rvecn(edi->sav.x_old, xstart, 0, edi->sav.nr); |
2760 | |
2761 | /* Make the fit to the REFERENCE structure, get translation and rotation */ |
2762 | fit_to_reference(xfit, fit_transvec, fit_rotmat, edi); |
2763 | |
2764 | /* Output how well we fit to the reference at the start */ |
2765 | translate_and_rotate(xfit, edi->sref.nr, fit_transvec, fit_rotmat); |
2766 | fprintf(stderrstderr, "ED: Initial RMSD from reference after fit = %f nm", |
2767 | rmsd_from_structure(xfit, &edi->sref)); |
2768 | if (EDstate->nED > 1) |
2769 | { |
2770 | fprintf(stderrstderr, " (ED group %c)", get_EDgroupChar(nr_edi, EDstate->nED)); |
2771 | } |
2772 | fprintf(stderrstderr, "\n"); |
2773 | |
2774 | /* Now apply the translation and rotation to the atoms on which ED sampling will be performed */ |
2775 | translate_and_rotate(xstart, edi->sav.nr, fit_transvec, fit_rotmat); |
2776 | |
2777 | /* calculate initial projections */ |
2778 | project(xstart, edi); |
2779 | |
2780 | /* For the target and origin structure both a reference (fit) and an |
2781 | * average structure can be provided in make_edi. If both structures |
2782 | * are the same, make_edi only stores one of them in the .edi file. |
2783 | * If they differ, first the fit and then the average structure is stored |
2784 | * in star (or sor), thus the number of entries in star/sor is |
2785 | * (n_fit + n_av) with n_fit the size of the fitting group and n_av |
2786 | * the size of the average group. */ |
2787 | |
2788 | /* process target structure, if required */ |
2789 | if (edi->star.nr > 0) |
2790 | { |
2791 | fprintf(stderrstderr, "ED: Fitting target structure to reference structure\n"); |
2792 | |
2793 | /* get translation & rotation for fit of target structure to reference structure */ |
2794 | fit_to_reference(edi->star.x, fit_transvec, fit_rotmat, edi); |
2795 | /* do the fit */ |
2796 | translate_and_rotate(edi->star.x, edi->star.nr, fit_transvec, fit_rotmat); |
2797 | if (edi->star.nr == edi->sav.nr) |
2798 | { |
2799 | avindex = 0; |
2800 | } |
2801 | else /* edi->star.nr = edi->sref.nr + edi->sav.nr */ |
2802 | { |
2803 | /* The last sav.nr indices of the target structure correspond to |
2804 | * the average structure, which must be projected */ |
2805 | avindex = edi->star.nr - edi->sav.nr; |
2806 | } |
2807 | rad_project(edi, &edi->star.x[avindex], &edi->vecs.radcon); |
2808 | } |
2809 | else |
2810 | { |
2811 | rad_project(edi, xstart, &edi->vecs.radcon); |
2812 | } |
2813 | |
2814 | /* process structure that will serve as origin of expansion circle */ |
2815 | if ( (eEDflood == ed->eEDtype) && (FALSE0 == edi->flood.bConstForce) ) |
2816 | { |
2817 | fprintf(stderrstderr, "ED: Setting center of flooding potential (0 = average structure)\n"); |
2818 | } |
2819 | |
2820 | if (edi->sori.nr > 0) |
2821 | { |
2822 | fprintf(stderrstderr, "ED: Fitting origin structure to reference structure\n"); |
2823 | |
2824 | /* fit this structure to reference structure */ |
2825 | fit_to_reference(edi->sori.x, fit_transvec, fit_rotmat, edi); |
2826 | /* do the fit */ |
2827 | translate_and_rotate(edi->sori.x, edi->sori.nr, fit_transvec, fit_rotmat); |
2828 | if (edi->sori.nr == edi->sav.nr) |
2829 | { |
2830 | avindex = 0; |
2831 | } |
2832 | else /* edi->sori.nr = edi->sref.nr + edi->sav.nr */ |
2833 | { |
2834 | /* For the projection, we need the last sav.nr indices of sori */ |
2835 | avindex = edi->sori.nr - edi->sav.nr; |
2836 | } |
2837 | |
2838 | rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radacc); |
2839 | rad_project(edi, &edi->sori.x[avindex], &edi->vecs.radfix); |
2840 | if ( (eEDflood == ed->eEDtype) && (FALSE0 == edi->flood.bConstForce) ) |
2841 | { |
2842 | fprintf(stderrstderr, "ED: The ORIGIN structure will define the flooding potential center.\n"); |
2843 | /* Set center of flooding potential to the ORIGIN structure */ |
2844 | rad_project(edi, &edi->sori.x[avindex], &edi->flood.vecs); |
2845 | /* We already know that no (moving) reference position was provided, |
2846 | * therefore we can overwrite refproj[0]*/ |
2847 | copyEvecReference(&edi->flood.vecs); |
2848 | } |
2849 | } |
2850 | else /* No origin structure given */ |
2851 | { |
2852 | rad_project(edi, xstart, &edi->vecs.radacc); |
2853 | rad_project(edi, xstart, &edi->vecs.radfix); |
2854 | if ( (eEDflood == ed->eEDtype) && (FALSE0 == edi->flood.bConstForce) ) |
2855 | { |
2856 | if (edi->flood.bHarmonic) |
2857 | { |
2858 | fprintf(stderrstderr, "ED: A (possibly changing) ref. projection will define the flooding potential center.\n"); |
2859 | for (i = 0; i < edi->flood.vecs.neig; i++) |
2860 | { |
2861 | edi->flood.vecs.refproj[i] = edi->flood.vecs.refproj0[i]; |
2862 | } |
2863 | } |
2864 | else |
2865 | { |
2866 | fprintf(stderrstderr, "ED: The AVERAGE structure will define the flooding potential center.\n"); |
2867 | /* Set center of flooding potential to the center of the covariance matrix, |
2868 | * i.e. the average structure, i.e. zero in the projected system */ |
2869 | for (i = 0; i < edi->flood.vecs.neig; i++) |
2870 | { |
2871 | edi->flood.vecs.refproj[i] = 0.0; |
2872 | } |
2873 | } |
2874 | } |
2875 | } |
2876 | /* For convenience, output the center of the flooding potential for the eigenvectors */ |
2877 | if ( (eEDflood == ed->eEDtype) && (FALSE0 == edi->flood.bConstForce) ) |
2878 | { |
2879 | for (i = 0; i < edi->flood.vecs.neig; i++) |
2880 | { |
2881 | fprintf(stdoutstdout, "ED: EV %d flooding potential center: %11.4e", edi->flood.vecs.ieig[i], edi->flood.vecs.refproj[i]); |
2882 | if (edi->flood.bHarmonic) |
2883 | { |
2884 | fprintf(stdoutstdout, " (adding %11.4e/timestep)", edi->flood.vecs.refprojslope[i]); |
2885 | } |
2886 | fprintf(stdoutstdout, "\n"); |
2887 | } |
2888 | } |
2889 | |
2890 | /* set starting projections for linsam */ |
2891 | rad_project(edi, xstart, &edi->vecs.linacc); |
2892 | rad_project(edi, xstart, &edi->vecs.linfix); |
2893 | |
2894 | /* Prepare for the next edi data set: */ |
2895 | edi = edi->next_edi; |
2896 | } |
2897 | /* Cleaning up on the master node: */ |
2898 | if (!EDstate->bFromCpt) |
2899 | { |
2900 | sfree(x_pbc)save_free("x_pbc", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2900, (x_pbc)); |
2901 | } |
2902 | sfree(xfit)save_free("xfit", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2902, (xfit)); |
2903 | sfree(xstart)save_free("xstart", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2903, (xstart)); |
2904 | |
2905 | } /* end of MASTER only section */ |
2906 | |
2907 | if (PAR(cr)((cr)->nnodes > 1)) |
2908 | { |
2909 | /* First let everybody know how many ED data sets to expect */ |
2910 | gmx_bcast(sizeof(EDstate->nED), &EDstate->nED, cr); |
2911 | /* Broadcast the essential dynamics / flooding data to all nodes */ |
2912 | broadcast_ed_data(cr, ed, EDstate->nED); |
2913 | } |
2914 | else |
2915 | { |
2916 | /* In the single-CPU case, point the local atom numbers pointers to the global |
2917 | * one, so that we can use the same notation in serial and parallel case: */ |
2918 | |
2919 | /* Loop over all ED data sets (usually only one, though) */ |
2920 | edi = ed->edpar; |
2921 | for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++) |
2922 | { |
2923 | edi->sref.anrs_loc = edi->sref.anrs; |
2924 | edi->sav.anrs_loc = edi->sav.anrs; |
2925 | edi->star.anrs_loc = edi->star.anrs; |
2926 | edi->sori.anrs_loc = edi->sori.anrs; |
2927 | /* For the same reason as above, make a dummy c_ind array: */ |
2928 | snew(edi->sav.c_ind, edi->sav.nr)(edi->sav.c_ind) = save_calloc("edi->sav.c_ind", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2928, (edi->sav.nr), sizeof(*(edi->sav.c_ind))); |
2929 | /* Initialize the array */ |
2930 | for (i = 0; i < edi->sav.nr; i++) |
2931 | { |
2932 | edi->sav.c_ind[i] = i; |
2933 | } |
2934 | /* In the general case we will need a different-sized array for the reference indices: */ |
2935 | if (!edi->bRefEqAv) |
2936 | { |
2937 | snew(edi->sref.c_ind, edi->sref.nr)(edi->sref.c_ind) = save_calloc("edi->sref.c_ind", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2937, (edi->sref.nr), sizeof(*(edi->sref.c_ind))); |
2938 | for (i = 0; i < edi->sref.nr; i++) |
2939 | { |
2940 | edi->sref.c_ind[i] = i; |
2941 | } |
2942 | } |
2943 | /* Point to the very same array in case of other structures: */ |
2944 | edi->star.c_ind = edi->sav.c_ind; |
2945 | edi->sori.c_ind = edi->sav.c_ind; |
2946 | /* In the serial case, the local number of atoms is the global one: */ |
2947 | edi->sref.nr_loc = edi->sref.nr; |
2948 | edi->sav.nr_loc = edi->sav.nr; |
2949 | edi->star.nr_loc = edi->star.nr; |
2950 | edi->sori.nr_loc = edi->sori.nr; |
2951 | |
2952 | /* An on we go to the next ED group */ |
2953 | edi = edi->next_edi; |
2954 | } |
2955 | } |
2956 | |
2957 | /* Allocate space for ED buffer variables */ |
2958 | /* Again, loop over ED data sets */ |
2959 | edi = ed->edpar; |
2960 | for (nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++) |
2961 | { |
2962 | /* Allocate space for ED buffer variables */ |
2963 | snew_bc(cr, edi->buf, 1){ if (!(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) {((edi->buf)) = save_calloc("(edi->buf)", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2963, ((1)), sizeof(*((edi->buf)))); }}; /* MASTER has already allocated edi->buf in init_edi() */ |
2964 | snew(edi->buf->do_edsam, 1)(edi->buf->do_edsam) = save_calloc("edi->buf->do_edsam" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2964, (1), sizeof(*(edi->buf->do_edsam))); |
2965 | |
2966 | /* Space for collective ED buffer variables */ |
2967 | |
2968 | /* Collective positions of atoms with the average indices */ |
2969 | snew(edi->buf->do_edsam->xcoll, edi->sav.nr)(edi->buf->do_edsam->xcoll) = save_calloc("edi->buf->do_edsam->xcoll" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2969, (edi->sav.nr), sizeof(*(edi->buf->do_edsam-> xcoll))); |
2970 | snew(edi->buf->do_edsam->shifts_xcoll, edi->sav.nr)(edi->buf->do_edsam->shifts_xcoll) = save_calloc("edi->buf->do_edsam->shifts_xcoll" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2970, (edi->sav.nr), sizeof(*(edi->buf->do_edsam-> shifts_xcoll))); /* buffer for xcoll shifts */ |
2971 | snew(edi->buf->do_edsam->extra_shifts_xcoll, edi->sav.nr)(edi->buf->do_edsam->extra_shifts_xcoll) = save_calloc ("edi->buf->do_edsam->extra_shifts_xcoll", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2971, (edi->sav.nr), sizeof(*(edi->buf->do_edsam-> extra_shifts_xcoll))); |
2972 | /* Collective positions of atoms with the reference indices */ |
2973 | if (!edi->bRefEqAv) |
2974 | { |
2975 | snew(edi->buf->do_edsam->xc_ref, edi->sref.nr)(edi->buf->do_edsam->xc_ref) = save_calloc("edi->buf->do_edsam->xc_ref" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2975, (edi->sref.nr), sizeof(*(edi->buf->do_edsam-> xc_ref))); |
2976 | snew(edi->buf->do_edsam->shifts_xc_ref, edi->sref.nr)(edi->buf->do_edsam->shifts_xc_ref) = save_calloc("edi->buf->do_edsam->shifts_xc_ref" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2976, (edi->sref.nr), sizeof(*(edi->buf->do_edsam-> shifts_xc_ref))); /* To store the shifts in */ |
2977 | snew(edi->buf->do_edsam->extra_shifts_xc_ref, edi->sref.nr)(edi->buf->do_edsam->extra_shifts_xc_ref) = save_calloc ("edi->buf->do_edsam->extra_shifts_xc_ref", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2977, (edi->sref.nr), sizeof(*(edi->buf->do_edsam-> extra_shifts_xc_ref))); |
2978 | } |
2979 | |
2980 | /* Get memory for flooding forces */ |
2981 | snew(edi->flood.forces_cartesian, edi->sav.nr)(edi->flood.forces_cartesian) = save_calloc("edi->flood.forces_cartesian" , "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 2981, (edi->sav.nr), sizeof(*(edi->flood.forces_cartesian ))); |
2982 | |
2983 | #ifdef DUMPEDI |
2984 | /* Dump it all into one file per process */ |
2985 | dump_edi(edi, cr, nr_edi); |
2986 | #endif |
2987 | |
2988 | /* Next ED group */ |
2989 | edi = edi->next_edi; |
2990 | } |
2991 | |
2992 | /* Flush the edo file so that the user can check some things |
2993 | * when the simulation has started */ |
2994 | if (ed->edo) |
2995 | { |
2996 | fflush(ed->edo); |
2997 | } |
2998 | } |
2999 | |
3000 | |
3001 | void do_edsam(t_inputrec *ir, |
3002 | gmx_int64_t step, |
3003 | t_commrec *cr, |
3004 | rvec xs[], |
3005 | rvec v[], |
3006 | matrix box, |
3007 | gmx_edsam_t ed) |
3008 | { |
3009 | int i, edinr, iupdate = 500; |
3010 | matrix rotmat; /* rotation matrix */ |
3011 | rvec transvec; /* translation vector */ |
3012 | rvec dv, dx, x_unsh; /* tmp vectors for velocity, distance, unshifted x coordinate */ |
3013 | real dt_1; /* 1/dt */ |
3014 | struct t_do_edsam *buf; |
3015 | t_edpar *edi; |
3016 | real rmsdev = -1; /* RMSD from reference structure prior to applying the constraints */ |
3017 | gmx_bool bSuppress = FALSE0; /* Write .xvg output file on master? */ |
3018 | |
3019 | |
3020 | /* Check if ED sampling has to be performed */ |
3021 | if (ed->eEDtype == eEDnone) |
3022 | { |
3023 | return; |
3024 | } |
3025 | |
3026 | /* Suppress output on first call of do_edsam if |
3027 | * two-step sd2 integrator is used */ |
3028 | if ( (ir->eI == eiSD2) && (v != NULL((void*)0)) ) |
3029 | { |
3030 | bSuppress = TRUE1; |
3031 | } |
3032 | |
3033 | dt_1 = 1.0/ir->delta_t; |
3034 | |
3035 | /* Loop over all ED groups (usually one) */ |
3036 | edi = ed->edpar; |
3037 | edinr = 0; |
3038 | while (edi != NULL((void*)0)) |
3039 | { |
3040 | edinr++; |
3041 | if (bNeedDoEdsam(edi)) |
3042 | { |
3043 | |
3044 | buf = edi->buf->do_edsam; |
3045 | |
3046 | if (ed->bFirst) |
3047 | { |
3048 | /* initialize radacc radius for slope criterion */ |
3049 | buf->oldrad = calc_radius(&edi->vecs.radacc); |
3050 | } |
3051 | |
3052 | /* Copy the positions into buf->xc* arrays and after ED |
3053 | * feed back corrections to the official positions */ |
3054 | |
3055 | /* Broadcast the ED positions such that every node has all of them |
3056 | * Every node contributes its local positions xs and stores it in |
3057 | * the collective buf->xcoll array. Note that for edinr > 1 |
3058 | * xs could already have been modified by an earlier ED */ |
3059 | |
3060 | communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr)((cr)->nnodes > 1) ? buf->bUpdateShifts : TRUE1, xs, |
3061 | edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box); |
3062 | |
3063 | /* Only assembly reference positions if their indices differ from the average ones */ |
3064 | if (!edi->bRefEqAv) |
3065 | { |
3066 | communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr)((cr)->nnodes > 1) ? buf->bUpdateShifts : TRUE1, xs, |
3067 | edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box); |
3068 | } |
3069 | |
3070 | /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions. |
3071 | * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices |
3072 | * set bUpdateShifts=TRUE in the parallel case. */ |
3073 | buf->bUpdateShifts = FALSE0; |
3074 | |
3075 | /* Now all nodes have all of the ED positions in edi->sav->xcoll, |
3076 | * as well as the indices in edi->sav.anrs */ |
3077 | |
3078 | /* Fit the reference indices to the reference structure */ |
3079 | if (edi->bRefEqAv) |
3080 | { |
3081 | fit_to_reference(buf->xcoll, transvec, rotmat, edi); |
3082 | } |
3083 | else |
3084 | { |
3085 | fit_to_reference(buf->xc_ref, transvec, rotmat, edi); |
3086 | } |
3087 | |
3088 | /* Now apply the translation and rotation to the ED structure */ |
3089 | translate_and_rotate(buf->xcoll, edi->sav.nr, transvec, rotmat); |
3090 | |
3091 | /* Find out how well we fit to the reference (just for output steps) */ |
3092 | if (do_per_step(step, edi->outfrq) && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
3093 | { |
3094 | if (edi->bRefEqAv) |
3095 | { |
3096 | /* Indices of reference and average structures are identical, |
3097 | * thus we can calculate the rmsd to SREF using xcoll */ |
3098 | rmsdev = rmsd_from_structure(buf->xcoll, &edi->sref); |
3099 | } |
3100 | else |
3101 | { |
3102 | /* We have to translate & rotate the reference atoms first */ |
3103 | translate_and_rotate(buf->xc_ref, edi->sref.nr, transvec, rotmat); |
3104 | rmsdev = rmsd_from_structure(buf->xc_ref, &edi->sref); |
3105 | } |
3106 | } |
3107 | |
3108 | /* update radsam references, when required */ |
3109 | if (do_per_step(step, edi->maxedsteps) && step >= edi->presteps) |
3110 | { |
3111 | project(buf->xcoll, edi); |
3112 | rad_project(edi, buf->xcoll, &edi->vecs.radacc); |
3113 | rad_project(edi, buf->xcoll, &edi->vecs.radfix); |
3114 | buf->oldrad = -1.e5; |
3115 | } |
3116 | |
3117 | /* update radacc references, when required */ |
3118 | if (do_per_step(step, iupdate) && step >= edi->presteps) |
3119 | { |
3120 | edi->vecs.radacc.radius = calc_radius(&edi->vecs.radacc); |
3121 | if (edi->vecs.radacc.radius - buf->oldrad < edi->slope) |
3122 | { |
3123 | project(buf->xcoll, edi); |
3124 | rad_project(edi, buf->xcoll, &edi->vecs.radacc); |
3125 | buf->oldrad = 0.0; |
3126 | } |
3127 | else |
3128 | { |
3129 | buf->oldrad = edi->vecs.radacc.radius; |
3130 | } |
3131 | } |
3132 | |
3133 | /* apply the constraints */ |
3134 | if (step >= edi->presteps && ed_constraints(ed->eEDtype, edi)) |
3135 | { |
3136 | /* ED constraints should be applied already in the first MD step |
3137 | * (which is step 0), therefore we pass step+1 to the routine */ |
3138 | ed_apply_constraints(buf->xcoll, edi, step+1 - ir->init_step); |
3139 | } |
3140 | |
3141 | /* write to edo, when required */ |
3142 | if (do_per_step(step, edi->outfrq)) |
3143 | { |
3144 | project(buf->xcoll, edi); |
3145 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && !bSuppress) |
3146 | { |
3147 | write_edo(edi, ed->edo, rmsdev); |
3148 | } |
3149 | } |
3150 | |
3151 | /* Copy back the positions unless monitoring only */ |
3152 | if (ed_constraints(ed->eEDtype, edi)) |
3153 | { |
3154 | /* remove fitting */ |
3155 | rmfit(edi->sav.nr, buf->xcoll, transvec, rotmat); |
3156 | |
3157 | /* Copy the ED corrected positions into the coordinate array */ |
3158 | /* Each node copies its local part. In the serial case, nat_loc is the |
3159 | * total number of ED atoms */ |
3160 | for (i = 0; i < edi->sav.nr_loc; i++) |
3161 | { |
3162 | /* Unshift local ED coordinate and store in x_unsh */ |
3163 | ed_unshift_single_coord(box, buf->xcoll[edi->sav.c_ind[i]], |
3164 | buf->shifts_xcoll[edi->sav.c_ind[i]], x_unsh); |
3165 | |
3166 | /* dx is the ED correction to the positions: */ |
3167 | rvec_sub(x_unsh, xs[edi->sav.anrs_loc[i]], dx); |
3168 | |
3169 | if (v != NULL((void*)0)) |
3170 | { |
3171 | /* dv is the ED correction to the velocity: */ |
3172 | svmul(dt_1, dx, dv); |
3173 | /* apply the velocity correction: */ |
3174 | rvec_inc(v[edi->sav.anrs_loc[i]], dv); |
3175 | } |
3176 | /* Finally apply the position correction due to ED: */ |
3177 | copy_rvec(x_unsh, xs[edi->sav.anrs_loc[i]]); |
3178 | } |
3179 | } |
3180 | } /* END of if ( bNeedDoEdsam(edi) ) */ |
3181 | |
3182 | /* Prepare for the next ED group */ |
3183 | edi = edi->next_edi; |
3184 | |
3185 | } /* END of loop over ED groups */ |
3186 | |
3187 | ed->bFirst = FALSE0; |
3188 | } |