1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
51 #include "chargegroup.h"
56 #include "gmx_fatal.h"
71 #include "gmx_random.h"
76 #include "gmx_wallcycle.h"
77 #include "mtop_util.h"
83 #include "gmx_x86_sse2.h"
87 static void global_max(t_commrec
*cr
,int *n
)
93 gmx_sumi(cr
->nnodes
,sum
,cr
);
94 for(i
=0; i
<cr
->nnodes
; i
++)
100 static void realloc_bins(double **bin
,int *nbin
,int nbin_new
)
104 if (nbin_new
!= *nbin
) {
105 srenew(*bin
,nbin_new
);
106 for(i
=*nbin
; i
<nbin_new
; i
++)
112 double do_tpi(FILE *fplog
,t_commrec
*cr
,
113 int nfile
, const t_filenm fnm
[],
114 const output_env_t oenv
, gmx_bool bVerbose
,gmx_bool bCompact
,
116 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
118 t_inputrec
*inputrec
,
119 gmx_mtop_t
*top_global
,t_fcdata
*fcd
,
122 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
125 int repl_ex_nst
, int repl_ex_nex
, int repl_ex_seed
,
127 real cpt_period
,real max_hours
,
128 const char *deviceOptions
,
130 gmx_runtime_t
*runtime
)
132 const char *TPI
="Test Particle Insertion";
134 gmx_groups_t
*groups
;
135 gmx_enerdata_t
*enerd
;
137 real lambda
,t
,temp
,beta
,drmax
,epot
;
138 double embU
,sum_embU
,*sum_UgembU
,V
,V_all
,VembU_all
;
141 gmx_bool bDispCorr
,bCharge
,bRFExcl
,bNotLastFrame
,bStateChanged
,bNS
,bOurStep
;
142 tensor force_vir
,shake_vir
,vir
,pres
;
143 int cg_tp
,a_tp0
,a_tp1
,ngid
,gid_tp
,nener
,e
;
145 rvec mu_tot
,x_init
,dx
,x_tp
;
146 int nnodes
,frame
,nsteps
,step
;
150 char *ptr
,*dump_pdb
,**leg
,str
[STRLEN
],str2
[STRLEN
];
151 double dbl
,dump_ener
;
154 real
*mass_cavity
=NULL
,mass_tot
;
156 double invbinw
,*bin
,refvolshift
,logV
,bUlogV
;
157 real dvdl
,prescorr
,enercorr
,dvdlcorr
;
158 gmx_bool bEnergyOutOfBounds
;
159 const char *tpid_leg
[2]={"direct","reweighted"};
161 /* Since there is no upper limit to the insertion energies,
162 * we need to set an upper limit for the distribution output.
164 real bU_bin_limit
= 50;
165 real bU_logV_bin_limit
= bU_bin_limit
+ 10;
169 top
= gmx_mtop_generate_local_top(top_global
,inputrec
);
171 groups
= &top_global
->groups
;
173 bCavity
= (inputrec
->eI
== eiTPIC
);
175 ptr
= getenv("GMX_TPIC_MASSES");
179 /* Read (multiple) masses from env var GMX_TPIC_MASSES,
180 * The center of mass of the last atoms is then used for TPIC.
183 while (sscanf(ptr
,"%lf%n",&dbl
,&i
) > 0) {
184 srenew(mass_cavity
,nat_cavity
+1);
185 mass_cavity
[nat_cavity
] = dbl
;
186 fprintf(fplog
,"mass[%d] = %f\n",
187 nat_cavity
+1,mass_cavity
[nat_cavity
]);
192 gmx_fatal(FARGS
,"Found %d masses in GMX_TPIC_MASSES",nat_cavity
);
197 init_em(fplog,TPI,inputrec,&lambda,nrnb,mu_tot,
198 state->box,fr,mdatoms,top,cr,nfile,fnm,NULL,NULL);*/
199 /* We never need full pbc for TPI */
201 /* Determine the temperature for the Boltzmann weighting */
202 temp
= inputrec
->opts
.ref_t
[0];
204 for(i
=1; (i
<inputrec
->opts
.ngtc
); i
++) {
205 if (inputrec
->opts
.ref_t
[i
] != temp
) {
206 fprintf(fplog
,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
207 fprintf(stderr
,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
211 "\n The temperature for test particle insertion is %.3f K\n\n",
214 beta
= 1.0/(BOLTZ
*temp
);
216 /* Number of insertions per frame */
217 nsteps
= inputrec
->nsteps
;
219 /* Use the same neighborlist with more insertions points
220 * in a sphere of radius drmax around the initial point
222 /* This should be a proper mdp parameter */
223 drmax
= inputrec
->rtpi
;
225 /* An environment variable can be set to dump all configurations
226 * to pdb with an insertion energy <= this value.
228 dump_pdb
= getenv("GMX_TPI_DUMP");
231 sscanf(dump_pdb
,"%lf",&dump_ener
);
233 atoms2md(top_global
,inputrec
,0,NULL
,0,top_global
->natoms
,mdatoms
);
234 update_mdatoms(mdatoms
,inputrec
->fepvals
->init_lambda
);
237 init_enerdata(groups
->grps
[egcENER
].nr
,inputrec
->fepvals
->n_lambda
,enerd
);
238 snew(f
,top_global
->natoms
);
240 /* Print to log file */
241 runtime_start(runtime
);
242 print_date_and_time(fplog
,cr
->nodeid
,
243 "Started Test Particle Insertion",runtime
);
244 wallcycle_start(wcycle
,ewcRUN
);
246 /* The last charge group is the group to be inserted */
247 cg_tp
= top
->cgs
.nr
- 1;
248 a_tp0
= top
->cgs
.index
[cg_tp
];
249 a_tp1
= top
->cgs
.index
[cg_tp
+1];
251 fprintf(debug
,"TPI cg %d, atoms %d-%d\n",cg_tp
,a_tp0
,a_tp1
);
252 if (a_tp1
- a_tp0
> 1 &&
253 (inputrec
->rlist
< inputrec
->rcoulomb
||
254 inputrec
->rlist
< inputrec
->rvdw
))
255 gmx_fatal(FARGS
,"Can not do TPI for multi-atom molecule with a twin-range cut-off");
256 snew(x_mol
,a_tp1
-a_tp0
);
258 bDispCorr
= (inputrec
->eDispCorr
!= edispcNO
);
260 for(i
=a_tp0
; i
<a_tp1
; i
++) {
261 /* Copy the coordinates of the molecule to be insterted */
262 copy_rvec(state
->x
[i
],x_mol
[i
-a_tp0
]);
263 /* Check if we need to print electrostatic energies */
264 bCharge
|= (mdatoms
->chargeA
[i
] != 0 ||
265 (mdatoms
->chargeB
&& mdatoms
->chargeB
[i
] != 0));
267 bRFExcl
= (bCharge
&& EEL_RF(fr
->eeltype
) && fr
->eeltype
!=eelRF_NEC
);
269 calc_cgcm(fplog
,cg_tp
,cg_tp
+1,&(top
->cgs
),state
->x
,fr
->cg_cm
);
271 if (norm(fr
->cg_cm
[cg_tp
]) > 0.5*inputrec
->rlist
&& fplog
) {
272 fprintf(fplog
, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
273 fprintf(stderr
,"WARNING: Your TPI molecule is not centered at 0,0,0\n");
276 /* Center the molecule to be inserted at zero */
277 for(i
=0; i
<a_tp1
-a_tp0
; i
++)
278 rvec_dec(x_mol
[i
],fr
->cg_cm
[cg_tp
]);
282 fprintf(fplog
,"\nWill insert %d atoms %s partial charges\n",
283 a_tp1
-a_tp0
,bCharge
? "with" : "without");
285 fprintf(fplog
,"\nWill insert %d times in each frame of %s\n",
286 nsteps
,opt2fn("-rerun",nfile
,fnm
));
291 if (inputrec
->nstlist
> 1)
293 if (drmax
==0 && a_tp1
-a_tp0
==1)
295 gmx_fatal(FARGS
,"Re-using the neighborlist %d times for insertions of a single atom in a sphere of radius %f does not make sense",inputrec
->nstlist
,drmax
);
299 fprintf(fplog
,"Will use the same neighborlist for %d insertions in a sphere of radius %f\n",inputrec
->nstlist
,drmax
);
307 fprintf(fplog
,"Will insert randomly in a sphere of radius %f around the center of the cavity\n",drmax
);
311 ngid
= groups
->grps
[egcENER
].nr
;
312 gid_tp
= GET_CGINFO_GID(fr
->cginfo
[cg_tp
]);
320 if (EEL_FULL(fr
->eeltype
))
323 snew(sum_UgembU
,nener
);
325 /* Initialize random generator */
326 tpi_rand
= gmx_rng_init(inputrec
->ld_seed
);
329 fp_tpi
= xvgropen(opt2fn("-tpi",nfile
,fnm
),
330 "TPI energies","Time (ps)",
331 "(kJ mol\\S-1\\N) / (nm\\S3\\N)",oenv
);
332 xvgr_subtitle(fp_tpi
,"f. are averages over one frame",oenv
);
335 sprintf(str
,"-kT log(<Ve\\S-\\betaU\\N>/<V>)");
336 leg
[e
++] = strdup(str
);
337 sprintf(str
,"f. -kT log<e\\S-\\betaU\\N>");
338 leg
[e
++] = strdup(str
);
339 sprintf(str
,"f. <e\\S-\\betaU\\N>");
340 leg
[e
++] = strdup(str
);
342 leg
[e
++] = strdup(str
);
343 sprintf(str
,"f. <Ue\\S-\\betaU\\N>");
344 leg
[e
++] = strdup(str
);
345 for(i
=0; i
<ngid
; i
++) {
346 sprintf(str
,"f. <U\\sVdW %s\\Ne\\S-\\betaU\\N>",
347 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
348 leg
[e
++] = strdup(str
);
351 sprintf(str
,"f. <U\\sdisp c\\Ne\\S-\\betaU\\N>");
352 leg
[e
++] = strdup(str
);
355 for(i
=0; i
<ngid
; i
++) {
356 sprintf(str
,"f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
357 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
358 leg
[e
++] = strdup(str
);
361 sprintf(str
,"f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
362 leg
[e
++] = strdup(str
);
364 if (EEL_FULL(fr
->eeltype
)) {
365 sprintf(str
,"f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
366 leg
[e
++] = strdup(str
);
369 xvgr_legend(fp_tpi
,4+nener
,(const char**)leg
,oenv
);
370 for(i
=0; i
<4+nener
; i
++)
382 bNotLastFrame
= read_first_frame(oenv
,&status
,opt2fn("-rerun",nfile
,fnm
),
383 &rerun_fr
,TRX_NEED_X
);
386 if (rerun_fr
.natoms
- (bCavity
? nat_cavity
: 0) !=
387 mdatoms
->nr
- (a_tp1
- a_tp0
))
388 gmx_fatal(FARGS
,"Number of atoms in trajectory (%d)%s "
389 "is not equal the number in the run input file (%d) "
390 "minus the number of atoms to insert (%d)\n",
391 rerun_fr
.natoms
,bCavity
? " minus one" : "",
392 mdatoms
->nr
,a_tp1
-a_tp0
);
394 refvolshift
= log(det(rerun_fr
.box
));
397 /* Make sure we don't detect SSE overflow generated before this point */
398 gmx_mm_check_and_reset_overflow();
401 while (bNotLastFrame
)
403 lambda
= rerun_fr
.lambda
;
407 for(e
=0; e
<nener
; e
++)
412 /* Copy the coordinates from the input trajectory */
413 for(i
=0; i
<rerun_fr
.natoms
; i
++)
415 copy_rvec(rerun_fr
.x
[i
],state
->x
[i
]);
417 copy_mat(rerun_fr
.box
,state
->box
);
422 bStateChanged
= TRUE
;
424 for(step
=0; step
<nsteps
; step
++)
426 /* In parallel all nodes generate all random configurations.
427 * In that way the result is identical to a single cpu tpi run.
431 /* Random insertion in the whole volume */
432 bNS
= (step
% inputrec
->nstlist
== 0);
435 /* Generate a random position in the box */
436 x_init
[XX
] = gmx_rng_uniform_real(tpi_rand
)*state
->box
[XX
][XX
];
437 x_init
[YY
] = gmx_rng_uniform_real(tpi_rand
)*state
->box
[YY
][YY
];
438 x_init
[ZZ
] = gmx_rng_uniform_real(tpi_rand
)*state
->box
[ZZ
][ZZ
];
440 if (inputrec
->nstlist
== 1)
442 copy_rvec(x_init
,x_tp
);
446 /* Generate coordinates within |dx|=drmax of x_init */
449 dx
[XX
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
450 dx
[YY
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
451 dx
[ZZ
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
453 while (norm2(dx
) > drmax
*drmax
);
454 rvec_add(x_init
,dx
,x_tp
);
459 /* Random insertion around a cavity location
460 * given by the last coordinate of the trajectory.
466 /* Copy the location of the cavity */
467 copy_rvec(rerun_fr
.x
[rerun_fr
.natoms
-1],x_init
);
471 /* Determine the center of mass of the last molecule */
474 for(i
=0; i
<nat_cavity
; i
++)
479 mass_cavity
[i
]*rerun_fr
.x
[rerun_fr
.natoms
-nat_cavity
+i
][d
];
481 mass_tot
+= mass_cavity
[i
];
485 x_init
[d
] /= mass_tot
;
489 /* Generate coordinates within |dx|=drmax of x_init */
492 dx
[XX
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
493 dx
[YY
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
494 dx
[ZZ
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
496 while (norm2(dx
) > drmax
*drmax
);
497 rvec_add(x_init
,dx
,x_tp
);
500 if (a_tp1
- a_tp0
== 1)
502 /* Insert a single atom, just copy the insertion location */
503 copy_rvec(x_tp
,state
->x
[a_tp0
]);
507 /* Copy the coordinates from the top file */
508 for(i
=a_tp0
; i
<a_tp1
; i
++)
510 copy_rvec(x_mol
[i
-a_tp0
],state
->x
[i
]);
512 /* Rotate the molecule randomly */
513 rotate_conf(a_tp1
-a_tp0
,state
->x
+a_tp0
,NULL
,
514 2*M_PI
*gmx_rng_uniform_real(tpi_rand
),
515 2*M_PI
*gmx_rng_uniform_real(tpi_rand
),
516 2*M_PI
*gmx_rng_uniform_real(tpi_rand
));
517 /* Shift to the insertion location */
518 for(i
=a_tp0
; i
<a_tp1
; i
++)
520 rvec_inc(state
->x
[i
],x_tp
);
524 /* Check if this insertion belongs to this node */
528 switch (inputrec
->eI
)
531 bOurStep
= ((step
/ inputrec
->nstlist
) % nnodes
== cr
->nodeid
);
534 bOurStep
= (step
% nnodes
== cr
->nodeid
);
537 gmx_fatal(FARGS
,"Unknown integrator %s",ei_names
[inputrec
->eI
]);
542 /* Clear some matrix variables */
543 clear_mat(force_vir
);
544 clear_mat(shake_vir
);
548 /* Set the charge group center of mass of the test particle */
549 copy_rvec(x_init
,fr
->cg_cm
[top
->cgs
.nr
-1]);
551 /* Calc energy (no forces) on new positions.
552 * Since we only need the intermolecular energy
553 * and the RF exclusion terms of the inserted molecule occur
554 * within a single charge group we can pass NULL for the graph.
555 * This also avoids shifts that would move charge groups
558 * Some checks above ensure than we can not have
559 * twin-range interactions together with nstlist > 1,
560 * therefore we do not need to remember the LR energies.
562 /* Make do_force do a single node force calculation */
564 do_force(fplog
,cr
,inputrec
,
565 step
,nrnb
,wcycle
,top
,top_global
,&top_global
->groups
,
566 state
->box
,state
->x
,&state
->hist
,
567 f
,force_vir
,mdatoms
,enerd
,fcd
,
569 NULL
,fr
,NULL
,mu_tot
,t
,NULL
,NULL
,FALSE
,
570 GMX_FORCE_NONBONDED
| GMX_FORCE_ENERGY
|
571 (bNS
? GMX_FORCE_DYNAMICBOX
| GMX_FORCE_NS
| GMX_FORCE_DOLR
: 0) |
572 (bStateChanged
? GMX_FORCE_STATECHANGED
: 0));
574 bStateChanged
= FALSE
;
577 /* Calculate long range corrections to pressure and energy */
578 calc_dispcorr(fplog
,inputrec
,fr
,step
,top_global
->natoms
,state
->box
,
579 lambda
,pres
,vir
,&prescorr
,&enercorr
,&dvdlcorr
);
580 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
581 enerd
->term
[F_DISPCORR
] = enercorr
;
582 enerd
->term
[F_EPOT
] += enercorr
;
583 enerd
->term
[F_PRES
] += prescorr
;
584 enerd
->term
[F_DVDL_VDW
] += dvdlcorr
;
586 epot
= enerd
->term
[F_EPOT
];
587 bEnergyOutOfBounds
= FALSE
;
589 /* With SSE the energy can overflow, check for this */
590 if (gmx_mm_check_and_reset_overflow())
594 fprintf(debug
,"Found an SSE overflow, assuming the energy is out of bounds\n");
596 bEnergyOutOfBounds
= TRUE
;
599 /* If the compiler doesn't optimize this check away
600 * we catch the NAN energies.
601 * The epot>GMX_REAL_MAX check catches inf values,
602 * which should nicely result in embU=0 through the exp below,
603 * but it does not hurt to check anyhow.
605 /* Non-bonded Interaction usually diverge at r=0.
606 * With tabulated interaction functions the first few entries
607 * should be capped in a consistent fashion between
608 * repulsion, dispersion and Coulomb to avoid accidental
609 * negative values in the total energy.
610 * The table generation code in tables.c does this.
611 * With user tbales the user should take care of this.
613 if (epot
!= epot
|| epot
> GMX_REAL_MAX
)
615 bEnergyOutOfBounds
= TRUE
;
617 if (bEnergyOutOfBounds
)
621 fprintf(debug
,"\n time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n",t
,step
,epot
);
627 embU
= exp(-beta
*epot
);
629 /* Determine the weighted energy contributions of each energy group */
631 sum_UgembU
[e
++] += epot
*embU
;
634 for(i
=0; i
<ngid
; i
++)
637 (enerd
->grpp
.ener
[egBHAMSR
][GID(i
,gid_tp
,ngid
)] +
638 enerd
->grpp
.ener
[egBHAMLR
][GID(i
,gid_tp
,ngid
)])*embU
;
643 for(i
=0; i
<ngid
; i
++)
646 (enerd
->grpp
.ener
[egLJSR
][GID(i
,gid_tp
,ngid
)] +
647 enerd
->grpp
.ener
[egLJLR
][GID(i
,gid_tp
,ngid
)])*embU
;
652 sum_UgembU
[e
++] += enerd
->term
[F_DISPCORR
]*embU
;
656 for(i
=0; i
<ngid
; i
++)
659 (enerd
->grpp
.ener
[egCOULSR
][GID(i
,gid_tp
,ngid
)] +
660 enerd
->grpp
.ener
[egCOULLR
][GID(i
,gid_tp
,ngid
)])*embU
;
664 sum_UgembU
[e
++] += enerd
->term
[F_RF_EXCL
]*embU
;
666 if (EEL_FULL(fr
->eeltype
))
668 sum_UgembU
[e
++] += enerd
->term
[F_COUL_RECIP
]*embU
;
673 if (embU
== 0 || beta
*epot
> bU_bin_limit
)
679 i
= (int)((bU_logV_bin_limit
680 - (beta
*epot
- logV
+ refvolshift
))*invbinw
688 realloc_bins(&bin
,&nbin
,i
+10);
695 fprintf(debug
,"TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
696 step
,epot
,x_tp
[XX
],x_tp
[YY
],x_tp
[ZZ
]);
699 if (dump_pdb
&& epot
<= dump_ener
)
701 sprintf(str
,"t%g_step%d.pdb",t
,step
);
702 sprintf(str2
,"t: %f step %d ener: %f",t
,step
,epot
);
703 write_sto_conf_mtop(str
,str2
,top_global
,state
->x
,state
->v
,
704 inputrec
->ePBC
,state
->box
);
711 /* When running in parallel sum the energies over the processes */
712 gmx_sumd(1, &sum_embU
, cr
);
713 gmx_sumd(nener
,sum_UgembU
,cr
);
718 VembU_all
+= V
*sum_embU
/nsteps
;
722 if (bVerbose
|| frame
%10==0 || frame
<10)
724 fprintf(stderr
,"mu %10.3e <mu> %10.3e\n",
725 -log(sum_embU
/nsteps
)/beta
,-log(VembU_all
/V_all
)/beta
);
728 fprintf(fp_tpi
,"%10.3f %12.5e %12.5e %12.5e %12.5e",
730 VembU_all
==0 ? 20/beta
: -log(VembU_all
/V_all
)/beta
,
731 sum_embU
==0 ? 20/beta
: -log(sum_embU
/nsteps
)/beta
,
733 for(e
=0; e
<nener
; e
++)
735 fprintf(fp_tpi
," %12.5e",sum_UgembU
[e
]/nsteps
);
737 fprintf(fp_tpi
,"\n");
741 bNotLastFrame
= read_next_frame(oenv
, status
,&rerun_fr
);
742 } /* End of the loop */
743 runtime_end(runtime
);
749 gmx_fio_fclose(fp_tpi
);
755 fprintf(fplog
," <V> = %12.5e nm^3\n",V_all
/frame
);
756 fprintf(fplog
," <mu> = %12.5e kJ/mol\n",-log(VembU_all
/V_all
)/beta
);
759 /* Write the Boltzmann factor histogram */
762 /* When running in parallel sum the bins over the processes */
765 realloc_bins(&bin
,&nbin
,i
);
766 gmx_sumd(nbin
,bin
,cr
);
770 fp_tpi
= xvgropen(opt2fn("-tpid",nfile
,fnm
),
771 "TPI energy distribution",
772 "\\betaU - log(V/<V>)","count",oenv
);
773 sprintf(str
,"number \\betaU > %g: %9.3e",bU_bin_limit
,bin
[0]);
774 xvgr_subtitle(fp_tpi
,str
,oenv
);
775 xvgr_legend(fp_tpi
,2,(const char **)tpid_leg
,oenv
);
776 for(i
=nbin
-1; i
>0; i
--)
778 bUlogV
= -i
/invbinw
+ bU_logV_bin_limit
- refvolshift
+ log(V_all
/frame
);
779 fprintf(fp_tpi
,"%6.2f %10d %12.5e\n",
782 bin
[i
]*exp(-bUlogV
)*V_all
/VembU_all
);
784 gmx_fio_fclose(fp_tpi
);
790 runtime
->nsteps_done
= frame
*inputrec
->nsteps
;