File: | gromacs/essentialdynamics/edsam.c |
Location: | line 2732, column 33 |
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 | } |