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! */
49 #include "thread_mpi.h"
52 /* The source code in this file should be thread-safe.
53 Please keep it that way. */
57 static gmx_bool bOverAllocDD = FALSE;
59 static tMPI_Thread_mutex_t over_alloc_mutex = TMPI_THREAD_MUTEX_INITIALIZER;
63 void set_over_alloc_dd(gmx_bool set)
66 tMPI_Thread_mutex_lock(&over_alloc_mutex);
67 /* we just make sure that we don't set this at the same time.
68 We don't worry too much about reading this rarely-set variable */
72 tMPI_Thread_mutex_unlock(&over_alloc_mutex);
76 int over_alloc_dd(int n)
80 return OVER_ALLOC_FAC*n + 100;
88 int gmx_large_int_to_int(gmx_large_int_t step, const char *warn)
94 if (warn != NULL && (step < INT_MIN || step > INT_MAX))
96 fprintf(stderr, "\nWARNING during %s:\n", warn);
97 fprintf(stderr, "step value ");
98 fprintf(stderr, gmx_large_int_pfmt, step);
99 fprintf(stderr, " does not fit in int, converted to %d\n\n", i);
105 char *gmx_step_str(gmx_large_int_t i, char *buf)
107 sprintf(buf, gmx_large_int_pfmt, i);
112 void init_block(t_block *block)
117 block->nalloc_index = 1;
118 snew(block->index, block->nalloc_index);
122 void init_blocka(t_blocka *block)
128 block->nalloc_index = 1;
129 snew(block->index, block->nalloc_index);
135 void init_atom(t_atoms *at)
145 at->atomtypeB = NULL;
149 void init_atomtypes(t_atomtypes *at)
154 at->atomnumber = NULL;
155 at->gb_radius = NULL;
159 void init_groups(gmx_groups_t *groups)
163 groups->ngrpname = 0;
164 groups->grpname = NULL;
165 for (g = 0; (g < egcNR); g++)
167 groups->grps[g].nm_ind = NULL;
168 groups->ngrpnr[g] = 0;
169 groups->grpnr[g] = NULL;
174 void init_mtop(gmx_mtop_t *mtop)
178 mtop->moltype = NULL;
180 mtop->molblock = NULL;
181 mtop->maxres_renum = 0;
183 init_groups(&mtop->groups);
184 init_block(&mtop->mols);
185 open_symtab(&mtop->symtab);
188 void init_top (t_topology *top)
193 init_atom (&(top->atoms));
194 init_atomtypes(&(top->atomtypes));
195 init_block(&top->cgs);
196 init_block(&top->mols);
197 init_blocka(&top->excls);
198 open_symtab(&top->symtab);
201 void init_inputrec(t_inputrec *ir)
203 memset(ir, 0, (size_t)sizeof(*ir));
204 snew(ir->fepvals, 1);
205 snew(ir->expandedvals, 1);
206 snew(ir->simtempvals, 1);
209 void stupid_fill_block(t_block *grp, int natom, gmx_bool bOneIndexGroup)
215 grp->nalloc_index = 2;
216 snew(grp->index, grp->nalloc_index);
218 grp->index[1] = natom;
223 grp->nalloc_index = natom+1;
224 snew(grp->index, grp->nalloc_index);
225 snew(grp->index, natom+1);
226 for (i = 0; (i <= natom); i++)
234 void stupid_fill_blocka(t_blocka *grp, int natom)
238 grp->nalloc_a = natom;
239 snew(grp->a, grp->nalloc_a);
240 for (i = 0; (i < natom); i++)
246 grp->nalloc_index = natom + 1;
247 snew(grp->index, grp->nalloc_index);
248 for (i = 0; (i <= natom); i++)
255 void copy_blocka(const t_blocka *src, t_blocka *dest)
260 dest->nalloc_index = dest->nr + 1;
261 snew(dest->index, dest->nalloc_index);
262 for (i = 0; i < dest->nr+1; i++)
264 dest->index[i] = src->index[i];
266 dest->nra = src->nra;
267 dest->nalloc_a = dest->nra + 1;
268 snew(dest->a, dest->nalloc_a);
269 for (i = 0; i < dest->nra+1; i++)
271 dest->a[i] = src->a[i];
275 void done_block(t_block *block)
279 block->nalloc_index = 0;
282 void done_blocka(t_blocka *block)
290 block->nalloc_index = 0;
294 void done_atom (t_atoms *at)
302 sfree(at->atomtypeB);
309 void done_atomtypes(t_atomtypes *atype)
312 sfree(atype->radius);
314 sfree(atype->surftens);
315 sfree(atype->atomnumber);
316 sfree(atype->gb_radius);
320 void done_moltype(gmx_moltype_t *molt)
324 done_atom(&molt->atoms);
325 done_block(&molt->cgs);
326 done_blocka(&molt->excls);
328 for (f = 0; f < F_NRE; f++)
330 sfree(molt->ilist[f].iatoms);
331 molt->ilist[f].nalloc = 0;
335 void done_molblock(gmx_molblock_t *molb)
337 if (molb->nposres_xA > 0)
339 molb->nposres_xA = 0;
340 free(molb->posres_xA);
342 if (molb->nposres_xB > 0)
344 molb->nposres_xB = 0;
345 free(molb->posres_xB);
349 void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
355 done_symtab(&mtop->symtab);
358 sfree(mtop->ffparams.functype);
359 sfree(mtop->ffparams.iparams);
361 for (i = 0; i < mtop->nmoltype; i++)
363 done_moltype(&mtop->moltype[i]);
365 sfree(mtop->moltype);
366 for (i = 0; i < mtop->nmolblock; i++)
368 done_molblock(&mtop->molblock[i]);
370 sfree(mtop->molblock);
371 done_block(&mtop->mols);
374 void done_top(t_topology *top)
378 sfree(top->idef.functype);
379 sfree(top->idef.iparams);
380 for (f = 0; f < F_NRE; ++f)
382 sfree(top->idef.il[f].iatoms);
383 top->idef.il[f].iatoms = NULL;
384 top->idef.il[f].nalloc = 0;
387 done_atom (&(top->atoms));
390 done_atomtypes(&(top->atomtypes));
392 done_symtab(&(top->symtab));
393 done_block(&(top->cgs));
394 done_block(&(top->mols));
395 done_blocka(&(top->excls));
398 static void done_pullgrp(t_pullgrp *pgrp)
401 sfree(pgrp->ind_loc);
403 sfree(pgrp->weight_loc);
406 static void done_pull(t_pull *pull)
410 for (i = 0; i < pull->ngrp+1; i++)
412 done_pullgrp(pull->grp);
413 done_pullgrp(pull->dyna);
417 void done_inputrec(t_inputrec *ir)
421 for (m = 0; (m < DIM); m++)
429 sfree(ir->ex[m].phi);
437 sfree(ir->et[m].phi);
441 sfree(ir->opts.nrdf);
442 sfree(ir->opts.ref_t);
443 sfree(ir->opts.annealing);
444 sfree(ir->opts.anneal_npoints);
445 sfree(ir->opts.anneal_time);
446 sfree(ir->opts.anneal_temp);
447 sfree(ir->opts.tau_t);
449 sfree(ir->opts.nFreeze);
450 sfree(ir->opts.QMmethod);
451 sfree(ir->opts.QMbasis);
452 sfree(ir->opts.QMcharge);
453 sfree(ir->opts.QMmult);
455 sfree(ir->opts.CASorbitals);
456 sfree(ir->opts.CASelectrons);
457 sfree(ir->opts.SAon);
458 sfree(ir->opts.SAoff);
459 sfree(ir->opts.SAsteps);
460 sfree(ir->opts.bOPT);
470 static void zero_ekinstate(ekinstate_t *eks)
475 eks->ekinh_old = NULL;
476 eks->ekinscalef_nhc = NULL;
477 eks->ekinscaleh_nhc = NULL;
478 eks->vscale_nhc = NULL;
483 void init_energyhistory(energyhistory_t * enerhist)
487 enerhist->ener_ave = NULL;
488 enerhist->ener_sum = NULL;
489 enerhist->ener_sum_sim = NULL;
490 enerhist->dht = NULL;
492 enerhist->nsteps = 0;
494 enerhist->nsteps_sim = 0;
495 enerhist->nsum_sim = 0;
497 enerhist->dht = NULL;
500 static void done_delta_h_history(delta_h_history_t *dht)
504 for (i = 0; i < dht->nndh; i++)
512 void done_energyhistory(energyhistory_t * enerhist)
514 sfree(enerhist->ener_ave);
515 sfree(enerhist->ener_sum);
516 sfree(enerhist->ener_sum_sim);
518 if (enerhist->dht != NULL)
520 done_delta_h_history(enerhist->dht);
521 sfree(enerhist->dht);
525 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
530 state->nnhpres = nnhpres;
531 state->nhchainlength = nhchainlength;
534 snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
535 snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
536 snew(state->therm_integral, state->ngtc);
537 for (i = 0; i < state->ngtc; i++)
539 for (j = 0; j < state->nhchainlength; j++)
541 state->nosehoover_xi[i*state->nhchainlength + j] = 0.0;
542 state->nosehoover_vxi[i*state->nhchainlength + j] = 0.0;
545 for (i = 0; i < state->ngtc; i++)
547 state->therm_integral[i] = 0.0;
552 state->nosehoover_xi = NULL;
553 state->nosehoover_vxi = NULL;
554 state->therm_integral = NULL;
557 if (state->nnhpres > 0)
559 snew(state->nhpres_xi, state->nhchainlength*nnhpres);
560 snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
561 for (i = 0; i < nnhpres; i++)
563 for (j = 0; j < state->nhchainlength; j++)
565 state->nhpres_xi[i*nhchainlength + j] = 0.0;
566 state->nhpres_vxi[i*nhchainlength + j] = 0.0;
572 state->nhpres_xi = NULL;
573 state->nhpres_vxi = NULL;
578 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
582 state->natoms = natoms;
586 snew(state->lambda, efptNR);
587 for (i = 0; i < efptNR; i++)
589 state->lambda[i] = 0;
592 clear_mat(state->box);
593 clear_mat(state->box_rel);
594 clear_mat(state->boxv);
595 clear_mat(state->pres_prev);
596 clear_mat(state->svir_prev);
597 clear_mat(state->fvir_prev);
598 init_gtc_state(state, ngtc, nnhpres, nhchainlength);
599 state->nalloc = state->natoms;
600 if (state->nalloc > 0)
602 snew(state->x, state->nalloc);
603 snew(state->v, state->nalloc);
613 zero_ekinstate(&state->ekinstate);
615 init_energyhistory(&state->enerhist);
617 init_df_history(&state->dfhist, nlambda, 0);
619 state->ddp_count = 0;
620 state->ddp_count_cg_gl = 0;
622 state->cg_gl_nalloc = 0;
625 void done_state(t_state *state)
627 if (state->nosehoover_xi)
629 sfree(state->nosehoover_xi);
652 state->cg_gl_nalloc = 0;
655 sfree(state->lambda);
659 sfree(state->nosehoover_xi);
660 sfree(state->nosehoover_vxi);
661 sfree(state->therm_integral);
665 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
669 for (d = YY; d <= ZZ; d++)
671 for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
673 /* We need to check if this box component is deformed
674 * or if deformation of another component might cause
675 * changes in this component due to box corrections.
677 if (ir->deform[d][d2] == 0 &&
678 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
679 (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
683 box_rel[d][d2] = b[d][d2]/b[XX][XX];
687 b[d][d2] = b[XX][XX]*box_rel[d][d2];
694 void set_box_rel(t_inputrec *ir, t_state *state)
696 /* Make sure the box obeys the restrictions before we fix the ratios */
697 correct_box(NULL, 0, state->box, NULL);
699 clear_mat(state->box_rel);
701 if (PRESERVE_SHAPE(*ir))
703 do_box_rel(ir, state->box_rel, state->box, TRUE);
707 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
709 if (PRESERVE_SHAPE(*ir))
711 do_box_rel(ir, box_rel, b, FALSE);
715 void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
721 srenew(atoms->atomname, atoms->nr+natom_extra);
722 srenew(atoms->atom, atoms->nr+natom_extra);
723 if (NULL != atoms->pdbinfo)
725 srenew(atoms->pdbinfo, atoms->nr+natom_extra);
727 if (NULL != atoms->atomtype)
729 srenew(atoms->atomtype, atoms->nr+natom_extra);
731 if (NULL != atoms->atomtypeB)
733 srenew(atoms->atomtypeB, atoms->nr+natom_extra);
735 for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
737 atoms->atomname[i] = NULL;
738 memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
739 if (NULL != atoms->pdbinfo)
741 memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
743 if (NULL != atoms->atomtype)
745 atoms->atomtype[i] = NULL;
747 if (NULL != atoms->atomtypeB)
749 atoms->atomtypeB[i] = NULL;
752 atoms->nr += natom_extra;
756 srenew(atoms->resinfo, atoms->nres+nres_extra);
757 for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
759 memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
761 atoms->nres += nres_extra;
765 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
769 snew(atoms->atomname, natoms);
770 atoms->atomtype = NULL;
771 atoms->atomtypeB = NULL;
772 snew(atoms->resinfo, natoms);
773 snew(atoms->atom, natoms);
776 snew(atoms->pdbinfo, natoms);
780 atoms->pdbinfo = NULL;
784 t_atoms *copy_t_atoms(t_atoms *src)
790 init_t_atoms(dst, src->nr, (NULL != src->pdbinfo));
792 if (NULL != src->atomname)
794 snew(dst->atomname, src->nr);
796 if (NULL != src->atomtype)
798 snew(dst->atomtype, src->nr);
800 if (NULL != src->atomtypeB)
802 snew(dst->atomtypeB, src->nr);
804 for (i = 0; (i < src->nr); i++)
806 dst->atom[i] = src->atom[i];
807 if (NULL != src->pdbinfo)
809 dst->pdbinfo[i] = src->pdbinfo[i];
811 if (NULL != src->atomname)
813 dst->atomname[i] = src->atomname[i];
815 if (NULL != src->atomtype)
817 dst->atomtype[i] = src->atomtype[i];
819 if (NULL != src->atomtypeB)
821 dst->atomtypeB[i] = src->atomtypeB[i];
824 dst->nres = src->nres;
825 for (i = 0; (i < src->nres); i++)
827 dst->resinfo[i] = src->resinfo[i];
832 void t_atoms_set_resinfo(t_atoms *atoms, int atom_ind, t_symtab *symtab,
833 const char *resname, int resnr, unsigned char ic,
834 int chainnum, char chainid)
838 ri = &atoms->resinfo[atoms->atom[atom_ind].resind];
839 ri->name = put_symtab(symtab, resname);
843 ri->chainnum = chainnum;
844 ri->chainid = chainid;
847 void free_t_atoms(t_atoms *atoms, gmx_bool bFreeNames)
851 if (bFreeNames && atoms->atomname != NULL)
853 for (i = 0; i < atoms->nr; i++)
855 if (atoms->atomname[i] != NULL)
857 sfree(*atoms->atomname[i]);
858 *atoms->atomname[i] = NULL;
862 if (bFreeNames && atoms->resinfo != NULL)
864 for (i = 0; i < atoms->nres; i++)
866 if (atoms->resinfo[i].name != NULL)
868 sfree(*atoms->resinfo[i].name);
869 *atoms->resinfo[i].name = NULL;
873 if (bFreeNames && atoms->atomtype != NULL)
875 for (i = 0; i < atoms->nr; i++)
877 if (atoms->atomtype[i] != NULL)
879 sfree(*atoms->atomtype[i]);
880 *atoms->atomtype[i] = NULL;
884 if (bFreeNames && atoms->atomtypeB != NULL)
886 for (i = 0; i < atoms->nr; i++)
888 if (atoms->atomtypeB[i] != NULL)
890 sfree(*atoms->atomtypeB[i]);
891 *atoms->atomtypeB[i] = NULL;
895 sfree(atoms->atomname);
896 sfree(atoms->atomtype);
897 sfree(atoms->atomtypeB);
898 sfree(atoms->resinfo);
900 sfree(atoms->pdbinfo);
903 atoms->atomname = NULL;
904 atoms->atomtype = NULL;
905 atoms->atomtypeB = NULL;
906 atoms->resinfo = NULL;
908 atoms->pdbinfo = NULL;
911 real max_cutoff(real cutoff1, real cutoff2)
913 if (cutoff1 == 0 || cutoff2 == 0)
919 return max(cutoff1, cutoff2);
923 extern void init_df_history(df_history_t *dfhist, int nlambda, real wl_delta)
928 dfhist->nlambda = nlambda;
929 dfhist->wl_delta = wl_delta;
930 snew(dfhist->sum_weights, dfhist->nlambda);
931 snew(dfhist->sum_dg, dfhist->nlambda);
932 snew(dfhist->sum_minvar, dfhist->nlambda);
933 snew(dfhist->sum_variance, dfhist->nlambda);
934 snew(dfhist->n_at_lam, dfhist->nlambda);
935 snew(dfhist->wl_histo, dfhist->nlambda);
937 /* allocate transition matrices here */
938 snew(dfhist->Tij, dfhist->nlambda);
939 snew(dfhist->Tij_empirical, dfhist->nlambda);
941 for (i = 0; i < dfhist->nlambda; i++)
943 snew(dfhist->Tij[i], dfhist->nlambda);
944 snew(dfhist->Tij_empirical[i], dfhist->nlambda);
947 snew(dfhist->accum_p, dfhist->nlambda);
948 snew(dfhist->accum_m, dfhist->nlambda);
949 snew(dfhist->accum_p2, dfhist->nlambda);
950 snew(dfhist->accum_m2, dfhist->nlambda);
952 for (i = 0; i < dfhist->nlambda; i++)
954 snew((dfhist->accum_p)[i], dfhist->nlambda);
955 snew((dfhist->accum_m)[i], dfhist->nlambda);
956 snew((dfhist->accum_p2)[i], dfhist->nlambda);
957 snew((dfhist->accum_m2)[i], dfhist->nlambda);
961 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
965 init_df_history(df_dest, df_source->nlambda, df_source->wl_delta);
966 df_dest->nlambda = df_source->nlambda;
967 df_dest->bEquil = df_source->bEquil;
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];
976 df_dest->accum_p[i] = df_source->accum_p[i];
977 df_dest->accum_m[i] = df_source->accum_m[i];
978 df_dest->accum_p2[i] = df_source->accum_p2[i];
979 df_dest->accum_m2[i] = df_source->accum_m2[i];
982 for (i = 0; i < df_dest->nlambda; i++)
984 for (j = 0; j < df_dest->nlambda; 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];