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.
37 /* This file is completely threadsafe - keep it that way! */
46 #include "thread_mpi/threads.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/pbcutil/pbc.h"
51 #include "gromacs/random/random.h"
52 #include "gromacs/topology/block.h"
53 #include "gromacs/topology/symtab.h"
54 #include "gromacs/utility/smalloc.h"
56 /* The source code in this file should be thread-safe.
57 Please keep it that way. */
61 static gmx_bool bOverAllocDD = FALSE;
62 static tMPI_Thread_mutex_t over_alloc_mutex = TMPI_THREAD_MUTEX_INITIALIZER;
65 void set_over_alloc_dd(gmx_bool set)
67 tMPI_Thread_mutex_lock(&over_alloc_mutex);
68 /* we just make sure that we don't set this at the same time.
69 We don't worry too much about reading this rarely-set variable */
71 tMPI_Thread_mutex_unlock(&over_alloc_mutex);
74 int over_alloc_dd(int n)
78 return OVER_ALLOC_FAC*n + 100;
86 int gmx_int64_to_int(gmx_int64_t step, const char *warn)
92 if (warn != NULL && (step < INT_MIN || step > INT_MAX))
94 fprintf(stderr, "\nWARNING during %s:\n", warn);
95 fprintf(stderr, "step value ");
96 fprintf(stderr, "%"GMX_PRId64, step);
97 fprintf(stderr, " does not fit in int, converted to %d\n\n", i);
103 char *gmx_step_str(gmx_int64_t i, char *buf)
105 sprintf(buf, "%"GMX_PRId64, i);
110 void init_atom(t_atoms *at)
120 at->atomtypeB = NULL;
124 void init_atomtypes(t_atomtypes *at)
129 at->atomnumber = NULL;
130 at->gb_radius = NULL;
134 void init_groups(gmx_groups_t *groups)
138 groups->ngrpname = 0;
139 groups->grpname = NULL;
140 for (g = 0; (g < egcNR); g++)
142 groups->grps[g].nm_ind = NULL;
143 groups->ngrpnr[g] = 0;
144 groups->grpnr[g] = NULL;
149 void init_mtop(gmx_mtop_t *mtop)
153 mtop->moltype = NULL;
155 mtop->molblock = NULL;
156 mtop->maxres_renum = 0;
158 init_groups(&mtop->groups);
159 init_block(&mtop->mols);
160 open_symtab(&mtop->symtab);
163 void init_top(t_topology *top)
168 init_atom (&(top->atoms));
169 init_atomtypes(&(top->atomtypes));
170 init_block(&top->cgs);
171 init_block(&top->mols);
172 init_blocka(&top->excls);
173 open_symtab(&top->symtab);
176 void init_inputrec(t_inputrec *ir)
178 memset(ir, 0, (size_t)sizeof(*ir));
179 snew(ir->fepvals, 1);
180 snew(ir->expandedvals, 1);
181 snew(ir->simtempvals, 1);
184 void done_atom (t_atoms *at)
192 sfree(at->atomtypeB);
199 void done_atomtypes(t_atomtypes *atype)
202 sfree(atype->radius);
204 sfree(atype->surftens);
205 sfree(atype->atomnumber);
206 sfree(atype->gb_radius);
210 void done_moltype(gmx_moltype_t *molt)
214 done_atom(&molt->atoms);
215 done_block(&molt->cgs);
216 done_blocka(&molt->excls);
218 for (f = 0; f < F_NRE; f++)
220 sfree(molt->ilist[f].iatoms);
221 molt->ilist[f].nalloc = 0;
225 void done_molblock(gmx_molblock_t *molb)
227 if (molb->nposres_xA > 0)
229 molb->nposres_xA = 0;
230 sfree(molb->posres_xA);
232 if (molb->nposres_xB > 0)
234 molb->nposres_xB = 0;
235 sfree(molb->posres_xB);
239 void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
245 done_symtab(&mtop->symtab);
248 sfree(mtop->ffparams.functype);
249 sfree(mtop->ffparams.iparams);
251 for (i = 0; i < mtop->nmoltype; i++)
253 done_moltype(&mtop->moltype[i]);
255 sfree(mtop->moltype);
256 for (i = 0; i < mtop->nmolblock; i++)
258 done_molblock(&mtop->molblock[i]);
260 sfree(mtop->molblock);
261 done_block(&mtop->mols);
264 void done_top(t_topology *top)
268 sfree(top->idef.functype);
269 sfree(top->idef.iparams);
270 for (f = 0; f < F_NRE; ++f)
272 sfree(top->idef.il[f].iatoms);
273 top->idef.il[f].iatoms = NULL;
274 top->idef.il[f].nalloc = 0;
277 done_atom (&(top->atoms));
280 done_atomtypes(&(top->atomtypes));
282 done_symtab(&(top->symtab));
283 done_block(&(top->cgs));
284 done_block(&(top->mols));
285 done_blocka(&(top->excls));
288 static void done_pull_group(t_pull_group *pgrp)
293 sfree(pgrp->ind_loc);
295 sfree(pgrp->weight_loc);
299 static void done_pull(t_pull *pull)
303 for (i = 0; i < pull->ngroup+1; i++)
305 done_pull_group(pull->group);
306 done_pull_group(pull->dyna);
310 void done_inputrec(t_inputrec *ir)
314 for (m = 0; (m < DIM); m++)
322 sfree(ir->ex[m].phi);
330 sfree(ir->et[m].phi);
334 sfree(ir->opts.nrdf);
335 sfree(ir->opts.ref_t);
336 sfree(ir->opts.annealing);
337 sfree(ir->opts.anneal_npoints);
338 sfree(ir->opts.anneal_time);
339 sfree(ir->opts.anneal_temp);
340 sfree(ir->opts.tau_t);
342 sfree(ir->opts.nFreeze);
343 sfree(ir->opts.QMmethod);
344 sfree(ir->opts.QMbasis);
345 sfree(ir->opts.QMcharge);
346 sfree(ir->opts.QMmult);
348 sfree(ir->opts.CASorbitals);
349 sfree(ir->opts.CASelectrons);
350 sfree(ir->opts.SAon);
351 sfree(ir->opts.SAoff);
352 sfree(ir->opts.SAsteps);
353 sfree(ir->opts.bOPT);
363 static void zero_history(history_t *hist)
365 hist->disre_initf = 0;
366 hist->ndisrepairs = 0;
367 hist->disre_rm3tav = NULL;
368 hist->orire_initf = 0;
369 hist->norire_Dtav = 0;
370 hist->orire_Dtav = NULL;
373 static void zero_ekinstate(ekinstate_t *eks)
378 eks->ekinh_old = NULL;
379 eks->ekinscalef_nhc = NULL;
380 eks->ekinscaleh_nhc = NULL;
381 eks->vscale_nhc = NULL;
386 static void init_swapstate(swapstate_t *swapstate)
390 swapstate->eSwapCoords = 0;
391 swapstate->nAverage = 0;
393 /* Ion/water position swapping */
394 for (ic = 0; ic < eCompNR; ic++)
396 for (ii = 0; ii < eIonNR; ii++)
398 swapstate->nat_req[ic][ii] = 0;
399 swapstate->nat_req_p[ic][ii] = NULL;
400 swapstate->inflow_netto[ic][ii] = 0;
401 swapstate->inflow_netto_p[ic][ii] = NULL;
402 swapstate->nat_past[ic][ii] = NULL;
403 swapstate->nat_past_p[ic][ii] = NULL;
404 swapstate->fluxfromAtoB[ic][ii] = 0;
405 swapstate->fluxfromAtoB_p[ic][ii] = NULL;
408 swapstate->fluxleak = NULL;
409 swapstate->nions = 0;
410 swapstate->comp_from = NULL;
411 swapstate->channel_label = NULL;
412 swapstate->bFromCpt = 0;
413 swapstate->nat[eChan0] = 0;
414 swapstate->nat[eChan1] = 0;
415 swapstate->xc_old_whole[eChan0] = NULL;
416 swapstate->xc_old_whole[eChan1] = NULL;
417 swapstate->xc_old_whole_p[eChan0] = NULL;
418 swapstate->xc_old_whole_p[eChan1] = NULL;
421 void init_energyhistory(energyhistory_t * enerhist)
425 enerhist->ener_ave = NULL;
426 enerhist->ener_sum = NULL;
427 enerhist->ener_sum_sim = NULL;
428 enerhist->dht = NULL;
430 enerhist->nsteps = 0;
432 enerhist->nsteps_sim = 0;
433 enerhist->nsum_sim = 0;
435 enerhist->dht = NULL;
438 static void done_delta_h_history(delta_h_history_t *dht)
442 for (i = 0; i < dht->nndh; i++)
450 void done_energyhistory(energyhistory_t * enerhist)
452 sfree(enerhist->ener_ave);
453 sfree(enerhist->ener_sum);
454 sfree(enerhist->ener_sum_sim);
456 if (enerhist->dht != NULL)
458 done_delta_h_history(enerhist->dht);
459 sfree(enerhist->dht);
463 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
468 state->nnhpres = nnhpres;
469 state->nhchainlength = nhchainlength;
472 snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
473 snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
474 snew(state->therm_integral, state->ngtc);
475 for (i = 0; i < state->ngtc; i++)
477 for (j = 0; j < state->nhchainlength; j++)
479 state->nosehoover_xi[i*state->nhchainlength + j] = 0.0;
480 state->nosehoover_vxi[i*state->nhchainlength + j] = 0.0;
483 for (i = 0; i < state->ngtc; i++)
485 state->therm_integral[i] = 0.0;
490 state->nosehoover_xi = NULL;
491 state->nosehoover_vxi = NULL;
492 state->therm_integral = NULL;
495 if (state->nnhpres > 0)
497 snew(state->nhpres_xi, state->nhchainlength*nnhpres);
498 snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
499 for (i = 0; i < nnhpres; i++)
501 for (j = 0; j < state->nhchainlength; j++)
503 state->nhpres_xi[i*nhchainlength + j] = 0.0;
504 state->nhpres_vxi[i*nhchainlength + j] = 0.0;
510 state->nhpres_xi = NULL;
511 state->nhpres_vxi = NULL;
516 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
520 state->natoms = natoms;
523 snew(state->lambda, efptNR);
524 for (i = 0; i < efptNR; i++)
526 state->lambda[i] = 0;
529 clear_mat(state->box);
530 clear_mat(state->box_rel);
531 clear_mat(state->boxv);
532 clear_mat(state->pres_prev);
533 clear_mat(state->svir_prev);
534 clear_mat(state->fvir_prev);
535 init_gtc_state(state, ngtc, nnhpres, nhchainlength);
536 state->nalloc = state->natoms;
537 if (state->nalloc > 0)
539 snew(state->x, state->nalloc);
540 snew(state->v, state->nalloc);
549 zero_history(&state->hist);
550 zero_ekinstate(&state->ekinstate);
551 init_energyhistory(&state->enerhist);
552 init_df_history(&state->dfhist, nlambda);
553 init_swapstate(&state->swapstate);
554 state->ddp_count = 0;
555 state->ddp_count_cg_gl = 0;
557 state->cg_gl_nalloc = 0;
560 void done_state(t_state *state)
583 state->cg_gl_nalloc = 0;
586 sfree(state->lambda);
590 sfree(state->nosehoover_xi);
591 sfree(state->nosehoover_vxi);
592 sfree(state->therm_integral);
596 t_state *serial_init_local_state(t_state *state_global)
599 t_state *state_local;
601 snew(state_local, 1);
603 /* Copy all the contents */
604 *state_local = *state_global;
605 snew(state_local->lambda, efptNR);
606 /* local storage for lambda */
607 for (i = 0; i < efptNR; i++)
609 state_local->lambda[i] = state_global->lambda[i];
615 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
619 for (d = YY; d <= ZZ; d++)
621 for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
623 /* We need to check if this box component is deformed
624 * or if deformation of another component might cause
625 * changes in this component due to box corrections.
627 if (ir->deform[d][d2] == 0 &&
628 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
629 (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
633 box_rel[d][d2] = b[d][d2]/b[XX][XX];
637 b[d][d2] = b[XX][XX]*box_rel[d][d2];
644 void set_box_rel(t_inputrec *ir, t_state *state)
646 /* Make sure the box obeys the restrictions before we fix the ratios */
647 correct_box(NULL, 0, state->box, NULL);
649 clear_mat(state->box_rel);
651 if (PRESERVE_SHAPE(*ir))
653 do_box_rel(ir, state->box_rel, state->box, TRUE);
657 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
659 if (PRESERVE_SHAPE(*ir))
661 do_box_rel(ir, box_rel, b, FALSE);
665 void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
671 srenew(atoms->atomname, atoms->nr+natom_extra);
672 srenew(atoms->atom, atoms->nr+natom_extra);
673 if (NULL != atoms->pdbinfo)
675 srenew(atoms->pdbinfo, atoms->nr+natom_extra);
677 if (NULL != atoms->atomtype)
679 srenew(atoms->atomtype, atoms->nr+natom_extra);
681 if (NULL != atoms->atomtypeB)
683 srenew(atoms->atomtypeB, atoms->nr+natom_extra);
685 for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
687 atoms->atomname[i] = NULL;
688 memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
689 if (NULL != atoms->pdbinfo)
691 memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
693 if (NULL != atoms->atomtype)
695 atoms->atomtype[i] = NULL;
697 if (NULL != atoms->atomtypeB)
699 atoms->atomtypeB[i] = NULL;
702 atoms->nr += natom_extra;
706 srenew(atoms->resinfo, atoms->nres+nres_extra);
707 for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
709 memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
711 atoms->nres += nres_extra;
715 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
719 snew(atoms->atomname, natoms);
720 atoms->atomtype = NULL;
721 atoms->atomtypeB = NULL;
722 snew(atoms->resinfo, natoms);
723 snew(atoms->atom, natoms);
726 snew(atoms->pdbinfo, natoms);
730 atoms->pdbinfo = NULL;
734 t_atoms *copy_t_atoms(t_atoms *src)
740 init_t_atoms(dst, src->nr, (NULL != src->pdbinfo));
742 if (NULL != src->atomname)
744 snew(dst->atomname, src->nr);
746 if (NULL != src->atomtype)
748 snew(dst->atomtype, src->nr);
750 if (NULL != src->atomtypeB)
752 snew(dst->atomtypeB, src->nr);
754 for (i = 0; (i < src->nr); i++)
756 dst->atom[i] = src->atom[i];
757 if (NULL != src->pdbinfo)
759 dst->pdbinfo[i] = src->pdbinfo[i];
761 if (NULL != src->atomname)
763 dst->atomname[i] = src->atomname[i];
765 if (NULL != src->atomtype)
767 dst->atomtype[i] = src->atomtype[i];
769 if (NULL != src->atomtypeB)
771 dst->atomtypeB[i] = src->atomtypeB[i];
774 dst->nres = src->nres;
775 for (i = 0; (i < src->nres); i++)
777 dst->resinfo[i] = src->resinfo[i];
782 void t_atoms_set_resinfo(t_atoms *atoms, int atom_ind, t_symtab *symtab,
783 const char *resname, int resnr, unsigned char ic,
784 int chainnum, char chainid)
788 ri = &atoms->resinfo[atoms->atom[atom_ind].resind];
789 ri->name = put_symtab(symtab, resname);
793 ri->chainnum = chainnum;
794 ri->chainid = chainid;
797 void free_t_atoms(t_atoms *atoms, gmx_bool bFreeNames)
801 if (bFreeNames && atoms->atomname != NULL)
803 for (i = 0; i < atoms->nr; i++)
805 if (atoms->atomname[i] != NULL)
807 sfree(*atoms->atomname[i]);
808 *atoms->atomname[i] = NULL;
812 if (bFreeNames && atoms->resinfo != NULL)
814 for (i = 0; i < atoms->nres; i++)
816 if (atoms->resinfo[i].name != NULL)
818 sfree(*atoms->resinfo[i].name);
819 *atoms->resinfo[i].name = NULL;
823 if (bFreeNames && atoms->atomtype != NULL)
825 for (i = 0; i < atoms->nr; i++)
827 if (atoms->atomtype[i] != NULL)
829 sfree(*atoms->atomtype[i]);
830 *atoms->atomtype[i] = NULL;
834 if (bFreeNames && atoms->atomtypeB != NULL)
836 for (i = 0; i < atoms->nr; i++)
838 if (atoms->atomtypeB[i] != NULL)
840 sfree(*atoms->atomtypeB[i]);
841 *atoms->atomtypeB[i] = NULL;
845 sfree(atoms->atomname);
846 sfree(atoms->atomtype);
847 sfree(atoms->atomtypeB);
848 sfree(atoms->resinfo);
850 sfree(atoms->pdbinfo);
853 atoms->atomname = NULL;
854 atoms->atomtype = NULL;
855 atoms->atomtypeB = NULL;
856 atoms->resinfo = NULL;
858 atoms->pdbinfo = NULL;
861 real max_cutoff(real cutoff1, real cutoff2)
863 if (cutoff1 == 0 || cutoff2 == 0)
869 return max(cutoff1, cutoff2);
873 void init_df_history(df_history_t *dfhist, int nlambda)
877 dfhist->nlambda = nlambda;
879 dfhist->wl_delta = 0;
883 snew(dfhist->sum_weights, dfhist->nlambda);
884 snew(dfhist->sum_dg, dfhist->nlambda);
885 snew(dfhist->sum_minvar, dfhist->nlambda);
886 snew(dfhist->sum_variance, dfhist->nlambda);
887 snew(dfhist->n_at_lam, dfhist->nlambda);
888 snew(dfhist->wl_histo, dfhist->nlambda);
890 /* allocate transition matrices here */
891 snew(dfhist->Tij, dfhist->nlambda);
892 snew(dfhist->Tij_empirical, dfhist->nlambda);
894 /* allocate accumulators for various transition matrix
895 free energy methods here */
896 snew(dfhist->accum_p, dfhist->nlambda);
897 snew(dfhist->accum_m, dfhist->nlambda);
898 snew(dfhist->accum_p2, dfhist->nlambda);
899 snew(dfhist->accum_m2, dfhist->nlambda);
901 for (i = 0; i < dfhist->nlambda; i++)
903 snew(dfhist->Tij[i], dfhist->nlambda);
904 snew(dfhist->Tij_empirical[i], dfhist->nlambda);
905 snew((dfhist->accum_p)[i], dfhist->nlambda);
906 snew((dfhist->accum_m)[i], dfhist->nlambda);
907 snew((dfhist->accum_p2)[i], dfhist->nlambda);
908 snew((dfhist->accum_m2)[i], dfhist->nlambda);
913 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
917 /* Currently, there should not be any difference in nlambda between the two,
918 but this is included for completeness for potential later functionality */
919 df_dest->nlambda = df_source->nlambda;
920 df_dest->bEquil = df_source->bEquil;
921 df_dest->wl_delta = df_source->wl_delta;
923 for (i = 0; i < df_dest->nlambda; i++)
925 df_dest->sum_weights[i] = df_source->sum_weights[i];
926 df_dest->sum_dg[i] = df_source->sum_dg[i];
927 df_dest->sum_minvar[i] = df_source->sum_minvar[i];
928 df_dest->sum_variance[i] = df_source->sum_variance[i];
929 df_dest->n_at_lam[i] = df_source->n_at_lam[i];
930 df_dest->wl_histo[i] = df_source->wl_histo[i];
933 for (i = 0; i < df_dest->nlambda; i++)
935 for (j = 0; j < df_dest->nlambda; j++)
937 df_dest->accum_p[i][j] = df_source->accum_p[i][j];
938 df_dest->accum_m[i][j] = df_source->accum_m[i][j];
939 df_dest->accum_p2[i][j] = df_source->accum_p2[i][j];
940 df_dest->accum_m2[i][j] = df_source->accum_m2[i][j];
941 df_dest->Tij[i][j] = df_source->Tij[i][j];
942 df_dest->Tij_empirical[i][j] = df_source->Tij_empirical[i][j];
947 void done_df_history(df_history_t *dfhist)
951 if (dfhist->nlambda > 0)
953 sfree(dfhist->n_at_lam);
954 sfree(dfhist->wl_histo);
955 sfree(dfhist->sum_weights);
956 sfree(dfhist->sum_dg);
957 sfree(dfhist->sum_minvar);
958 sfree(dfhist->sum_variance);
960 for (i = 0; i < dfhist->nlambda; i++)
962 sfree(dfhist->Tij[i]);
963 sfree(dfhist->Tij_empirical[i]);
964 sfree(dfhist->accum_p[i]);
965 sfree(dfhist->accum_m[i]);
966 sfree(dfhist->accum_p2[i]);
967 sfree(dfhist->accum_m2[i]);
972 dfhist->wl_delta = 0;