?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Redefine the default boolean type to gmx_bool.
[gromacs.git] / src / kernel / repl_ex.c
bloba1cfc99367c4c8050a22a9c15afe56a0b2b2c5be
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <math.h>
41 #include "repl_ex.h"
42 #include "network.h"
43 #include "random.h"
44 #include "smalloc.h"
45 #include "physics.h"
46 #include "copyrite.h"
47 #include "macros.h"
48 #include "vec.h"
49 #include "names.h"
50 #include "mvdata.h"
51 #include "domdec.h"
52 #include "partdec.h"
53
54 typedef struct gmx_repl_ex {
55 int repl;
56 int nrepl;
57 real temp;
58 int type;
59 real *q;
60 gmx_bool bNPT;
61 real *pres;
62 int *ind;
63 int nst;
64 int seed;
65 int nattempt[2];
66 real *prob_sum;
67 int *nexchange;
68 } t_gmx_repl_ex;
69
70 enum { ereTEMP, ereLAMBDA, ereNR };
71 const char *erename[ereNR] = { "temperature", "lambda" };
72
73 static void repl_quantity(FILE *fplog,const gmx_multisim_t *ms,
74 struct gmx_repl_ex *re,int ere,real q)
75 {
76 real *qall;
77 gmx_bool bDiff;
78 int s;
79
80 snew(qall,ms->nsim);
81 qall[re->repl] = q;
82 gmx_sum_sim(ms->nsim,qall,ms);
83
84 bDiff = FALSE;
85 for(s=1; s<ms->nsim; s++)
86 if (qall[s] != qall[0])
87 bDiff = TRUE;
88
89 if (bDiff) {
90 if (re->type >= 0 && re->type < ereNR) {
91 gmx_fatal(FARGS,"For replica exchange both %s and %s differ",
92 erename[re->type],erename[ere]);
93 } else {
94 /* Set the replica exchange type and quantities */
95 re->type = ere;
96 snew(re->q,re->nrepl);
97 for(s=0; s<ms->nsim; s++)
98 re->q[s] = qall[s];
99 }
100 }
101
102 sfree(qall);
103 }
104
105 gmx_repl_ex_t init_replica_exchange(FILE *fplog,
106 const gmx_multisim_t *ms,
107 const t_state *state,
108 const t_inputrec *ir,
109 int nst,int init_seed)
110 {
111 real temp,pres;
112 int i,j,k;
113 struct gmx_repl_ex *re;
114
115 fprintf(fplog,"\nInitializing Replica Exchange\n");
116
117 if (ms == NULL || ms->nsim == 1)
118 gmx_fatal(FARGS,"Nothing to exchange with only one replica, maybe you forgot to set the -multi option of mdrun?");
119
120 snew(re,1);
121
122 re->repl = ms->sim;
123 re->nrepl = ms->nsim;
124
125 fprintf(fplog,"Repl There are %d replicas:\n",re->nrepl);
126
127 check_multi_int(fplog,ms,state->natoms,"the number of atoms");
128 check_multi_int(fplog,ms,ir->eI,"the integrator");
129 check_multi_int(fplog,ms,ir->init_step+ir->nsteps,"init_step+nsteps");
130 check_multi_int(fplog,ms,(ir->init_step+nst-1)/nst,
131 "first exchange step: init_step/-replex");
132 check_multi_int(fplog,ms,ir->etc,"the temperature coupling");
133 check_multi_int(fplog,ms,ir->opts.ngtc,
134 "the number of temperature coupling groups");
135 check_multi_int(fplog,ms,ir->epc,"the pressure coupling");
136 check_multi_int(fplog,ms,ir->efep,"free energy");
137
138 re->temp = ir->opts.ref_t[0];
139 for(i=1; (i<ir->opts.ngtc); i++) {
140 if (ir->opts.ref_t[i] != re->temp) {
141 fprintf(fplog,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
142 fprintf(stderr,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
143 }
144 }
145
146 re->type = -1;
147 for(i=0; i<ereNR; i++)
148 {
149 switch (i)
150 {
151 case ereTEMP:
152 repl_quantity(fplog,ms,re,i,re->temp);
153 break;
154 case ereLAMBDA:
155 if (ir->efep != efepNO)
156 {
157 repl_quantity(fplog,ms,re,i,ir->init_lambda);
158 }
159 break;
160 default:
161 gmx_incons("Unknown replica exchange quantity");
162 }
163 }
164 if (re->type == -1)
165 {
166 gmx_fatal(FARGS,"The properties of the %d systems are all the same, there is nothing to exchange",re->nrepl);
167 }
168
169 switch (re->type)
170 {
171 case ereTEMP:
172 please_cite(fplog,"Hukushima96a");
173 if (ir->epc != epcNO)
174 {
175 re->bNPT = TRUE;
176 fprintf(fplog,"Repl Using Constant Pressure REMD.\n");
177 please_cite(fplog,"Okabe2001a");
178 }
179 if (ir->etc == etcBERENDSEN)
180 {
181 gmx_fatal(FARGS,"REMD with the %s thermostat does not produce correct potential energy distributions, consider using the %s thermostat instead",
182 ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
183 }
184 break;
185 case ereLAMBDA:
186 if (ir->delta_lambda != 0)
187 {
188 gmx_fatal(FARGS,"delta_lambda is not zero");
189 }
190 break;
191 }
192
193 if (re->bNPT) {
194 snew(re->pres,re->nrepl);
195 if (ir->epct == epctSURFACETENSION) {
196 pres = ir->ref_p[ZZ][ZZ];
197 } else {
198 pres = 0;
199 j = 0;
200 for(i=0; i<DIM; i++)
201 if (ir->compress[i][i] != 0) {
202 pres += ir->ref_p[i][i];
203 j++;
204 }
205 pres /= j;
206 }
207 re->pres[re->repl] = pres;
208 gmx_sum_sim(re->nrepl,re->pres,ms);
209 }
210
211 snew(re->ind,re->nrepl);
212 /* Make an index for increasing temperature order */
213 for(i=0; i<re->nrepl; i++)
214 re->ind[i] = i;
215 for(i=0; i<re->nrepl; i++) {
216 for(j=i+1; j<re->nrepl; j++) {
217 if (re->q[re->ind[j]] < re->q[re->ind[i]]) {
218 k = re->ind[i];
219 re->ind[i] = re->ind[j];
220 re->ind[j] = k;
221 } else if (re->q[re->ind[j]] == re->q[re->ind[i]]) {
222 gmx_fatal(FARGS,"Two replicas have identical %ss",erename[re->type]);
223 }
224 }
225 }
226 fprintf(fplog,"Repl ");
227 for(i=0; i<re->nrepl; i++)
228 fprintf(fplog," %3d ",re->ind[i]);
229 switch (re->type) {
230 case ereTEMP:
231 fprintf(fplog,"\nRepl T");
232 for(i=0; i<re->nrepl; i++)
233 fprintf(fplog," %5.1f",re->q[re->ind[i]]);
234 break;
235 case ereLAMBDA:
236 fprintf(fplog,"\nRepl l");
237 for(i=0; i<re->nrepl; i++)
238 fprintf(fplog," %5.3f",re->q[re->ind[i]]);
239 break;
240 default:
241 gmx_incons("Unknown replica exchange quantity");
242 }
243 if (re->bNPT) {
244 fprintf(fplog,"\nRepl p");
245 for(i=0; i<re->nrepl; i++)
246 {
247 fprintf(fplog," %5.2f",re->pres[re->ind[i]]);
248 }
249
250 for(i=0; i<re->nrepl; i++)
251 {
252 if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i-1]]))
253 {
254 gmx_fatal(FARGS,"The reference pressure decreases with increasing temperature");
255 }
256 }
257 }
258 fprintf(fplog,"\nRepl ");
259
260 re->nst = nst;
261 if (init_seed == -1) {
262 if (MASTERSIM(ms))
263 re->seed = make_seed();
264 else
265 re->seed = 0;
266 gmx_sumi_sim(1,&(re->seed),ms);
267 } else {
268 re->seed = init_seed;
269 }
270 fprintf(fplog,"\nRepl exchange interval: %d\n",re->nst);
271 fprintf(fplog,"\nRepl random seed: %d\n",re->seed);
272
273 re->nattempt[0] = 0;
274 re->nattempt[1] = 0;
275 snew(re->prob_sum,re->nrepl);
276 snew(re->nexchange,re->nrepl);
277
278 fprintf(fplog,
279 "Repl below: x=exchange, pr=probability\n");
280
281 return re;
282 }
283
284 static void exchange_reals(const gmx_multisim_t *ms,int b,real *v,int n)
285 {
286 real *buf;
287 int i;
288
289 if (v) {
290 snew(buf,n);
291 #ifdef GMX_MPI
292 /*
293 MPI_Sendrecv(v, n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
294 buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
295 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
296 */
297 {
298 MPI_Request mpi_req;
299
300 MPI_Isend(v,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
301 ms->mpi_comm_masters,&mpi_req);
302 MPI_Recv(buf,n*sizeof(real),MPI_BYTE,MSRANK(ms,b),0,
303 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
304 MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
305 }
306 #endif
307 for(i=0; i<n; i++)
308 v[i] = buf[i];
309 sfree(buf);
310 }
311 }
312
313 static void exchange_doubles(const gmx_multisim_t *ms,int b,double *v,int n)
314 {
315 double *buf;
316 int i;
317
318 if (v) {
319 snew(buf,n);
320 #ifdef GMX_MPI
321 /*
322 MPI_Sendrecv(v, n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
323 buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
324 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
325 */
326 {
327 MPI_Request mpi_req;
328
329 MPI_Isend(v,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
330 ms->mpi_comm_masters,&mpi_req);
331 MPI_Recv(buf,n*sizeof(double),MPI_BYTE,MSRANK(ms,b),0,
332 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
333 MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
334 }
335 #endif
336 for(i=0; i<n; i++)
337 v[i] = buf[i];
338 sfree(buf);
339 }
340 }
341
342 static void exchange_rvecs(const gmx_multisim_t *ms,int b,rvec *v,int n)
343 {
344 rvec *buf;
345 int i;
346
347 if (v) {
348 snew(buf,n);
349 #ifdef GMX_MPI
350 /*
351 MPI_Sendrecv(v[0], n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
352 buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
353 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
354 */
355 {
356 MPI_Request mpi_req;
357
358 MPI_Isend(v[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
359 ms->mpi_comm_masters,&mpi_req);
360 MPI_Recv(buf[0],n*sizeof(rvec),MPI_BYTE,MSRANK(ms,b),0,
361 ms->mpi_comm_masters,MPI_STATUS_IGNORE);
362 MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
363 }
364 #endif
365 for(i=0; i<n; i++)
366 copy_rvec(buf[i],v[i]);
367 sfree(buf);
368 }
369 }
370
371 static void exchange_state(const gmx_multisim_t *ms,int b,t_state *state)
372 {
373 /* When t_state changes, this code should be updated. */
374 int ngtc,nnhpres;
375 ngtc = state->ngtc * state->nhchainlength;
376 nnhpres = state->nnhpres* state->nhchainlength;
377 exchange_rvecs(ms,b,state->box,DIM);
378 exchange_rvecs(ms,b,state->box_rel,DIM);
379 exchange_rvecs(ms,b,state->boxv,DIM);
380 exchange_reals(ms,b,&(state->veta),1);
381 exchange_reals(ms,b,&(state->vol0),1);
382 exchange_rvecs(ms,b,state->svir_prev,DIM);
383 exchange_rvecs(ms,b,state->fvir_prev,DIM);
384 exchange_rvecs(ms,b,state->pres_prev,DIM);
385 exchange_doubles(ms,b,state->nosehoover_xi,ngtc);
386 exchange_doubles(ms,b,state->nosehoover_vxi,ngtc);
387 exchange_doubles(ms,b,state->nhpres_xi,nnhpres);
388 exchange_doubles(ms,b,state->nhpres_vxi,nnhpres);
389 exchange_doubles(ms,b,state->therm_integral,state->ngtc);
390 exchange_rvecs(ms,b,state->x,state->natoms);
391 exchange_rvecs(ms,b,state->v,state->natoms);
392 exchange_rvecs(ms,b,state->sd_X,state->natoms);
393 }
394
395 static void copy_rvecs(rvec *s,rvec *d,int n)
396 {
397 int i;
398
399 if (d != NULL)
400 {
401 for(i=0; i<n; i++)
402 {
403 copy_rvec(s[i],d[i]);
404 }
405 }
406 }
407
408 static void copy_doubles(const double *s,double *d,int n)
409 {
410 int i;
411
412 if (d != NULL)
413 {
414 for(i=0; i<n; i++)
415 {
416 d[i] = s[i];
417 }
418 }
419 }
420
421 #define scopy_rvecs(v,n) copy_rvecs(state->v,state_local->v,n);
422 #define scopy_doubles(v,n) copy_doubles(state->v,state_local->v,n);
423
424 static void copy_state_nonatomdata(t_state *state,t_state *state_local)
425 {
426 /* When t_state changes, this code should be updated. */
427 int ngtc,nnhpres;
428 ngtc = state->ngtc * state->nhchainlength;
429 nnhpres = state->nnhpres* state->nhchainlength;
430 scopy_rvecs(box,DIM);
431 scopy_rvecs(box_rel,DIM);
432 scopy_rvecs(boxv,DIM);
433 state_local->veta = state->veta;
434 state_local->vol0 = state->vol0;
435 scopy_rvecs(svir_prev,DIM);
436 scopy_rvecs(fvir_prev,DIM);
437 scopy_rvecs(pres_prev,DIM);
438 scopy_doubles(nosehoover_xi,ngtc);
439 scopy_doubles(nosehoover_vxi,ngtc);
440 scopy_doubles(nhpres_xi,nnhpres);
441 scopy_doubles(nhpres_vxi,nnhpres);
442 scopy_doubles(therm_integral,state->ngtc);
443 scopy_rvecs(x,state->natoms);
444 scopy_rvecs(v,state->natoms);
445 scopy_rvecs(sd_X,state->natoms);
446 }
447
448 static void scale_velocities(t_state *state,real fac)
449 {
450 int i;
451
452 if (state->v)
453 for(i=0; i<state->natoms; i++)
454 svmul(fac,state->v[i],state->v[i]);
455 }
456
457 static void pd_collect_state(const t_commrec *cr,t_state *state)
458 {
459 int shift;
460
461 if (debug)
462 fprintf(debug,"Collecting state before exchange\n");
463 shift = cr->nnodes - cr->npmenodes - 1;
464 move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->x,NULL,shift,NULL);
465 if (state->v)
466 move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->v,NULL,shift,NULL);
467 if (state->sd_X)
468 move_rvecs(cr,FALSE,FALSE,GMX_LEFT,GMX_RIGHT,state->sd_X,NULL,shift,NULL);
469 }
470
471 static void print_ind(FILE *fplog,const char *leg,int n,int *ind,gmx_bool *bEx)
472 {
473 int i;
474
475 fprintf(fplog,"Repl %2s %2d",leg,ind[0]);
476 for(i=1; i<n; i++) {
477 fprintf(fplog," %c %2d",(bEx!=0 && bEx[i]) ? 'x' : ' ',ind[i]);
478 }
479 fprintf(fplog,"\n");
480 }
481
482 static void print_prob(FILE *fplog,const char *leg,int n,real *prob)
483 {
484 int i;
485 char buf[8];
486
487 fprintf(fplog,"Repl %2s ",leg);
488 for(i=1; i<n; i++) {
489 if (prob[i] >= 0) {
490 sprintf(buf,"%4.2f",prob[i]);
491 fprintf(fplog," %3s",buf[0]=='1' ? "1.0" : buf+1);
492 } else {
493 fprintf(fplog," ");
494 }
495 }
496 fprintf(fplog,"\n");
497 }
498
499 static void print_count(FILE *fplog,const char *leg,int n,int *count)
500 {
501 int i;
502
503 fprintf(fplog,"Repl %2s ",leg);
504 for(i=1; i<n; i++) {
505 fprintf(fplog," %4d",count[i]);
506 }
507 fprintf(fplog,"\n");
508 }
509
510 static int get_replica_exchange(FILE *fplog,const gmx_multisim_t *ms,
511 struct gmx_repl_ex *re,real *ener,real vol,
512 int step,real time)
513 {
514 int m,i,a,b;
515 real *Epot=NULL,*Vol=NULL,*dvdl=NULL,*prob;
516 real ediff=0,delta=0,dpV=0,betaA=0,betaB=0;
517 gmx_bool *bEx,bPrint;
518 int exchange;
519
520 fprintf(fplog,"Replica exchange at step %d time %g\n",step,time);
521
522 switch (re->type) {
523 case ereTEMP:
524 snew(Epot,re->nrepl);
525 snew(Vol,re->nrepl);
526 Epot[re->repl] = ener[F_EPOT];
527 Vol[re->repl] = vol;
528 gmx_sum_sim(re->nrepl,Epot,ms);
529 gmx_sum_sim(re->nrepl,Vol,ms);
530 break;
531 case ereLAMBDA:
532 snew(dvdl,re->nrepl);
533 dvdl[re->repl] = ener[F_DVDL];
534 gmx_sum_sim(re->nrepl,dvdl,ms);
535 break;
536 }
537
538 snew(bEx,re->nrepl);
539 snew(prob,re->nrepl);
540
541 exchange = -1;
542 m = (step / re->nst) % 2;
543 for(i=1; i<re->nrepl; i++) {
544 a = re->ind[i-1];
545 b = re->ind[i];
546 bPrint = (re->repl==a || re->repl==b);
547 if (i % 2 == m) {
548 switch (re->type) {
549 case ereTEMP:
550 /* Use equations from:
551 * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
552 */
553 ediff = Epot[b] - Epot[a];
554 betaA = 1.0/(re->q[a]*BOLTZ);
555 betaB = 1.0/(re->q[b]*BOLTZ);
556 delta = (betaA - betaB)*ediff;
557 break;
558 case ereLAMBDA:
559 /* Here we exchange based on a linear extrapolation of dV/dlambda.
560 * We would like to have the real energies
561 * from foreign lambda calculations.
562 */
563 ediff = (dvdl[a] - dvdl[b])*(re->q[b] - re->q[a]);
564 delta = ediff/(BOLTZ*re->temp);
565 break;
566 default:
567 gmx_incons("Unknown replica exchange quantity");
568 }
569 if (bPrint)
570 fprintf(fplog,"Repl %d <-> %d dE = %10.3e",a,b,delta);
571 if (re->bNPT) {
572 dpV = (betaA*re->pres[a]-betaB*re->pres[b])*(Vol[b]-Vol[a])/PRESFAC;
573 if (bPrint)
574 fprintf(fplog," dpV = %10.3e d = %10.3e",dpV,delta + dpV);
575 delta += dpV;
576 }
577 if (bPrint)
578 fprintf(fplog,"\n");
579 if (delta <= 0) {
580 prob[i] = 1;
581 bEx[i] = TRUE;
582 } else {
583 if (delta > 100)
584 prob[i] = 0;
585 else
586 prob[i] = exp(-delta);
587 bEx[i] = (rando(&(re->seed)) < prob[i]);
588 }
589 re->prob_sum[i] += prob[i];
590 if (bEx[i]) {
591 if (a == re->repl) {
592 exchange = b;
593 } else if (b == re->repl) {
594 exchange = a;
595 }
596 re->nexchange[i]++;
597 }
598 } else {
599 prob[i] = -1;
600 bEx[i] = FALSE;
601 }
602 }
603 print_ind(fplog,"ex",re->nrepl,re->ind,bEx);
604 print_prob(fplog,"pr",re->nrepl,prob);
605 fprintf(fplog,"\n");
606
607 sfree(bEx);
608 sfree(prob);
609 sfree(Epot);
610 sfree(Vol);
611 sfree(dvdl);
612
613 re->nattempt[m]++;
614
615 return exchange;
616 }
617
618 static void write_debug_x(t_state *state)
619 {
620 int i;
621
622 if (debug) {
623 for(i=0; i<state->natoms; i+=10)
624 fprintf(debug,"dx %5d %10.5f %10.5f %10.5f\n",i,state->x[i][XX],state->x[i][YY],state->x[i][ZZ]);
625 }
626 }
627
628 gmx_bool replica_exchange(FILE *fplog,const t_commrec *cr,struct gmx_repl_ex *re,
629 t_state *state,real *ener,
630 t_state *state_local,
631 int step,real time)
632 {
633 gmx_multisim_t *ms;
634 int exchange=-1,shift;
635 gmx_bool bExchanged=FALSE;
636
637 ms = cr->ms;
638
639 if (MASTER(cr))
640 {
641 exchange = get_replica_exchange(fplog,ms,re,ener,det(state->box),
642 step,time);
643 bExchanged = (exchange >= 0);
644 }
645
646 if (PAR(cr))
647 {
648 #ifdef GMX_MPI
649 MPI_Bcast(&bExchanged,sizeof(gmx_bool),MPI_BYTE,MASTERRANK(cr),
650 cr->mpi_comm_mygroup);
651 #endif
652 }
653
654 if (bExchanged)
655 {
656 /* Exchange the states */
657
658 if (PAR(cr))
659 {
660 /* Collect the global state on the master node */
661 if (DOMAINDECOMP(cr))
662 {
663 dd_collect_state(cr->dd,state_local,state);
664 }
665 else
666 {
667 pd_collect_state(cr,state);
668 }
669 }
670
671 if (MASTER(cr))
672 {
673 /* Exchange the global states between the master nodes */
674 if (debug)
675 {
676 fprintf(debug,"Exchanging %d with %d\n",ms->sim,exchange);
677 }
678 exchange_state(ms,exchange,state);
679
680 if (re->type == ereTEMP)
681 {
682 scale_velocities(state,sqrt(re->q[ms->sim]/re->q[exchange]));
683 }
684 }
685
686 /* With domain decomposition the global state is distributed later */
687 if (!DOMAINDECOMP(cr))
688 {
689 /* Copy the global state to the local state data structure */
690 copy_state_nonatomdata(state,state_local);
691
692 if (PAR(cr))
693 {
694 bcast_state(cr,state,FALSE);
695 }
696 }
697 }
698
699 return bExchanged;
700 }
701
702 void print_replica_exchange_statistics(FILE *fplog,struct gmx_repl_ex *re)
703 {
704 real *prob;
705 int i;
706
707 fprintf(fplog,"\nReplica exchange statistics\n");
708 fprintf(fplog,"Repl %d attempts, %d odd, %d even\n",
709 re->nattempt[0]+re->nattempt[1],re->nattempt[1],re->nattempt[0]);
710
711 snew(prob,re->nrepl);
712
713 fprintf(fplog,"Repl average probabilities:\n");
714 for(i=1; i<re->nrepl; i++) {
715 if (re->nattempt[i%2] == 0)
716 prob[i] = 0;
717 else
718 prob[i] = re->prob_sum[i]/re->nattempt[i%2];
719 }
720 print_ind(fplog,"",re->nrepl,re->ind,NULL);
721 print_prob(fplog,"",re->nrepl,prob);
722
723 fprintf(fplog,"Repl number of exchanges:\n");
724 print_ind(fplog,"",re->nrepl,re->ind,NULL);
725 print_count(fplog,"",re->nrepl,re->nexchange);
726
727 fprintf(fplog,"Repl average number of exchanges:\n");
728 for(i=1; i<re->nrepl; i++) {
729 if (re->nattempt[i%2] == 0)
730 prob[i] = 0;
731 else
732 prob[i] = ((real)re->nexchange[i])/re->nattempt[i%2];
733 }
734 print_ind(fplog,"",re->nrepl,re->ind,NULL);
735 print_prob(fplog,"",re->nrepl,prob);
736
737 sfree(prob);
738
739 fprintf(fplog,"\n");
740 }
The ultimate molecular dynamics simulation package