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 * GROningen Mixture of Alchemy and Childrens' Stories
36 /* This file is completely threadsafe - keep it that way! */
48 #include "gromacs/legacyheaders/thread_mpi/threads.h"
50 /* The source code in this file should be thread-safe.
51 Please keep it that way. */
55 static gmx_bool bOverAllocDD = FALSE;
56 static tMPI_Thread_mutex_t over_alloc_mutex = TMPI_THREAD_MUTEX_INITIALIZER;
59 void set_over_alloc_dd(gmx_bool set)
61 tMPI_Thread_mutex_lock(&over_alloc_mutex);
62 /* we just make sure that we don't set this at the same time.
63 We don't worry too much about reading this rarely-set variable */
65 tMPI_Thread_mutex_unlock(&over_alloc_mutex);
68 int over_alloc_dd(int n)
72 return OVER_ALLOC_FAC*n + 100;
80 int gmx_large_int_to_int(gmx_large_int_t step, const char *warn)
86 if (warn != NULL && (step < INT_MIN || step > INT_MAX))
88 fprintf(stderr, "\nWARNING during %s:\n", warn);
89 fprintf(stderr, "step value ");
90 fprintf(stderr, gmx_large_int_pfmt, step);
91 fprintf(stderr, " does not fit in int, converted to %d\n\n", i);
97 char *gmx_step_str(gmx_large_int_t i, char *buf)
99 sprintf(buf, gmx_large_int_pfmt, i);
104 void init_block(t_block *block)
109 block->nalloc_index = 1;
110 snew(block->index, block->nalloc_index);
114 void init_blocka(t_blocka *block)
120 block->nalloc_index = 1;
121 snew(block->index, block->nalloc_index);
127 void init_atom(t_atoms *at)
137 at->atomtypeB = NULL;
141 void init_atomtypes(t_atomtypes *at)
146 at->atomnumber = NULL;
147 at->gb_radius = NULL;
151 void init_groups(gmx_groups_t *groups)
155 groups->ngrpname = 0;
156 groups->grpname = NULL;
157 for (g = 0; (g < egcNR); g++)
159 groups->grps[g].nm_ind = NULL;
160 groups->ngrpnr[g] = 0;
161 groups->grpnr[g] = NULL;
166 void init_mtop(gmx_mtop_t *mtop)
170 mtop->moltype = NULL;
172 mtop->molblock = NULL;
173 mtop->maxres_renum = 0;
175 init_groups(&mtop->groups);
176 init_block(&mtop->mols);
177 open_symtab(&mtop->symtab);
180 void init_top (t_topology *top)
185 init_atom (&(top->atoms));
186 init_atomtypes(&(top->atomtypes));
187 init_block(&top->cgs);
188 init_block(&top->mols);
189 init_blocka(&top->excls);
190 open_symtab(&top->symtab);
193 void init_inputrec(t_inputrec *ir)
195 memset(ir, 0, (size_t)sizeof(*ir));
196 snew(ir->fepvals, 1);
197 snew(ir->expandedvals, 1);
198 snew(ir->simtempvals, 1);
201 void stupid_fill_block(t_block *grp, int natom, gmx_bool bOneIndexGroup)
207 grp->nalloc_index = 2;
208 snew(grp->index, grp->nalloc_index);
210 grp->index[1] = natom;
215 grp->nalloc_index = natom+1;
216 snew(grp->index, grp->nalloc_index);
217 snew(grp->index, natom+1);
218 for (i = 0; (i <= natom); i++)
226 void stupid_fill_blocka(t_blocka *grp, int natom)
230 grp->nalloc_a = natom;
231 snew(grp->a, grp->nalloc_a);
232 for (i = 0; (i < natom); i++)
238 grp->nalloc_index = natom + 1;
239 snew(grp->index, grp->nalloc_index);
240 for (i = 0; (i <= natom); i++)
247 void copy_blocka(const t_blocka *src, t_blocka *dest)
252 dest->nalloc_index = dest->nr + 1;
253 snew(dest->index, dest->nalloc_index);
254 for (i = 0; i < dest->nr+1; i++)
256 dest->index[i] = src->index[i];
258 dest->nra = src->nra;
259 dest->nalloc_a = dest->nra + 1;
260 snew(dest->a, dest->nalloc_a);
261 for (i = 0; i < dest->nra+1; i++)
263 dest->a[i] = src->a[i];
267 void done_block(t_block *block)
271 block->nalloc_index = 0;
274 void done_blocka(t_blocka *block)
282 block->nalloc_index = 0;
286 void done_atom (t_atoms *at)
294 sfree(at->atomtypeB);
301 void done_atomtypes(t_atomtypes *atype)
304 sfree(atype->radius);
306 sfree(atype->surftens);
307 sfree(atype->atomnumber);
308 sfree(atype->gb_radius);
312 void done_moltype(gmx_moltype_t *molt)
316 done_atom(&molt->atoms);
317 done_block(&molt->cgs);
318 done_blocka(&molt->excls);
320 for (f = 0; f < F_NRE; f++)
322 sfree(molt->ilist[f].iatoms);
323 molt->ilist[f].nalloc = 0;
327 void done_molblock(gmx_molblock_t *molb)
329 if (molb->nposres_xA > 0)
331 molb->nposres_xA = 0;
332 free(molb->posres_xA);
334 if (molb->nposres_xB > 0)
336 molb->nposres_xB = 0;
337 free(molb->posres_xB);
341 void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
347 done_symtab(&mtop->symtab);
350 sfree(mtop->ffparams.functype);
351 sfree(mtop->ffparams.iparams);
353 for (i = 0; i < mtop->nmoltype; i++)
355 done_moltype(&mtop->moltype[i]);
357 sfree(mtop->moltype);
358 for (i = 0; i < mtop->nmolblock; i++)
360 done_molblock(&mtop->molblock[i]);
362 sfree(mtop->molblock);
363 done_block(&mtop->mols);
366 void done_top(t_topology *top)
370 sfree(top->idef.functype);
371 sfree(top->idef.iparams);
372 for (f = 0; f < F_NRE; ++f)
374 sfree(top->idef.il[f].iatoms);
375 top->idef.il[f].iatoms = NULL;
376 top->idef.il[f].nalloc = 0;
379 done_atom (&(top->atoms));
382 done_atomtypes(&(top->atomtypes));
384 done_symtab(&(top->symtab));
385 done_block(&(top->cgs));
386 done_block(&(top->mols));
387 done_blocka(&(top->excls));
390 static void done_pullgrp(t_pullgrp *pgrp)
393 sfree(pgrp->ind_loc);
395 sfree(pgrp->weight_loc);
398 static void done_pull(t_pull *pull)
402 for (i = 0; i < pull->ngrp+1; i++)
404 done_pullgrp(pull->grp);
405 done_pullgrp(pull->dyna);
409 void done_inputrec(t_inputrec *ir)
413 for (m = 0; (m < DIM); m++)
421 sfree(ir->ex[m].phi);
429 sfree(ir->et[m].phi);
433 sfree(ir->opts.nrdf);
434 sfree(ir->opts.ref_t);
435 sfree(ir->opts.annealing);
436 sfree(ir->opts.anneal_npoints);
437 sfree(ir->opts.anneal_time);
438 sfree(ir->opts.anneal_temp);
439 sfree(ir->opts.tau_t);
441 sfree(ir->opts.nFreeze);
442 sfree(ir->opts.QMmethod);
443 sfree(ir->opts.QMbasis);
444 sfree(ir->opts.QMcharge);
445 sfree(ir->opts.QMmult);
447 sfree(ir->opts.CASorbitals);
448 sfree(ir->opts.CASelectrons);
449 sfree(ir->opts.SAon);
450 sfree(ir->opts.SAoff);
451 sfree(ir->opts.SAsteps);
452 sfree(ir->opts.bOPT);
462 static void zero_history(history_t *hist)
464 hist->disre_initf = 0;
465 hist->ndisrepairs = 0;
466 hist->disre_rm3tav = NULL;
467 hist->orire_initf = 0;
468 hist->norire_Dtav = 0;
469 hist->orire_Dtav = NULL;
472 static void zero_ekinstate(ekinstate_t *eks)
477 eks->ekinh_old = NULL;
478 eks->ekinscalef_nhc = NULL;
479 eks->ekinscaleh_nhc = NULL;
480 eks->vscale_nhc = NULL;
485 void init_energyhistory(energyhistory_t * enerhist)
489 enerhist->ener_ave = NULL;
490 enerhist->ener_sum = NULL;
491 enerhist->ener_sum_sim = NULL;
492 enerhist->dht = NULL;
494 enerhist->nsteps = 0;
496 enerhist->nsteps_sim = 0;
497 enerhist->nsum_sim = 0;
499 enerhist->dht = NULL;
502 static void done_delta_h_history(delta_h_history_t *dht)
506 for (i = 0; i < dht->nndh; i++)
514 void done_energyhistory(energyhistory_t * enerhist)
516 sfree(enerhist->ener_ave);
517 sfree(enerhist->ener_sum);
518 sfree(enerhist->ener_sum_sim);
520 if (enerhist->dht != NULL)
522 done_delta_h_history(enerhist->dht);
523 sfree(enerhist->dht);
527 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
532 state->nnhpres = nnhpres;
533 state->nhchainlength = nhchainlength;
536 snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
537 snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
538 snew(state->therm_integral, state->ngtc);
539 for (i = 0; i < state->ngtc; i++)
541 for (j = 0; j < state->nhchainlength; j++)
543 state->nosehoover_xi[i*state->nhchainlength + j] = 0.0;
544 state->nosehoover_vxi[i*state->nhchainlength + j] = 0.0;
547 for (i = 0; i < state->ngtc; i++)
549 state->therm_integral[i] = 0.0;
554 state->nosehoover_xi = NULL;
555 state->nosehoover_vxi = NULL;
556 state->therm_integral = NULL;
559 if (state->nnhpres > 0)
561 snew(state->nhpres_xi, state->nhchainlength*nnhpres);
562 snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
563 for (i = 0; i < nnhpres; i++)
565 for (j = 0; j < state->nhchainlength; j++)
567 state->nhpres_xi[i*nhchainlength + j] = 0.0;
568 state->nhpres_vxi[i*nhchainlength + j] = 0.0;
574 state->nhpres_xi = NULL;
575 state->nhpres_vxi = NULL;
580 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
584 state->natoms = natoms;
588 snew(state->lambda, efptNR);
589 for (i = 0; i < efptNR; i++)
591 state->lambda[i] = 0;
594 clear_mat(state->box);
595 clear_mat(state->box_rel);
596 clear_mat(state->boxv);
597 clear_mat(state->pres_prev);
598 clear_mat(state->svir_prev);
599 clear_mat(state->fvir_prev);
600 init_gtc_state(state, ngtc, nnhpres, nhchainlength);
601 state->nalloc = state->natoms;
602 if (state->nalloc > 0)
604 snew(state->x, state->nalloc);
605 snew(state->v, state->nalloc);
614 zero_history(&state->hist);
615 zero_ekinstate(&state->ekinstate);
616 init_energyhistory(&state->enerhist);
617 init_df_history(&state->dfhist,nlambda);
618 state->ddp_count = 0;
619 state->ddp_count_cg_gl = 0;
621 state->cg_gl_nalloc = 0;
624 void done_state(t_state *state)
647 state->cg_gl_nalloc = 0;
650 sfree(state->lambda);
654 sfree(state->nosehoover_xi);
655 sfree(state->nosehoover_vxi);
656 sfree(state->therm_integral);
660 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
664 for (d = YY; d <= ZZ; d++)
666 for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
668 /* We need to check if this box component is deformed
669 * or if deformation of another component might cause
670 * changes in this component due to box corrections.
672 if (ir->deform[d][d2] == 0 &&
673 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
674 (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
678 box_rel[d][d2] = b[d][d2]/b[XX][XX];
682 b[d][d2] = b[XX][XX]*box_rel[d][d2];
689 void set_box_rel(t_inputrec *ir, t_state *state)
691 /* Make sure the box obeys the restrictions before we fix the ratios */
692 correct_box(NULL, 0, state->box, NULL);
694 clear_mat(state->box_rel);
696 if (PRESERVE_SHAPE(*ir))
698 do_box_rel(ir, state->box_rel, state->box, TRUE);
702 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
704 if (PRESERVE_SHAPE(*ir))
706 do_box_rel(ir, box_rel, b, FALSE);
710 void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
716 srenew(atoms->atomname, atoms->nr+natom_extra);
717 srenew(atoms->atom, atoms->nr+natom_extra);
718 if (NULL != atoms->pdbinfo)
720 srenew(atoms->pdbinfo, atoms->nr+natom_extra);
722 if (NULL != atoms->atomtype)
724 srenew(atoms->atomtype, atoms->nr+natom_extra);
726 if (NULL != atoms->atomtypeB)
728 srenew(atoms->atomtypeB, atoms->nr+natom_extra);
730 for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
732 atoms->atomname[i] = NULL;
733 memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
734 if (NULL != atoms->pdbinfo)
736 memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
738 if (NULL != atoms->atomtype)
740 atoms->atomtype[i] = NULL;
742 if (NULL != atoms->atomtypeB)
744 atoms->atomtypeB[i] = NULL;
747 atoms->nr += natom_extra;
751 srenew(atoms->resinfo, atoms->nres+nres_extra);
752 for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
754 memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
756 atoms->nres += nres_extra;
760 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
764 snew(atoms->atomname, natoms);
765 atoms->atomtype = NULL;
766 atoms->atomtypeB = NULL;
767 snew(atoms->resinfo, natoms);
768 snew(atoms->atom, natoms);
771 snew(atoms->pdbinfo, natoms);
775 atoms->pdbinfo = NULL;
779 t_atoms *copy_t_atoms(t_atoms *src)
785 init_t_atoms(dst, src->nr, (NULL != src->pdbinfo));
787 if (NULL != src->atomname)
789 snew(dst->atomname, src->nr);
791 if (NULL != src->atomtype)
793 snew(dst->atomtype, src->nr);
795 if (NULL != src->atomtypeB)
797 snew(dst->atomtypeB, src->nr);
799 for (i = 0; (i < src->nr); i++)
801 dst->atom[i] = src->atom[i];
802 if (NULL != src->pdbinfo)
804 dst->pdbinfo[i] = src->pdbinfo[i];
806 if (NULL != src->atomname)
808 dst->atomname[i] = src->atomname[i];
810 if (NULL != src->atomtype)
812 dst->atomtype[i] = src->atomtype[i];
814 if (NULL != src->atomtypeB)
816 dst->atomtypeB[i] = src->atomtypeB[i];
819 dst->nres = src->nres;
820 for (i = 0; (i < src->nres); i++)
822 dst->resinfo[i] = src->resinfo[i];
827 void t_atoms_set_resinfo(t_atoms *atoms, int atom_ind, t_symtab *symtab,
828 const char *resname, int resnr, unsigned char ic,
829 int chainnum, char chainid)
833 ri = &atoms->resinfo[atoms->atom[atom_ind].resind];
834 ri->name = put_symtab(symtab, resname);
838 ri->chainnum = chainnum;
839 ri->chainid = chainid;
842 void free_t_atoms(t_atoms *atoms, gmx_bool bFreeNames)
846 if (bFreeNames && atoms->atomname != NULL)
848 for (i = 0; i < atoms->nr; i++)
850 if (atoms->atomname[i] != NULL)
852 sfree(*atoms->atomname[i]);
853 *atoms->atomname[i] = NULL;
857 if (bFreeNames && atoms->resinfo != NULL)
859 for (i = 0; i < atoms->nres; i++)
861 if (atoms->resinfo[i].name != NULL)
863 sfree(*atoms->resinfo[i].name);
864 *atoms->resinfo[i].name = NULL;
868 if (bFreeNames && atoms->atomtype != NULL)
870 for (i = 0; i < atoms->nr; i++)
872 if (atoms->atomtype[i] != NULL)
874 sfree(*atoms->atomtype[i]);
875 *atoms->atomtype[i] = NULL;
879 if (bFreeNames && atoms->atomtypeB != NULL)
881 for (i = 0; i < atoms->nr; i++)
883 if (atoms->atomtypeB[i] != NULL)
885 sfree(*atoms->atomtypeB[i]);
886 *atoms->atomtypeB[i] = NULL;
890 sfree(atoms->atomname);
891 sfree(atoms->atomtype);
892 sfree(atoms->atomtypeB);
893 sfree(atoms->resinfo);
895 sfree(atoms->pdbinfo);
898 atoms->atomname = NULL;
899 atoms->atomtype = NULL;
900 atoms->atomtypeB = NULL;
901 atoms->resinfo = NULL;
903 atoms->pdbinfo = NULL;
906 real max_cutoff(real cutoff1, real cutoff2)
908 if (cutoff1 == 0 || cutoff2 == 0)
914 return max(cutoff1, cutoff2);
918 void init_df_history(df_history_t *dfhist, int nlambda)
922 dfhist->nlambda = nlambda;
924 dfhist->wl_delta = 0;
928 snew(dfhist->sum_weights, dfhist->nlambda);
929 snew(dfhist->sum_dg, dfhist->nlambda);
930 snew(dfhist->sum_minvar, dfhist->nlambda);
931 snew(dfhist->sum_variance, dfhist->nlambda);
932 snew(dfhist->n_at_lam, dfhist->nlambda);
933 snew(dfhist->wl_histo, dfhist->nlambda);
935 /* allocate transition matrices here */
936 snew(dfhist->Tij, dfhist->nlambda);
937 snew(dfhist->Tij_empirical, dfhist->nlambda);
939 /* allocate accumulators for various transition matrix
940 free energy methods here */
941 snew(dfhist->accum_p, dfhist->nlambda);
942 snew(dfhist->accum_m, dfhist->nlambda);
943 snew(dfhist->accum_p2, dfhist->nlambda);
944 snew(dfhist->accum_m2, dfhist->nlambda);
946 for (i = 0; i < dfhist->nlambda; i++)
948 snew(dfhist->Tij[i], dfhist->nlambda);
949 snew(dfhist->Tij_empirical[i], dfhist->nlambda);
950 snew((dfhist->accum_p)[i], dfhist->nlambda);
951 snew((dfhist->accum_m)[i], dfhist->nlambda);
952 snew((dfhist->accum_p2)[i], dfhist->nlambda);
953 snew((dfhist->accum_m2)[i], dfhist->nlambda);
958 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
962 /* Currently, there should not be any difference in nlambda between the two,
963 but this is included for completeness for potential later functionality */
964 df_dest->nlambda = df_source->nlambda;
965 df_dest->bEquil = df_source->bEquil;
966 df_dest->wl_delta = df_source->wl_delta;
968 for (i = 0; i < df_dest->nlambda; i++)
970 df_dest->sum_weights[i] = df_source->sum_weights[i];
971 df_dest->sum_dg[i] = df_source->sum_dg[i];
972 df_dest->sum_minvar[i] = df_source->sum_minvar[i];
973 df_dest->sum_variance[i] = df_source->sum_variance[i];
974 df_dest->n_at_lam[i] = df_source->n_at_lam[i];
975 df_dest->wl_histo[i] = df_source->wl_histo[i];
978 for (i = 0; i < df_dest->nlambda; i++)
980 for (j = 0; j < df_dest->nlambda; j++)
982 df_dest->accum_p[i][j] = df_source->accum_p[i][j];
983 df_dest->accum_m[i][j] = df_source->accum_m[i][j];
984 df_dest->accum_p2[i][j] = df_source->accum_p2[i][j];
985 df_dest->accum_m2[i][j] = df_source->accum_m2[i][j];
986 df_dest->Tij[i][j] = df_source->Tij[i][j];
987 df_dest->Tij_empirical[i][j] = df_source->Tij_empirical[i][j];
992 void done_df_history(df_history_t *dfhist)
996 if (dfhist->nlambda > 0)
998 sfree(dfhist->n_at_lam);
999 sfree(dfhist->wl_histo);
1000 sfree(dfhist->sum_weights);
1001 sfree(dfhist->sum_dg);
1002 sfree(dfhist->sum_minvar);
1003 sfree(dfhist->sum_variance);
1005 for (i=0;i<dfhist->nlambda;i++)
1007 sfree(dfhist->Tij[i]);
1008 sfree(dfhist->Tij_empirical[i]);
1009 sfree(dfhist->accum_p[i]);
1010 sfree(dfhist->accum_m[i]);
1011 sfree(dfhist->accum_p2[i]);
1012 sfree(dfhist->accum_m2[i]);
1016 dfhist->nlambda = 0;
1017 dfhist->wl_delta = 0;