2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
45 #include "gromacs/utility/smalloc.h"
47 #include "chargegroup.h"
51 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/math/vec.h"
63 #include "gromacs/random/random.h"
64 #include "gromacs/math/units.h"
65 #include "gromacs/fileio/xvgr.h"
68 #include "gromacs/topology/mtop_util.h"
70 #include "gromacs/gmxlib/conformation-utilities.h"
72 #include "gromacs/legacyheaders/types/commrec.h"
73 #include "gromacs/fileio/confio.h"
74 #include "gromacs/fileio/gmxfio.h"
75 #include "gromacs/fileio/trxio.h"
76 #include "gromacs/timing/wallcycle.h"
77 #include "gromacs/timing/walltime_accounting.h"
79 static void global_max(t_commrec
*cr
, int *n
)
83 snew(sum
, cr
->nnodes
);
85 gmx_sumi(cr
->nnodes
, sum
, cr
);
86 for (i
= 0; i
< cr
->nnodes
; i
++)
94 static void realloc_bins(double **bin
, int *nbin
, int nbin_new
)
98 if (nbin_new
!= *nbin
)
100 srenew(*bin
, nbin_new
);
101 for (i
= *nbin
; i
< nbin_new
; i
++)
109 double do_tpi(FILE *fplog
, t_commrec
*cr
,
110 int nfile
, const t_filenm fnm
[],
111 const output_env_t oenv
, gmx_bool bVerbose
, gmx_bool gmx_unused bCompact
,
112 int gmx_unused nstglobalcomm
,
113 gmx_vsite_t gmx_unused
*vsite
, gmx_constr_t gmx_unused constr
,
114 int gmx_unused stepout
,
115 t_inputrec
*inputrec
,
116 gmx_mtop_t
*top_global
, t_fcdata
*fcd
,
119 t_nrnb
*nrnb
, gmx_wallcycle_t wcycle
,
120 gmx_edsam_t gmx_unused ed
,
122 int gmx_unused repl_ex_nst
, int gmx_unused repl_ex_nex
, int gmx_unused repl_ex_seed
,
123 gmx_membed_t gmx_unused membed
,
124 real gmx_unused cpt_period
, real gmx_unused max_hours
,
125 const char gmx_unused
*deviceOptions
,
126 int gmx_unused imdport
,
127 unsigned long gmx_unused Flags
,
128 gmx_walltime_accounting_t walltime_accounting
)
130 const char *TPI
= "Test Particle Insertion";
132 gmx_groups_t
*groups
;
133 gmx_enerdata_t
*enerd
;
135 real lambda
, t
, temp
, beta
, drmax
, epot
;
136 double embU
, sum_embU
, *sum_UgembU
, V
, V_all
, VembU_all
;
139 gmx_bool bDispCorr
, bCharge
, bRFExcl
, bNotLastFrame
, bStateChanged
, bNS
;
140 tensor force_vir
, shake_vir
, vir
, pres
;
141 int cg_tp
, a_tp0
, a_tp1
, ngid
, gid_tp
, nener
, e
;
143 rvec mu_tot
, x_init
, dx
, x_tp
;
145 gmx_int64_t frame_step_prev
, frame_step
;
146 gmx_int64_t nsteps
, stepblocksize
= 0, step
;
147 gmx_int64_t rnd_count_stride
, rnd_count
;
152 char *ptr
, *dump_pdb
, **leg
, str
[STRLEN
], str2
[STRLEN
];
153 double dbl
, dump_ener
;
155 int nat_cavity
= 0, d
;
156 real
*mass_cavity
= NULL
, mass_tot
;
158 double invbinw
, *bin
, refvolshift
, logV
, bUlogV
;
159 real dvdl
, prescorr
, enercorr
, dvdlcorr
;
160 gmx_bool bEnergyOutOfBounds
;
161 const char *tpid_leg
[2] = {"direct", "reweighted"};
163 /* Since there is no upper limit to the insertion energies,
164 * we need to set an upper limit for the distribution output.
166 real bU_bin_limit
= 50;
167 real bU_logV_bin_limit
= bU_bin_limit
+ 10;
171 top
= gmx_mtop_generate_local_top(top_global
, inputrec
);
173 groups
= &top_global
->groups
;
175 bCavity
= (inputrec
->eI
== eiTPIC
);
178 ptr
= getenv("GMX_TPIC_MASSES");
185 /* Read (multiple) masses from env var GMX_TPIC_MASSES,
186 * The center of mass of the last atoms is then used for TPIC.
189 while (sscanf(ptr
, "%lf%n", &dbl
, &i
) > 0)
191 srenew(mass_cavity
, nat_cavity
+1);
192 mass_cavity
[nat_cavity
] = dbl
;
193 fprintf(fplog
, "mass[%d] = %f\n",
194 nat_cavity
+1, mass_cavity
[nat_cavity
]);
200 gmx_fatal(FARGS
, "Found %d masses in GMX_TPIC_MASSES", nat_cavity
);
206 init_em(fplog,TPI,inputrec,&lambda,nrnb,mu_tot,
207 state->box,fr,mdatoms,top,cr,nfile,fnm,NULL,NULL);*/
208 /* We never need full pbc for TPI */
210 /* Determine the temperature for the Boltzmann weighting */
211 temp
= inputrec
->opts
.ref_t
[0];
214 for (i
= 1; (i
< inputrec
->opts
.ngtc
); i
++)
216 if (inputrec
->opts
.ref_t
[i
] != temp
)
218 fprintf(fplog
, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
219 fprintf(stderr
, "\nWARNING: The temperatures of the different temperature coupling groups are not identical\n\n");
223 "\n The temperature for test particle insertion is %.3f K\n\n",
226 beta
= 1.0/(BOLTZ
*temp
);
228 /* Number of insertions per frame */
229 nsteps
= inputrec
->nsteps
;
231 /* Use the same neighborlist with more insertions points
232 * in a sphere of radius drmax around the initial point
234 /* This should be a proper mdp parameter */
235 drmax
= inputrec
->rtpi
;
237 /* An environment variable can be set to dump all configurations
238 * to pdb with an insertion energy <= this value.
240 dump_pdb
= getenv("GMX_TPI_DUMP");
244 sscanf(dump_pdb
, "%lf", &dump_ener
);
247 atoms2md(top_global
, inputrec
, 0, NULL
, top_global
->natoms
, mdatoms
);
248 update_mdatoms(mdatoms
, inputrec
->fepvals
->init_lambda
);
251 init_enerdata(groups
->grps
[egcENER
].nr
, inputrec
->fepvals
->n_lambda
, enerd
);
252 snew(f
, top_global
->natoms
);
254 /* Print to log file */
255 walltime_accounting_start(walltime_accounting
);
256 wallcycle_start(wcycle
, ewcRUN
);
257 print_start(fplog
, cr
, walltime_accounting
, "Test Particle Insertion");
259 /* The last charge group is the group to be inserted */
260 cg_tp
= top
->cgs
.nr
- 1;
261 a_tp0
= top
->cgs
.index
[cg_tp
];
262 a_tp1
= top
->cgs
.index
[cg_tp
+1];
265 fprintf(debug
, "TPI cg %d, atoms %d-%d\n", cg_tp
, a_tp0
, a_tp1
);
267 if (a_tp1
- a_tp0
> 1 &&
268 (inputrec
->rlist
< inputrec
->rcoulomb
||
269 inputrec
->rlist
< inputrec
->rvdw
))
271 gmx_fatal(FARGS
, "Can not do TPI for multi-atom molecule with a twin-range cut-off");
273 snew(x_mol
, a_tp1
-a_tp0
);
275 bDispCorr
= (inputrec
->eDispCorr
!= edispcNO
);
277 for (i
= a_tp0
; i
< a_tp1
; i
++)
279 /* Copy the coordinates of the molecule to be insterted */
280 copy_rvec(state
->x
[i
], x_mol
[i
-a_tp0
]);
281 /* Check if we need to print electrostatic energies */
282 bCharge
|= (mdatoms
->chargeA
[i
] != 0 ||
283 (mdatoms
->chargeB
&& mdatoms
->chargeB
[i
] != 0));
285 bRFExcl
= (bCharge
&& EEL_RF(fr
->eeltype
) && fr
->eeltype
!= eelRF_NEC
);
287 calc_cgcm(fplog
, cg_tp
, cg_tp
+1, &(top
->cgs
), state
->x
, fr
->cg_cm
);
290 if (norm(fr
->cg_cm
[cg_tp
]) > 0.5*inputrec
->rlist
&& fplog
)
292 fprintf(fplog
, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
293 fprintf(stderr
, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
298 /* Center the molecule to be inserted at zero */
299 for (i
= 0; i
< a_tp1
-a_tp0
; i
++)
301 rvec_dec(x_mol
[i
], fr
->cg_cm
[cg_tp
]);
307 fprintf(fplog
, "\nWill insert %d atoms %s partial charges\n",
308 a_tp1
-a_tp0
, bCharge
? "with" : "without");
310 fprintf(fplog
, "\nWill insert %d times in each frame of %s\n",
311 (int)nsteps
, opt2fn("-rerun", nfile
, fnm
));
316 if (inputrec
->nstlist
> 1)
318 if (drmax
== 0 && a_tp1
-a_tp0
== 1)
320 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
);
324 fprintf(fplog
, "Will use the same neighborlist for %d insertions in a sphere of radius %f\n", inputrec
->nstlist
, drmax
);
332 fprintf(fplog
, "Will insert randomly in a sphere of radius %f around the center of the cavity\n", drmax
);
336 ngid
= groups
->grps
[egcENER
].nr
;
337 gid_tp
= GET_CGINFO_GID(fr
->cginfo
[cg_tp
]);
350 if (EEL_FULL(fr
->eeltype
))
355 snew(sum_UgembU
, nener
);
357 /* Copy the random seed set by the user */
358 seed
= inputrec
->ld_seed
;
359 /* We use the frame step number as one random counter.
360 * The second counter use the insertion (step) count. But we
361 * need multiple random numbers per insertion. This number is
362 * not fixed, since we generate random locations in a sphere
363 * by putting locations in a cube and some of these fail.
364 * A count of 20 is already extremely unlikely, so 10000 is
365 * a safe margin for random numbers per insertion.
367 rnd_count_stride
= 10000;
371 fp_tpi
= xvgropen(opt2fn("-tpi", nfile
, fnm
),
372 "TPI energies", "Time (ps)",
373 "(kJ mol\\S-1\\N) / (nm\\S3\\N)", oenv
);
374 xvgr_subtitle(fp_tpi
, "f. are averages over one frame", oenv
);
377 sprintf(str
, "-kT log(<Ve\\S-\\betaU\\N>/<V>)");
378 leg
[e
++] = strdup(str
);
379 sprintf(str
, "f. -kT log<e\\S-\\betaU\\N>");
380 leg
[e
++] = strdup(str
);
381 sprintf(str
, "f. <e\\S-\\betaU\\N>");
382 leg
[e
++] = strdup(str
);
383 sprintf(str
, "f. V");
384 leg
[e
++] = strdup(str
);
385 sprintf(str
, "f. <Ue\\S-\\betaU\\N>");
386 leg
[e
++] = strdup(str
);
387 for (i
= 0; i
< ngid
; i
++)
389 sprintf(str
, "f. <U\\sVdW %s\\Ne\\S-\\betaU\\N>",
390 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
391 leg
[e
++] = strdup(str
);
395 sprintf(str
, "f. <U\\sdisp c\\Ne\\S-\\betaU\\N>");
396 leg
[e
++] = strdup(str
);
400 for (i
= 0; i
< ngid
; i
++)
402 sprintf(str
, "f. <U\\sCoul %s\\Ne\\S-\\betaU\\N>",
403 *(groups
->grpname
[groups
->grps
[egcENER
].nm_ind
[i
]]));
404 leg
[e
++] = strdup(str
);
408 sprintf(str
, "f. <U\\sRF excl\\Ne\\S-\\betaU\\N>");
409 leg
[e
++] = strdup(str
);
411 if (EEL_FULL(fr
->eeltype
))
413 sprintf(str
, "f. <U\\sCoul recip\\Ne\\S-\\betaU\\N>");
414 leg
[e
++] = strdup(str
);
417 xvgr_legend(fp_tpi
, 4+nener
, (const char**)leg
, oenv
);
418 for (i
= 0; i
< 4+nener
; i
++)
432 /* Avoid frame step numbers <= -1 */
433 frame_step_prev
= -1;
435 bNotLastFrame
= read_first_frame(oenv
, &status
, opt2fn("-rerun", nfile
, fnm
),
436 &rerun_fr
, TRX_NEED_X
);
439 if (rerun_fr
.natoms
- (bCavity
? nat_cavity
: 0) !=
440 mdatoms
->nr
- (a_tp1
- a_tp0
))
442 gmx_fatal(FARGS
, "Number of atoms in trajectory (%d)%s "
443 "is not equal the number in the run input file (%d) "
444 "minus the number of atoms to insert (%d)\n",
445 rerun_fr
.natoms
, bCavity
? " minus one" : "",
446 mdatoms
->nr
, a_tp1
-a_tp0
);
449 refvolshift
= log(det(rerun_fr
.box
));
451 switch (inputrec
->eI
)
454 stepblocksize
= inputrec
->nstlist
;
460 gmx_fatal(FARGS
, "Unknown integrator %s", ei_names
[inputrec
->eI
]);
464 /* Make sure we don't detect SIMD overflow generated before this point */
465 gmx_simd_check_and_reset_overflow();
468 while (bNotLastFrame
)
470 frame_step
= rerun_fr
.step
;
471 if (frame_step
<= frame_step_prev
)
473 /* We don't have step number in the trajectory file,
474 * or we have constant or decreasing step numbers.
475 * Ensure we have increasing step numbers, since we use
476 * the step numbers as a counter for random numbers.
478 frame_step
= frame_step_prev
+ 1;
480 frame_step_prev
= frame_step
;
482 lambda
= rerun_fr
.lambda
;
486 for (e
= 0; e
< nener
; e
++)
491 /* Copy the coordinates from the input trajectory */
492 for (i
= 0; i
< rerun_fr
.natoms
; i
++)
494 copy_rvec(rerun_fr
.x
[i
], state
->x
[i
]);
496 copy_mat(rerun_fr
.box
, state
->box
);
501 bStateChanged
= TRUE
;
504 step
= cr
->nodeid
*stepblocksize
;
505 while (step
< nsteps
)
507 /* Initialize the second counter for random numbers using
508 * the insertion step index. This ensures that we get
509 * the same random numbers independently of how many
510 * MPI ranks we use. Also for the same seed, we get
511 * the same initial random sequence for different nsteps.
513 rnd_count
= step
*rnd_count_stride
;
517 /* Random insertion in the whole volume */
518 bNS
= (step
% inputrec
->nstlist
== 0);
521 /* Generate a random position in the box */
522 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
);
523 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
+2);
524 for (d
= 0; d
< DIM
; d
++)
526 x_init
[d
] = rnd
[d
]*state
->box
[d
][d
];
529 if (inputrec
->nstlist
== 1)
531 copy_rvec(x_init
, x_tp
);
535 /* Generate coordinates within |dx|=drmax of x_init */
538 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
);
539 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
+2);
540 for (d
= 0; d
< DIM
; d
++)
542 dx
[d
] = (2*rnd
[d
] - 1)*drmax
;
545 while (norm2(dx
) > drmax
*drmax
);
546 rvec_add(x_init
, dx
, x_tp
);
551 /* Random insertion around a cavity location
552 * given by the last coordinate of the trajectory.
558 /* Copy the location of the cavity */
559 copy_rvec(rerun_fr
.x
[rerun_fr
.natoms
-1], x_init
);
563 /* Determine the center of mass of the last molecule */
566 for (i
= 0; i
< nat_cavity
; i
++)
568 for (d
= 0; d
< DIM
; d
++)
571 mass_cavity
[i
]*rerun_fr
.x
[rerun_fr
.natoms
-nat_cavity
+i
][d
];
573 mass_tot
+= mass_cavity
[i
];
575 for (d
= 0; d
< DIM
; d
++)
577 x_init
[d
] /= mass_tot
;
581 /* Generate coordinates within |dx|=drmax of x_init */
584 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
);
585 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
+2);
586 for (d
= 0; d
< DIM
; d
++)
588 dx
[d
] = (2*rnd
[d
] - 1)*drmax
;
591 while (norm2(dx
) > drmax
*drmax
);
592 rvec_add(x_init
, dx
, x_tp
);
595 if (a_tp1
- a_tp0
== 1)
597 /* Insert a single atom, just copy the insertion location */
598 copy_rvec(x_tp
, state
->x
[a_tp0
]);
602 /* Copy the coordinates from the top file */
603 for (i
= a_tp0
; i
< a_tp1
; i
++)
605 copy_rvec(x_mol
[i
-a_tp0
], state
->x
[i
]);
607 /* Rotate the molecule randomly */
608 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
);
609 gmx_rng_cycle_2uniform(frame_step
, rnd_count
++, seed
, RND_SEED_TPI
, rnd
+2);
610 rotate_conf(a_tp1
-a_tp0
, state
->x
+a_tp0
, NULL
,
614 /* Shift to the insertion location */
615 for (i
= a_tp0
; i
< a_tp1
; i
++)
617 rvec_inc(state
->x
[i
], x_tp
);
621 /* Clear some matrix variables */
622 clear_mat(force_vir
);
623 clear_mat(shake_vir
);
627 /* Set the charge group center of mass of the test particle */
628 copy_rvec(x_init
, fr
->cg_cm
[top
->cgs
.nr
-1]);
630 /* Calc energy (no forces) on new positions.
631 * Since we only need the intermolecular energy
632 * and the RF exclusion terms of the inserted molecule occur
633 * within a single charge group we can pass NULL for the graph.
634 * This also avoids shifts that would move charge groups
637 * Some checks above ensure than we can not have
638 * twin-range interactions together with nstlist > 1,
639 * therefore we do not need to remember the LR energies.
641 /* Make do_force do a single node force calculation */
643 do_force(fplog
, cr
, inputrec
,
644 step
, nrnb
, wcycle
, top
, &top_global
->groups
,
645 state
->box
, state
->x
, &state
->hist
,
646 f
, force_vir
, mdatoms
, enerd
, fcd
,
648 NULL
, fr
, NULL
, mu_tot
, t
, NULL
, NULL
, FALSE
,
649 GMX_FORCE_NONBONDED
| GMX_FORCE_ENERGY
|
650 (bNS
? GMX_FORCE_DYNAMICBOX
| GMX_FORCE_NS
| GMX_FORCE_DO_LR
: 0) |
651 (bStateChanged
? GMX_FORCE_STATECHANGED
: 0));
653 bStateChanged
= FALSE
;
656 /* Calculate long range corrections to pressure and energy */
657 calc_dispcorr(fplog
, inputrec
, fr
, step
, top_global
->natoms
, state
->box
,
658 lambda
, pres
, vir
, &prescorr
, &enercorr
, &dvdlcorr
);
659 /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
660 enerd
->term
[F_DISPCORR
] = enercorr
;
661 enerd
->term
[F_EPOT
] += enercorr
;
662 enerd
->term
[F_PRES
] += prescorr
;
663 enerd
->term
[F_DVDL_VDW
] += dvdlcorr
;
665 epot
= enerd
->term
[F_EPOT
];
666 bEnergyOutOfBounds
= FALSE
;
667 #ifdef GMX_SIMD_X86_SSE2_OR_HIGHER
668 /* With SSE the energy can overflow, check for this */
669 if (gmx_mm_check_and_reset_overflow())
673 fprintf(debug
, "Found an SSE overflow, assuming the energy is out of bounds\n");
675 bEnergyOutOfBounds
= TRUE
;
678 /* If the compiler doesn't optimize this check away
679 * we catch the NAN energies.
680 * The epot>GMX_REAL_MAX check catches inf values,
681 * which should nicely result in embU=0 through the exp below,
682 * but it does not hurt to check anyhow.
684 /* Non-bonded Interaction usually diverge at r=0.
685 * With tabulated interaction functions the first few entries
686 * should be capped in a consistent fashion between
687 * repulsion, dispersion and Coulomb to avoid accidental
688 * negative values in the total energy.
689 * The table generation code in tables.c does this.
690 * With user tbales the user should take care of this.
692 if (epot
!= epot
|| epot
> GMX_REAL_MAX
)
694 bEnergyOutOfBounds
= TRUE
;
696 if (bEnergyOutOfBounds
)
700 fprintf(debug
, "\n time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n", t
, (int)step
, epot
);
706 embU
= exp(-beta
*epot
);
708 /* Determine the weighted energy contributions of each energy group */
710 sum_UgembU
[e
++] += epot
*embU
;
713 for (i
= 0; i
< ngid
; i
++)
716 (enerd
->grpp
.ener
[egBHAMSR
][GID(i
, gid_tp
, ngid
)] +
717 enerd
->grpp
.ener
[egBHAMLR
][GID(i
, gid_tp
, ngid
)])*embU
;
722 for (i
= 0; i
< ngid
; i
++)
725 (enerd
->grpp
.ener
[egLJSR
][GID(i
, gid_tp
, ngid
)] +
726 enerd
->grpp
.ener
[egLJLR
][GID(i
, gid_tp
, ngid
)])*embU
;
731 sum_UgembU
[e
++] += enerd
->term
[F_DISPCORR
]*embU
;
735 for (i
= 0; i
< ngid
; i
++)
738 (enerd
->grpp
.ener
[egCOULSR
][GID(i
, gid_tp
, ngid
)] +
739 enerd
->grpp
.ener
[egCOULLR
][GID(i
, gid_tp
, ngid
)])*embU
;
743 sum_UgembU
[e
++] += enerd
->term
[F_RF_EXCL
]*embU
;
745 if (EEL_FULL(fr
->eeltype
))
747 sum_UgembU
[e
++] += enerd
->term
[F_COUL_RECIP
]*embU
;
752 if (embU
== 0 || beta
*epot
> bU_bin_limit
)
758 i
= (int)((bU_logV_bin_limit
759 - (beta
*epot
- logV
+ refvolshift
))*invbinw
767 realloc_bins(&bin
, &nbin
, i
+10);
774 fprintf(debug
, "TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
775 (int)step
, epot
, x_tp
[XX
], x_tp
[YY
], x_tp
[ZZ
]);
778 if (dump_pdb
&& epot
<= dump_ener
)
780 sprintf(str
, "t%g_step%d.pdb", t
, (int)step
);
781 sprintf(str2
, "t: %f step %d ener: %f", t
, (int)step
, epot
);
782 write_sto_conf_mtop(str
, str2
, top_global
, state
->x
, state
->v
,
783 inputrec
->ePBC
, state
->box
);
787 if ((step
/stepblocksize
) % cr
->nnodes
!= cr
->nodeid
)
789 /* Skip all steps assigned to the other MPI ranks */
790 step
+= (cr
->nnodes
- 1)*stepblocksize
;
796 /* When running in parallel sum the energies over the processes */
797 gmx_sumd(1, &sum_embU
, cr
);
798 gmx_sumd(nener
, sum_UgembU
, cr
);
803 VembU_all
+= V
*sum_embU
/nsteps
;
807 if (bVerbose
|| frame
%10 == 0 || frame
< 10)
809 fprintf(stderr
, "mu %10.3e <mu> %10.3e\n",
810 -log(sum_embU
/nsteps
)/beta
, -log(VembU_all
/V_all
)/beta
);
813 fprintf(fp_tpi
, "%10.3f %12.5e %12.5e %12.5e %12.5e",
815 VembU_all
== 0 ? 20/beta
: -log(VembU_all
/V_all
)/beta
,
816 sum_embU
== 0 ? 20/beta
: -log(sum_embU
/nsteps
)/beta
,
818 for (e
= 0; e
< nener
; e
++)
820 fprintf(fp_tpi
, " %12.5e", sum_UgembU
[e
]/nsteps
);
822 fprintf(fp_tpi
, "\n");
826 bNotLastFrame
= read_next_frame(oenv
, status
, &rerun_fr
);
827 } /* End of the loop */
828 walltime_accounting_end(walltime_accounting
);
834 gmx_fio_fclose(fp_tpi
);
839 fprintf(fplog
, "\n");
840 fprintf(fplog
, " <V> = %12.5e nm^3\n", V_all
/frame
);
841 fprintf(fplog
, " <mu> = %12.5e kJ/mol\n", -log(VembU_all
/V_all
)/beta
);
844 /* Write the Boltzmann factor histogram */
847 /* When running in parallel sum the bins over the processes */
850 realloc_bins(&bin
, &nbin
, i
);
851 gmx_sumd(nbin
, bin
, cr
);
855 fp_tpi
= xvgropen(opt2fn("-tpid", nfile
, fnm
),
856 "TPI energy distribution",
857 "\\betaU - log(V/<V>)", "count", oenv
);
858 sprintf(str
, "number \\betaU > %g: %9.3e", bU_bin_limit
, bin
[0]);
859 xvgr_subtitle(fp_tpi
, str
, oenv
);
860 xvgr_legend(fp_tpi
, 2, (const char **)tpid_leg
, oenv
);
861 for (i
= nbin
-1; i
> 0; i
--)
863 bUlogV
= -i
/invbinw
+ bU_logV_bin_limit
- refvolshift
+ log(V_all
/frame
);
864 fprintf(fp_tpi
, "%6.2f %10d %12.5e\n",
867 bin
[i
]*exp(-bUlogV
)*V_all
/VembU_all
);
869 gmx_fio_fclose(fp_tpi
);
875 walltime_accounting_set_nsteps_done(walltime_accounting
, frame
*inputrec
->nsteps
);