File: | gromacs/essentialdynamics/edsam.c |
Location: | line 2747, column 13 |
Description: | Dereference of null pointer |
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); | |||||
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 | } |