/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c

Bug Summary

File:	gromacs/essentialdynamics/edsam.c
Location:	line 2100, column 9
Description:	Value stored to 'rad' is never read

Annotated Source Code

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, 2DIM)(loc->omega) = save_calloc("loc->omega", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 581, (23), sizeof(*(loc->omega)));
582	snew(loc->om, 2DIM)(loc->om) = save_calloc("loc->om", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 582, (23), sizeof(*(loc->om)));
583	for (i = 0; i < 2*DIM3; i++)
584	{
585	snew(loc->omega[i], 2DIM)(loc->omega[i]) = save_calloc("loc->omega[i]", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 585, (23), sizeof(*(loc->omega[i])));
586	snew(loc->om[i], 2DIM)(loc->om[i]) = save_calloc("loc->om[i]", "/home/alexxy/Develop/gromacs/src/gromacs/essentialdynamics/edsam.c" , 586, (23), 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.5edi->flood.Eflsum; /* minus sign because Efl is negative, if restrain is on. */
847	}
848	else
849	{
850	Vfl = edi->flood.Efl != 0 ? edi->flood.Eflexp(-edi->flood.kT/2/edi->flood.Efl/edi->flood.alpha2sum) : 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.alpha2energyedi->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.tauedi->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]-txbox[XX0][XX0]-tybox[YY1][XX0]-tz*box[ZZ2][XX0];
1945	xu[YY1] = x[YY1]-tybox[YY1][YY1]-tzbox[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;
	Value stored to 'rad' is never read
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	}