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"
82 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_X86_64_SSE2) )
83 #if defined(GMX_DOUBLE)
84 #include "gmx_sse2_double.h"
86 #include "gmx_sse2_single.h"
91 static void global_max(t_commrec
*cr
,int *n
)
97 gmx_sumi(cr
->nnodes
,sum
,cr
);
98 for(i
=0; i
<cr
->nnodes
; i
++)
104 static void realloc_bins(double **bin
,int *nbin
,int nbin_new
)
108 if (nbin_new
!= *nbin
) {
109 srenew(*bin
,nbin_new
);
110 for(i
=*nbin
; i
<nbin_new
; i
++)
116 double do_tpi(FILE *fplog
,t_commrec
*cr
,
117 int nfile
, const t_filenm fnm
[],
118 const output_env_t oenv
, gmx_bool bVerbose
,gmx_bool bCompact
,
120 gmx_vsite_t
*vsite
,gmx_constr_t constr
,
122 t_inputrec
*inputrec
,
123 gmx_mtop_t
*top_global
,t_fcdata
*fcd
,
126 t_nrnb
*nrnb
,gmx_wallcycle_t wcycle
,
129 int repl_ex_nst
,int repl_ex_seed
,
130 real cpt_period
,real max_hours
,
131 const char *deviceOptions
,
133 gmx_runtime_t
*runtime
)
135 const char *TPI
="Test Particle Insertion";
137 gmx_groups_t
*groups
;
138 gmx_enerdata_t
*enerd
;
140 real lambda
,t
,temp
,beta
,drmax
,epot
;
141 double embU
,sum_embU
,*sum_UgembU
,V
,V_all
,VembU_all
;
144 gmx_bool bDispCorr
,bCharge
,bRFExcl
,bNotLastFrame
,bStateChanged
,bNS
,bOurStep
;
145 tensor force_vir
,shake_vir
,vir
,pres
;
146 int cg_tp
,a_tp0
,a_tp1
,ngid
,gid_tp
,nener
,e
;
148 rvec mu_tot
,x_init
,dx
,x_tp
;
149 int nnodes
,frame
,nsteps
,step
;
153 char *ptr
,*dump_pdb
,**leg
,str
[STRLEN
],str2
[STRLEN
];
154 double dbl
,dump_ener
;
157 real
*mass_cavity
=NULL
,mass_tot
;
159 double invbinw
,*bin
,refvolshift
,logV
,bUlogV
;
160 real dvdl
,prescorr
,enercorr
,dvdlcorr
;
161 gmx_bool bEnergyOutOfBounds
;
162 const char *tpid_leg
[2]={"direct","reweighted"};
164 /* Since there is no upper limit to the insertion energies,
165 * we need to set an upper limit for the distribution output.
167 real bU_bin_limit
= 50;
168 real bU_logV_bin_limit
= bU_bin_limit
+ 10;
172 top
= gmx_mtop_generate_local_top(top_global
,inputrec
);
174 groups
= &top_global
->groups
;
176 bCavity
= (inputrec
->eI
== eiTPIC
);
178 ptr
= getenv("GMX_TPIC_MASSES");
182 /* Read (multiple) masses from env var GMX_TPIC_MASSES,
183 * The center of mass of the last atoms is then used for TPIC.
186 while (sscanf(ptr
,"%lf%n",&dbl
,&i
) > 0) {
187 srenew(mass_cavity
,nat_cavity
+1);
188 mass_cavity
[nat_cavity
] = dbl
;
189 fprintf(fplog
,"mass[%d] = %f\n",
190 nat_cavity
+1,mass_cavity
[nat_cavity
]);
195 gmx_fatal(FARGS
,"Found %d masses in GMX_TPIC_MASSES",nat_cavity
);
200 init_em(fplog,TPI,inputrec,&lambda,nrnb,mu_tot,
201 state->box,fr,mdatoms,top,cr,nfile,fnm,NULL,NULL);*/
202 /* We never need full pbc for TPI */
204 /* Determine the temperature for the Boltzmann weighting */
205 temp
= inputrec
->opts
.ref_t
[0];
207 for(i
=1; (i
<inputrec
->opts
.ngtc
); i
++) {
208 if (inputrec
->opts
.ref_t
[i
] != temp
) {
209 fprintf(fplog
,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
210 fprintf(stderr
,"\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
214 "\n The temperature for test particle insertion is %.3f K\n\n",
217 beta
= 1.0/(BOLTZ
*temp
);
219 /* Number of insertions per frame */
220 nsteps
= inputrec
->nsteps
;
222 /* Use the same neighborlist with more insertions points
223 * in a sphere of radius drmax around the initial point
225 /* This should be a proper mdp parameter */
226 drmax
= inputrec
->rtpi
;
228 /* An environment variable can be set to dump all configurations
229 * to pdb with an insertion energy <= this value.
231 dump_pdb
= getenv("GMX_TPI_DUMP");
234 sscanf(dump_pdb
,"%lf",&dump_ener
);
236 atoms2md(top_global
,inputrec
,0,NULL
,0,top_global
->natoms
,mdatoms
);
237 update_mdatoms(mdatoms
,inputrec
->init_lambda
);
240 init_enerdata(groups
->grps
[egcENER
].nr
,inputrec
->n_flambda
,enerd
);
241 snew(f
,top_global
->natoms
);
243 /* Print to log file */
244 runtime_start(runtime
);
245 print_date_and_time(fplog
,cr
->nodeid
,
246 "Started Test Particle Insertion",runtime
);
247 wallcycle_start(wcycle
,ewcRUN
);
249 /* The last charge group is the group to be inserted */
250 cg_tp
= top
->cgs
.nr
- 1;
251 a_tp0
= top
->cgs
.index
[cg_tp
];
252 a_tp1
= top
->cgs
.index
[cg_tp
+1];
254 fprintf(debug
,"TPI cg %d, atoms %d-%d\n",cg_tp
,a_tp0
,a_tp1
);
255 if (a_tp1
- a_tp0
> 1 &&
256 (inputrec
->rlist
< inputrec
->rcoulomb
||
257 inputrec
->rlist
< inputrec
->rvdw
))
258 gmx_fatal(FARGS
,"Can not do TPI for multi-atom molecule with a twin-range cut-off");
259 snew(x_mol
,a_tp1
-a_tp0
);
261 bDispCorr
= (inputrec
->eDispCorr
!= edispcNO
);
263 for(i
=a_tp0
; i
<a_tp1
; i
++) {
264 /* Copy the coordinates of the molecule to be insterted */
265 copy_rvec(state
->x
[i
],x_mol
[i
-a_tp0
]);
266 /* Check if we need to print electrostatic energies */
267 bCharge
|= (mdatoms
->chargeA
[i
] != 0 ||
268 (mdatoms
->chargeB
&& mdatoms
->chargeB
[i
] != 0));
270 bRFExcl
= (bCharge
&& EEL_RF(fr
->eeltype
) && fr
->eeltype
!=eelRF_NEC
);
272 calc_cgcm(fplog
,cg_tp
,cg_tp
+1,&(top
->cgs
),state
->x
,fr
->cg_cm
);
274 if (norm(fr
->cg_cm
[cg_tp
]) > 0.5*inputrec
->rlist
&& fplog
) {
275 fprintf(fplog
, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
276 fprintf(stderr
,"WARNING: Your TPI molecule is not centered at 0,0,0\n");
279 /* Center the molecule to be inserted at zero */
280 for(i
=0; i
<a_tp1
-a_tp0
; i
++)
281 rvec_dec(x_mol
[i
],fr
->cg_cm
[cg_tp
]);
285 fprintf(fplog
,"\nWill insert %d atoms %s partial charges\n",
286 a_tp1
-a_tp0
,bCharge
? "with" : "without");
288 fprintf(fplog
,"\nWill insert %d times in each frame of %s\n",
289 nsteps
,opt2fn("-rerun",nfile
,fnm
));
294 if (inputrec
->nstlist
> 1)
296 if (drmax
==0 && a_tp1
-a_tp0
==1)
298 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
);
302 fprintf(fplog
,"Will use the same neighborlist for %d insertions in a sphere of radius %f\n",inputrec
->nstlist
,drmax
);
310 fprintf(fplog
,"Will insert randomly in a sphere of radius %f around the center of the cavity\n",drmax
);
314 ngid
= groups
->grps
[egcENER
].nr
;
315 gid_tp
= GET_CGINFO_GID(fr
->cginfo
[cg_tp
]);
323 if (EEL_FULL(fr
->eeltype
))
326 snew(sum_UgembU
,nener
);
328 /* Initialize random generator */
329 tpi_rand
= gmx_rng_init(inputrec
->ld_seed
);
332 fp_tpi
= xvgropen(opt2fn("-tpi",nfile
,fnm
),
333 "TPI energies","Time (ps)",
334 "(kJ mol\\S-1\\N) / (nm\\S3\\N)",oenv
);
335 xvgr_subtitle(fp_tpi
,"f. are averages over one frame",oenv
);
338 sprintf(str
,"-kT log(<Ve\\S-\\betaU\\N>/<V>)");
339 leg
[e
++] = strdup(str
);
340 sprintf(str
,"f. -kT log<e\\S-\\betaU\\N>");
341 leg
[e
++] = strdup(str
);
342 sprintf(str
,"f. <e\\S-\\betaU\\N>");
343 leg
[e
++] = strdup(str
);
345 leg
[e
++] = strdup(str
);
346 sprintf(str
,"f. <Ue\\S-\\betaU\\N>");
347 leg
[e
++] = strdup(str
);
348 for(i
=0; i
<ngid
; i
++) {
349 sprintf(str
,"f. <U\\sVdW %s\\Ne\\S-\\betaU\\N>",
350 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
351 leg
[e
++] = strdup(str
);
354 sprintf(str
,"f. <U\\sdisp c\\Ne\\S-\\betaU\\N>");
355 leg
[e
++] = strdup(str
);
358 for(i
=0; i
<ngid
; i
++) {
359 sprintf(str
,"f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
360 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
361 leg
[e
++] = strdup(str
);
364 sprintf(str
,"f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
365 leg
[e
++] = strdup(str
);
367 if (EEL_FULL(fr
->eeltype
)) {
368 sprintf(str
,"f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
369 leg
[e
++] = strdup(str
);
372 xvgr_legend(fp_tpi
,4+nener
,(const char**)leg
,oenv
);
373 for(i
=0; i
<4+nener
; i
++)
385 bNotLastFrame
= read_first_frame(oenv
,&status
,opt2fn("-rerun",nfile
,fnm
),
386 &rerun_fr
,TRX_NEED_X
);
389 if (rerun_fr
.natoms
- (bCavity
? nat_cavity
: 0) !=
390 mdatoms
->nr
- (a_tp1
- a_tp0
))
391 gmx_fatal(FARGS
,"Number of atoms in trajectory (%d)%s "
392 "is not equal the number in the run input file (%d) "
393 "minus the number of atoms to insert (%d)\n",
394 rerun_fr
.natoms
,bCavity
? " minus one" : "",
395 mdatoms
->nr
,a_tp1
-a_tp0
);
397 refvolshift
= log(det(rerun_fr
.box
));
399 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_X86_64_SSE2) )
400 /* Make sure we don't detect SSE overflow generated before this point */
401 gmx_mm_check_and_reset_overflow();
404 while (bNotLastFrame
)
406 lambda
= rerun_fr
.lambda
;
410 for(e
=0; e
<nener
; e
++)
415 /* Copy the coordinates from the input trajectory */
416 for(i
=0; i
<rerun_fr
.natoms
; i
++)
418 copy_rvec(rerun_fr
.x
[i
],state
->x
[i
]);
421 V
= det(rerun_fr
.box
);
424 bStateChanged
= TRUE
;
426 for(step
=0; step
<nsteps
; step
++)
428 /* In parallel all nodes generate all random configurations.
429 * In that way the result is identical to a single cpu tpi run.
433 /* Random insertion in the whole volume */
434 bNS
= (step
% inputrec
->nstlist
== 0);
437 /* Generate a random position in the box */
438 x_init
[XX
] = gmx_rng_uniform_real(tpi_rand
)*state
->box
[XX
][XX
];
439 x_init
[YY
] = gmx_rng_uniform_real(tpi_rand
)*state
->box
[YY
][YY
];
440 x_init
[ZZ
] = gmx_rng_uniform_real(tpi_rand
)*state
->box
[ZZ
][ZZ
];
442 if (inputrec
->nstlist
== 1)
444 copy_rvec(x_init
,x_tp
);
448 /* Generate coordinates within |dx|=drmax of x_init */
451 dx
[XX
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
452 dx
[YY
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
453 dx
[ZZ
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
455 while (norm2(dx
) > drmax
*drmax
);
456 rvec_add(x_init
,dx
,x_tp
);
461 /* Random insertion around a cavity location
462 * given by the last coordinate of the trajectory.
468 /* Copy the location of the cavity */
469 copy_rvec(rerun_fr
.x
[rerun_fr
.natoms
-1],x_init
);
473 /* Determine the center of mass of the last molecule */
476 for(i
=0; i
<nat_cavity
; i
++)
481 mass_cavity
[i
]*rerun_fr
.x
[rerun_fr
.natoms
-nat_cavity
+i
][d
];
483 mass_tot
+= mass_cavity
[i
];
487 x_init
[d
] /= mass_tot
;
491 /* Generate coordinates within |dx|=drmax of x_init */
494 dx
[XX
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
495 dx
[YY
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
496 dx
[ZZ
] = (2*gmx_rng_uniform_real(tpi_rand
) - 1)*drmax
;
498 while (norm2(dx
) > drmax
*drmax
);
499 rvec_add(x_init
,dx
,x_tp
);
502 if (a_tp1
- a_tp0
== 1)
504 /* Insert a single atom, just copy the insertion location */
505 copy_rvec(x_tp
,state
->x
[a_tp0
]);
509 /* Copy the coordinates from the top file */
510 for(i
=a_tp0
; i
<a_tp1
; i
++)
512 copy_rvec(x_mol
[i
-a_tp0
],state
->x
[i
]);
514 /* Rotate the molecule randomly */
515 rotate_conf(a_tp1
-a_tp0
,state
->x
+a_tp0
,NULL
,
516 2*M_PI
*gmx_rng_uniform_real(tpi_rand
),
517 2*M_PI
*gmx_rng_uniform_real(tpi_rand
),
518 2*M_PI
*gmx_rng_uniform_real(tpi_rand
));
519 /* Shift to the insertion location */
520 for(i
=a_tp0
; i
<a_tp1
; i
++)
522 rvec_inc(state
->x
[i
],x_tp
);
526 /* Check if this insertion belongs to this node */
530 switch (inputrec
->eI
)
533 bOurStep
= ((step
/ inputrec
->nstlist
) % nnodes
== cr
->nodeid
);
536 bOurStep
= (step
% nnodes
== cr
->nodeid
);
539 gmx_fatal(FARGS
,"Unknown integrator %s",ei_names
[inputrec
->eI
]);
544 /* Clear some matrix variables */
545 clear_mat(force_vir
);
546 clear_mat(shake_vir
);
550 /* Set the charge group center of mass of the test particle */
551 copy_rvec(x_init
,fr
->cg_cm
[top
->cgs
.nr
-1]);
553 /* Calc energy (no forces) on new positions.
554 * Since we only need the intermolecular energy
555 * and the RF exclusion terms of the inserted molecule occur
556 * within a single charge group we can pass NULL for the graph.
557 * This also avoids shifts that would move charge groups
560 * Some checks above ensure than we can not have
561 * twin-range interactions together with nstlist > 1,
562 * therefore we do not need to remember the LR energies.
564 /* Make do_force do a single node force calculation */
566 do_force(fplog
,cr
,inputrec
,
567 step
,nrnb
,wcycle
,top
,top_global
,&top_global
->groups
,
568 rerun_fr
.box
,state
->x
,&state
->hist
,
569 f
,force_vir
,mdatoms
,enerd
,fcd
,
570 lambda
,NULL
,fr
,NULL
,mu_tot
,t
,NULL
,NULL
,FALSE
,
571 GMX_FORCE_NONBONDED
|
572 (bNS
? GMX_FORCE_NS
| GMX_FORCE_DOLR
: 0) |
573 (bStateChanged
? GMX_FORCE_STATECHANGED
: 0));
575 bStateChanged
= FALSE
;
578 /* Calculate long range corrections to pressure and energy */
579 calc_dispcorr(fplog
,inputrec
,fr
,step
,top_global
->natoms
,rerun_fr
.box
,
580 lambda
,pres
,vir
,&prescorr
,&enercorr
,&dvdlcorr
);
581 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
582 enerd
->term
[F_DISPCORR
] = enercorr
;
583 enerd
->term
[F_EPOT
] += enercorr
;
584 enerd
->term
[F_PRES
] += prescorr
;
585 enerd
->term
[F_DVDL
] += dvdlcorr
;
587 epot
= enerd
->term
[F_EPOT
];
588 bEnergyOutOfBounds
= FALSE
;
589 #if ( defined(GMX_IA32_SSE) || defined(GMX_X86_64_SSE) || defined(GMX_X86_64_SSE2) )
590 /* With SSE the energy can overflow, check for this */
591 if (gmx_mm_check_and_reset_overflow())
595 fprintf(debug
,"Found an SSE overflow, assuming the energy is out of bounds\n");
597 bEnergyOutOfBounds
= TRUE
;
600 /* If the compiler doesn't optimize this check away
601 * we catch the NAN energies.
602 * The epot>GMX_REAL_MAX check catches inf values,
603 * which should nicely result in embU=0 through the exp below,
604 * but it does not hurt to check anyhow.
606 /* Non-bonded Interaction usually diverge at r=0.
607 * With tabulated interaction functions the first few entries
608 * should be capped in a consistent fashion between
609 * repulsion, dispersion and Coulomb to avoid accidental
610 * negative values in the total energy.
611 * The table generation code in tables.c does this.
612 * With user tbales the user should take care of this.
614 if (epot
!= epot
|| epot
> GMX_REAL_MAX
)
616 bEnergyOutOfBounds
= TRUE
;
618 if (bEnergyOutOfBounds
)
622 fprintf(debug
,"\n time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n",t
,step
,epot
);
628 embU
= exp(-beta
*epot
);
630 /* Determine the weighted energy contributions of each energy group */
632 sum_UgembU
[e
++] += epot
*embU
;
635 for(i
=0; i
<ngid
; i
++)
638 (enerd
->grpp
.ener
[egBHAMSR
][GID(i
,gid_tp
,ngid
)] +
639 enerd
->grpp
.ener
[egBHAMLR
][GID(i
,gid_tp
,ngid
)])*embU
;
644 for(i
=0; i
<ngid
; i
++)
647 (enerd
->grpp
.ener
[egLJSR
][GID(i
,gid_tp
,ngid
)] +
648 enerd
->grpp
.ener
[egLJLR
][GID(i
,gid_tp
,ngid
)])*embU
;
653 sum_UgembU
[e
++] += enerd
->term
[F_DISPCORR
]*embU
;
657 for(i
=0; i
<ngid
; i
++)
660 (enerd
->grpp
.ener
[egCOULSR
][GID(i
,gid_tp
,ngid
)] +
661 enerd
->grpp
.ener
[egCOULLR
][GID(i
,gid_tp
,ngid
)])*embU
;
665 sum_UgembU
[e
++] += enerd
->term
[F_RF_EXCL
]*embU
;
667 if (EEL_FULL(fr
->eeltype
))
669 sum_UgembU
[e
++] += enerd
->term
[F_COUL_RECIP
]*embU
;
674 if (embU
== 0 || beta
*epot
> bU_bin_limit
)
680 i
= (int)((bU_logV_bin_limit
681 - (beta
*epot
- logV
+ refvolshift
))*invbinw
689 realloc_bins(&bin
,&nbin
,i
+10);
696 fprintf(debug
,"TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
697 step
,epot
,x_tp
[XX
],x_tp
[YY
],x_tp
[ZZ
]);
700 if (dump_pdb
&& epot
<= dump_ener
)
702 sprintf(str
,"t%g_step%d.pdb",t
,step
);
703 sprintf(str2
,"t: %f step %d ener: %f",t
,step
,epot
);
704 write_sto_conf_mtop(str
,str2
,top_global
,state
->x
,state
->v
,
705 inputrec
->ePBC
,state
->box
);
712 /* When running in parallel sum the energies over the processes */
713 gmx_sumd(1, &sum_embU
, cr
);
714 gmx_sumd(nener
,sum_UgembU
,cr
);
719 VembU_all
+= V
*sum_embU
/nsteps
;
723 if (bVerbose
|| frame
%10==0 || frame
<10)
725 fprintf(stderr
,"mu %10.3e <mu> %10.3e\n",
726 -log(sum_embU
/nsteps
)/beta
,-log(VembU_all
/V_all
)/beta
);
729 fprintf(fp_tpi
,"%10.3f %12.5e %12.5e %12.5e %12.5e",
731 VembU_all
==0 ? 20/beta
: -log(VembU_all
/V_all
)/beta
,
732 sum_embU
==0 ? 20/beta
: -log(sum_embU
/nsteps
)/beta
,
734 for(e
=0; e
<nener
; e
++)
736 fprintf(fp_tpi
," %12.5e",sum_UgembU
[e
]/nsteps
);
738 fprintf(fp_tpi
,"\n");
742 bNotLastFrame
= read_next_frame(oenv
, status
,&rerun_fr
);
743 } /* End of the loop */
744 runtime_end(runtime
);
750 gmx_fio_fclose(fp_tpi
);
756 fprintf(fplog
," <V> = %12.5e nm^3\n",V_all
/frame
);
757 fprintf(fplog
," <mu> = %12.5e kJ/mol\n",-log(VembU_all
/V_all
)/beta
);
760 /* Write the Boltzmann factor histogram */
763 /* When running in parallel sum the bins over the processes */
766 realloc_bins(&bin
,&nbin
,i
);
767 gmx_sumd(nbin
,bin
,cr
);
771 fp_tpi
= xvgropen(opt2fn("-tpid",nfile
,fnm
),
772 "TPI energy distribution",
773 "\\betaU - log(V/<V>)","count",oenv
);
774 sprintf(str
,"number \\betaU > %g: %9.3e",bU_bin_limit
,bin
[0]);
775 xvgr_subtitle(fp_tpi
,str
,oenv
);
776 xvgr_legend(fp_tpi
,2,(const char **)tpid_leg
,oenv
);
777 for(i
=nbin
-1; i
>0; i
--)
779 bUlogV
= -i
/invbinw
+ bU_logV_bin_limit
- refvolshift
+ log(V_all
/frame
);
780 fprintf(fp_tpi
,"%6.2f %10d %12.5e\n",
783 bin
[i
]*exp(-bUlogV
)*V_all
/VembU_all
);
785 gmx_fio_fclose(fp_tpi
);
791 runtime
->nsteps_done
= frame
*inputrec
->nsteps
;