src/gromacs/gmxlib/typedefs.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 /* This file is completely threadsafe - keep it that way! */
  38 #include "typedefs.h"
  39
  40 #ifdef HAVE_CONFIG_H
  41 #include <config.h>
  42 #endif
  43
  44 #include <string.h>
  45
  46 #include "thread_mpi/threads.h"
  47
  48 #include "symtab.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "macros.h"
  52 #include "gromacs/random/random.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 /* The source code in this file should be thread-safe.
  56       Please keep it that way. */
  57
  58
  59
  60 static gmx_bool            bOverAllocDD     = FALSE;
  61 static tMPI_Thread_mutex_t over_alloc_mutex = TMPI_THREAD_MUTEX_INITIALIZER;
  62
  63
  64 void set_over_alloc_dd(gmx_bool set)
  65 {
  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 */
  69     bOverAllocDD = set;
  70     tMPI_Thread_mutex_unlock(&over_alloc_mutex);
  71 }
  72
  73 int over_alloc_dd(int n)
  74 {
  75     if (bOverAllocDD)
  76     {
  77         return OVER_ALLOC_FAC*n + 100;
  78     }
  79     else
  80     {
  81         return n;
  82     }
  83 }
  84
  85 int gmx_int64_to_int(gmx_int64_t step, const char *warn)
  86 {
  87     int i;
  88
  89     i = (int)step;
  90
  91     if (warn != NULL && (step < INT_MIN || step > INT_MAX))
  92     {
  93         fprintf(stderr, "\nWARNING during %s:\n", warn);
  94         fprintf(stderr, "step value ");
  95         fprintf(stderr, "%"GMX_PRId64, step);
  96         fprintf(stderr, " does not fit in int, converted to %d\n\n", i);
  97     }
  98
  99     return i;
 100 }
 101
 102 char *gmx_step_str(gmx_int64_t i, char *buf)
 103 {
 104     sprintf(buf, "%"GMX_PRId64, i);
 105
 106     return buf;
 107 }
 108
 109 void init_block(t_block *block)
 110 {
 111     int i;
 112
 113     block->nr           = 0;
 114     block->nalloc_index = 1;
 115     snew(block->index, block->nalloc_index);
 116     block->index[0]     = 0;
 117 }
 118
 119 void init_blocka(t_blocka *block)
 120 {
 121     int i;
 122
 123     block->nr           = 0;
 124     block->nra          = 0;
 125     block->nalloc_index = 1;
 126     snew(block->index, block->nalloc_index);
 127     block->index[0]     = 0;
 128     block->nalloc_a     = 0;
 129     block->a            = NULL;
 130 }
 131
 132 void init_atom(t_atoms *at)
 133 {
 134     int i;
 135
 136     at->nr        = 0;
 137     at->nres      = 0;
 138     at->atom      = NULL;
 139     at->resinfo   = NULL;
 140     at->atomname  = NULL;
 141     at->atomtype  = NULL;
 142     at->atomtypeB = NULL;
 143     at->pdbinfo   = NULL;
 144 }
 145
 146 void init_atomtypes(t_atomtypes *at)
 147 {
 148     at->nr         = 0;
 149     at->radius     = NULL;
 150     at->vol        = NULL;
 151     at->atomnumber = NULL;
 152     at->gb_radius  = NULL;
 153     at->S_hct      = NULL;
 154 }
 155
 156 void init_groups(gmx_groups_t *groups)
 157 {
 158     int g;
 159
 160     groups->ngrpname = 0;
 161     groups->grpname  = NULL;
 162     for (g = 0; (g < egcNR); g++)
 163     {
 164         groups->grps[g].nm_ind = NULL;
 165         groups->ngrpnr[g]      = 0;
 166         groups->grpnr[g]       = NULL;
 167     }
 168
 169 }
 170
 171 void init_mtop(gmx_mtop_t *mtop)
 172 {
 173     mtop->name         = NULL;
 174     mtop->nmoltype     = 0;
 175     mtop->moltype      = NULL;
 176     mtop->nmolblock    = 0;
 177     mtop->molblock     = NULL;
 178     mtop->maxres_renum = 0;
 179     mtop->maxresnr     = -1;
 180     init_groups(&mtop->groups);
 181     init_block(&mtop->mols);
 182     open_symtab(&mtop->symtab);
 183 }
 184
 185 void init_top(t_topology *top)
 186 {
 187     int i;
 188
 189     top->name = NULL;
 190     init_atom (&(top->atoms));
 191     init_atomtypes(&(top->atomtypes));
 192     init_block(&top->cgs);
 193     init_block(&top->mols);
 194     init_blocka(&top->excls);
 195     open_symtab(&top->symtab);
 196 }
 197
 198 void init_inputrec(t_inputrec *ir)
 199 {
 200     memset(ir, 0, (size_t)sizeof(*ir));
 201     snew(ir->fepvals, 1);
 202     snew(ir->expandedvals, 1);
 203     snew(ir->simtempvals, 1);
 204 }
 205
 206 void stupid_fill_block(t_block *grp, int natom, gmx_bool bOneIndexGroup)
 207 {
 208     int i;
 209
 210     if (bOneIndexGroup)
 211     {
 212         grp->nalloc_index = 2;
 213         snew(grp->index, grp->nalloc_index);
 214         grp->index[0] = 0;
 215         grp->index[1] = natom;
 216         grp->nr       = 1;
 217     }
 218     else
 219     {
 220         grp->nalloc_index = natom+1;
 221         snew(grp->index, grp->nalloc_index);
 222         snew(grp->index, natom+1);
 223         for (i = 0; (i <= natom); i++)
 224         {
 225             grp->index[i] = i;
 226         }
 227         grp->nr = natom;
 228     }
 229 }
 230
 231 void stupid_fill_blocka(t_blocka *grp, int natom)
 232 {
 233     int i;
 234
 235     grp->nalloc_a = natom;
 236     snew(grp->a, grp->nalloc_a);
 237     for (i = 0; (i < natom); i++)
 238     {
 239         grp->a[i] = i;
 240     }
 241     grp->nra = natom;
 242
 243     grp->nalloc_index = natom + 1;
 244     snew(grp->index, grp->nalloc_index);
 245     for (i = 0; (i <= natom); i++)
 246     {
 247         grp->index[i] = i;
 248     }
 249     grp->nr = natom;
 250 }
 251
 252 void copy_blocka(const t_blocka *src, t_blocka *dest)
 253 {
 254     int i;
 255
 256     dest->nr           = src->nr;
 257     dest->nalloc_index = dest->nr + 1;
 258     snew(dest->index, dest->nalloc_index);
 259     for (i = 0; i < dest->nr+1; i++)
 260     {
 261         dest->index[i] = src->index[i];
 262     }
 263     dest->nra      = src->nra;
 264     dest->nalloc_a = dest->nra + 1;
 265     snew(dest->a, dest->nalloc_a);
 266     for (i = 0; i < dest->nra+1; i++)
 267     {
 268         dest->a[i] = src->a[i];
 269     }
 270 }
 271
 272 void done_block(t_block *block)
 273 {
 274     block->nr    = 0;
 275     sfree(block->index);
 276     block->nalloc_index = 0;
 277 }
 278
 279 void done_blocka(t_blocka *block)
 280 {
 281     block->nr    = 0;
 282     block->nra   = 0;
 283     sfree(block->index);
 284     sfree(block->a);
 285     block->index        = NULL;
 286     block->a            = NULL;
 287     block->nalloc_index = 0;
 288     block->nalloc_a     = 0;
 289 }
 290
 291 void done_atom (t_atoms *at)
 292 {
 293     at->nr       = 0;
 294     at->nres     = 0;
 295     sfree(at->atom);
 296     sfree(at->resinfo);
 297     sfree(at->atomname);
 298     sfree(at->atomtype);
 299     sfree(at->atomtypeB);
 300     if (at->pdbinfo)
 301     {
 302         sfree(at->pdbinfo);
 303     }
 304 }
 305
 306 void done_atomtypes(t_atomtypes *atype)
 307 {
 308     atype->nr = 0;
 309     sfree(atype->radius);
 310     sfree(atype->vol);
 311     sfree(atype->surftens);
 312     sfree(atype->atomnumber);
 313     sfree(atype->gb_radius);
 314     sfree(atype->S_hct);
 315 }
 316
 317 void done_moltype(gmx_moltype_t *molt)
 318 {
 319     int f;
 320
 321     done_atom(&molt->atoms);
 322     done_block(&molt->cgs);
 323     done_blocka(&molt->excls);
 324
 325     for (f = 0; f < F_NRE; f++)
 326     {
 327         sfree(molt->ilist[f].iatoms);
 328         molt->ilist[f].nalloc = 0;
 329     }
 330 }
 331
 332 void done_molblock(gmx_molblock_t *molb)
 333 {
 334     if (molb->nposres_xA > 0)
 335     {
 336         molb->nposres_xA = 0;
 337         sfree(molb->posres_xA);
 338     }
 339     if (molb->nposres_xB > 0)
 340     {
 341         molb->nposres_xB = 0;
 342         sfree(molb->posres_xB);
 343     }
 344 }
 345
 346 void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
 347 {
 348     int i;
 349
 350     if (bDoneSymtab)
 351     {
 352         done_symtab(&mtop->symtab);
 353     }
 354
 355     sfree(mtop->ffparams.functype);
 356     sfree(mtop->ffparams.iparams);
 357
 358     for (i = 0; i < mtop->nmoltype; i++)
 359     {
 360         done_moltype(&mtop->moltype[i]);
 361     }
 362     sfree(mtop->moltype);
 363     for (i = 0; i < mtop->nmolblock; i++)
 364     {
 365         done_molblock(&mtop->molblock[i]);
 366     }
 367     sfree(mtop->molblock);
 368     done_block(&mtop->mols);
 369 }
 370
 371 void done_top(t_topology *top)
 372 {
 373     int f;
 374
 375     sfree(top->idef.functype);
 376     sfree(top->idef.iparams);
 377     for (f = 0; f < F_NRE; ++f)
 378     {
 379         sfree(top->idef.il[f].iatoms);
 380         top->idef.il[f].iatoms = NULL;
 381         top->idef.il[f].nalloc = 0;
 382     }
 383
 384     done_atom (&(top->atoms));
 385
 386     /* For GB */
 387     done_atomtypes(&(top->atomtypes));
 388
 389     done_symtab(&(top->symtab));
 390     done_block(&(top->cgs));
 391     done_block(&(top->mols));
 392     done_blocka(&(top->excls));
 393 }
 394
 395 static void done_pull_group(t_pull_group *pgrp)
 396 {
 397     if (pgrp->nat > 0)
 398     {
 399         sfree(pgrp->ind);
 400         sfree(pgrp->ind_loc);
 401         sfree(pgrp->weight);
 402         sfree(pgrp->weight_loc);
 403     }
 404 }
 405
 406 static void done_pull(t_pull *pull)
 407 {
 408     int i;
 409
 410     for (i = 0; i < pull->ngroup+1; i++)
 411     {
 412         done_pull_group(pull->group);
 413         done_pull_group(pull->dyna);
 414     }
 415 }
 416
 417 void done_inputrec(t_inputrec *ir)
 418 {
 419     int m;
 420
 421     for (m = 0; (m < DIM); m++)
 422     {
 423         if (ir->ex[m].a)
 424         {
 425             sfree(ir->ex[m].a);
 426         }
 427         if (ir->ex[m].phi)
 428         {
 429             sfree(ir->ex[m].phi);
 430         }
 431         if (ir->et[m].a)
 432         {
 433             sfree(ir->et[m].a);
 434         }
 435         if (ir->et[m].phi)
 436         {
 437             sfree(ir->et[m].phi);
 438         }
 439     }
 440
 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);
 448     sfree(ir->opts.acc);
 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);
 454     sfree(ir->opts.bSH);
 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);
 461     sfree(ir->opts.bTS);
 462
 463     if (ir->pull)
 464     {
 465         done_pull(ir->pull);
 466         sfree(ir->pull);
 467     }
 468 }
 469
 470 static void zero_history(history_t *hist)
 471 {
 472     hist->disre_initf  = 0;
 473     hist->ndisrepairs  = 0;
 474     hist->disre_rm3tav = NULL;
 475     hist->orire_initf  = 0;
 476     hist->norire_Dtav  = 0;
 477     hist->orire_Dtav   = NULL;
 478 }
 479
 480 static void zero_ekinstate(ekinstate_t *eks)
 481 {
 482     eks->ekin_n         = 0;
 483     eks->ekinh          = NULL;
 484     eks->ekinf          = NULL;
 485     eks->ekinh_old      = NULL;
 486     eks->ekinscalef_nhc = NULL;
 487     eks->ekinscaleh_nhc = NULL;
 488     eks->vscale_nhc     = NULL;
 489     eks->dekindl        = 0;
 490     eks->mvcos          = 0;
 491 }
 492
 493 static void init_swapstate(swapstate_t *swapstate)
 494 {
 495     int ii, ic;
 496
 497     swapstate->eSwapCoords = 0;
 498     swapstate->nAverage    = 0;
 499
 500     /* Ion/water position swapping */
 501     for (ic = 0; ic < eCompNR; ic++)
 502     {
 503         for (ii = 0; ii < eIonNR; ii++)
 504         {
 505             swapstate->nat_req[ic][ii]        = 0;
 506             swapstate->nat_req_p[ic][ii]      = NULL;
 507             swapstate->inflow_netto[ic][ii]   = 0;
 508             swapstate->inflow_netto_p[ic][ii] = NULL;
 509             swapstate->nat_past[ic][ii]       = NULL;
 510             swapstate->nat_past_p[ic][ii]     = NULL;
 511             swapstate->fluxfromAtoB[ic][ii]   = 0;
 512             swapstate->fluxfromAtoB_p[ic][ii] = NULL;
 513         }
 514     }
 515     swapstate->fluxleak               = NULL;
 516     swapstate->nions                  = 0;
 517     swapstate->comp_from              = NULL;
 518     swapstate->channel_label          = NULL;
 519     swapstate->bFromCpt               = 0;
 520     swapstate->nat[eChan0]            = 0;
 521     swapstate->nat[eChan1]            = 0;
 522     swapstate->xc_old_whole[eChan0]   = NULL;
 523     swapstate->xc_old_whole[eChan1]   = NULL;
 524     swapstate->xc_old_whole_p[eChan0] = NULL;
 525     swapstate->xc_old_whole_p[eChan1] = NULL;
 526 }
 527
 528 void init_energyhistory(energyhistory_t * enerhist)
 529 {
 530     enerhist->nener = 0;
 531
 532     enerhist->ener_ave     = NULL;
 533     enerhist->ener_sum     = NULL;
 534     enerhist->ener_sum_sim = NULL;
 535     enerhist->dht          = NULL;
 536
 537     enerhist->nsteps     = 0;
 538     enerhist->nsum       = 0;
 539     enerhist->nsteps_sim = 0;
 540     enerhist->nsum_sim   = 0;
 541
 542     enerhist->dht = NULL;
 543 }
 544
 545 static void done_delta_h_history(delta_h_history_t *dht)
 546 {
 547     int i;
 548
 549     for (i = 0; i < dht->nndh; i++)
 550     {
 551         sfree(dht->dh[i]);
 552     }
 553     sfree(dht->dh);
 554     sfree(dht->ndh);
 555 }
 556
 557 void done_energyhistory(energyhistory_t * enerhist)
 558 {
 559     sfree(enerhist->ener_ave);
 560     sfree(enerhist->ener_sum);
 561     sfree(enerhist->ener_sum_sim);
 562
 563     if (enerhist->dht != NULL)
 564     {
 565         done_delta_h_history(enerhist->dht);
 566         sfree(enerhist->dht);
 567     }
 568 }
 569
 570 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
 571 {
 572     int i, j;
 573
 574     state->ngtc          = ngtc;
 575     state->nnhpres       = nnhpres;
 576     state->nhchainlength = nhchainlength;
 577     if (state->ngtc > 0)
 578     {
 579         snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
 580         snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
 581         snew(state->therm_integral, state->ngtc);
 582         for (i = 0; i < state->ngtc; i++)
 583         {
 584             for (j = 0; j < state->nhchainlength; j++)
 585             {
 586                 state->nosehoover_xi[i*state->nhchainlength + j]   = 0.0;
 587                 state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
 588             }
 589         }
 590         for (i = 0; i < state->ngtc; i++)
 591         {
 592             state->therm_integral[i]  = 0.0;
 593         }
 594     }
 595     else
 596     {
 597         state->nosehoover_xi  = NULL;
 598         state->nosehoover_vxi = NULL;
 599         state->therm_integral = NULL;
 600     }
 601
 602     if (state->nnhpres > 0)
 603     {
 604         snew(state->nhpres_xi, state->nhchainlength*nnhpres);
 605         snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
 606         for (i = 0; i < nnhpres; i++)
 607         {
 608             for (j = 0; j < state->nhchainlength; j++)
 609             {
 610                 state->nhpres_xi[i*nhchainlength + j]   = 0.0;
 611                 state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
 612             }
 613         }
 614     }
 615     else
 616     {
 617         state->nhpres_xi  = NULL;
 618         state->nhpres_vxi = NULL;
 619     }
 620 }
 621
 622
 623 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
 624 {
 625     int i;
 626
 627     state->natoms = natoms;
 628     state->flags  = 0;
 629     state->lambda = 0;
 630     snew(state->lambda, efptNR);
 631     for (i = 0; i < efptNR; i++)
 632     {
 633         state->lambda[i] = 0;
 634     }
 635     state->veta   = 0;
 636     clear_mat(state->box);
 637     clear_mat(state->box_rel);
 638     clear_mat(state->boxv);
 639     clear_mat(state->pres_prev);
 640     clear_mat(state->svir_prev);
 641     clear_mat(state->fvir_prev);
 642     init_gtc_state(state, ngtc, nnhpres, nhchainlength);
 643     state->nalloc = state->natoms;
 644     if (state->nalloc > 0)
 645     {
 646         snew(state->x, state->nalloc);
 647         snew(state->v, state->nalloc);
 648     }
 649     else
 650     {
 651         state->x = NULL;
 652         state->v = NULL;
 653     }
 654     state->sd_X = NULL;
 655     state->cg_p = NULL;
 656     zero_history(&state->hist);
 657     zero_ekinstate(&state->ekinstate);
 658     init_energyhistory(&state->enerhist);
 659     init_df_history(&state->dfhist, nlambda);
 660     init_swapstate(&state->swapstate);
 661     state->ddp_count       = 0;
 662     state->ddp_count_cg_gl = 0;
 663     state->cg_gl           = NULL;
 664     state->cg_gl_nalloc    = 0;
 665 }
 666
 667 void done_state(t_state *state)
 668 {
 669     if (state->x)
 670     {
 671         sfree(state->x);
 672     }
 673     if (state->v)
 674     {
 675         sfree(state->v);
 676     }
 677     if (state->sd_X)
 678     {
 679         sfree(state->sd_X);
 680     }
 681     if (state->cg_p)
 682     {
 683         sfree(state->cg_p);
 684     }
 685     state->nalloc = 0;
 686     if (state->cg_gl)
 687     {
 688         sfree(state->cg_gl);
 689     }
 690     state->cg_gl_nalloc = 0;
 691     if (state->lambda)
 692     {
 693         sfree(state->lambda);
 694     }
 695     if (state->ngtc > 0)
 696     {
 697         sfree(state->nosehoover_xi);
 698         sfree(state->nosehoover_vxi);
 699         sfree(state->therm_integral);
 700     }
 701 }
 702
 703 t_state *serial_init_local_state(t_state *state_global)
 704 {
 705     int      i;
 706     t_state *state_local;
 707
 708     snew(state_local, 1);
 709
 710     /* Copy all the contents */
 711     *state_local = *state_global;
 712     snew(state_local->lambda, efptNR);
 713     /* local storage for lambda */
 714     for (i = 0; i < efptNR; i++)
 715     {
 716         state_local->lambda[i] = state_global->lambda[i];
 717     }
 718
 719     return state_local;
 720 }
 721
 722 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
 723 {
 724     int d, d2;
 725
 726     for (d = YY; d <= ZZ; d++)
 727     {
 728         for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
 729         {
 730             /* We need to check if this box component is deformed
 731              * or if deformation of another component might cause
 732              * changes in this component due to box corrections.
 733              */
 734             if (ir->deform[d][d2] == 0 &&
 735                 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
 736                   (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
 737             {
 738                 if (bInit)
 739                 {
 740                     box_rel[d][d2] = b[d][d2]/b[XX][XX];
 741                 }
 742                 else
 743                 {
 744                     b[d][d2] = b[XX][XX]*box_rel[d][d2];
 745                 }
 746             }
 747         }
 748     }
 749 }
 750
 751 void set_box_rel(t_inputrec *ir, t_state *state)
 752 {
 753     /* Make sure the box obeys the restrictions before we fix the ratios */
 754     correct_box(NULL, 0, state->box, NULL);
 755
 756     clear_mat(state->box_rel);
 757
 758     if (PRESERVE_SHAPE(*ir))
 759     {
 760         do_box_rel(ir, state->box_rel, state->box, TRUE);
 761     }
 762 }
 763
 764 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
 765 {
 766     if (PRESERVE_SHAPE(*ir))
 767     {
 768         do_box_rel(ir, box_rel, b, FALSE);
 769     }
 770 }
 771
 772 void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
 773 {
 774     int i;
 775
 776     if (natom_extra > 0)
 777     {
 778         srenew(atoms->atomname, atoms->nr+natom_extra);
 779         srenew(atoms->atom, atoms->nr+natom_extra);
 780         if (NULL != atoms->pdbinfo)
 781         {
 782             srenew(atoms->pdbinfo, atoms->nr+natom_extra);
 783         }
 784         if (NULL != atoms->atomtype)
 785         {
 786             srenew(atoms->atomtype, atoms->nr+natom_extra);
 787         }
 788         if (NULL != atoms->atomtypeB)
 789         {
 790             srenew(atoms->atomtypeB, atoms->nr+natom_extra);
 791         }
 792         for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
 793         {
 794             atoms->atomname[i] = NULL;
 795             memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
 796             if (NULL != atoms->pdbinfo)
 797             {
 798                 memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
 799             }
 800             if (NULL != atoms->atomtype)
 801             {
 802                 atoms->atomtype[i] = NULL;
 803             }
 804             if (NULL != atoms->atomtypeB)
 805             {
 806                 atoms->atomtypeB[i] = NULL;
 807             }
 808         }
 809         atoms->nr += natom_extra;
 810     }
 811     if (nres_extra > 0)
 812     {
 813         srenew(atoms->resinfo, atoms->nres+nres_extra);
 814         for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
 815         {
 816             memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
 817         }
 818         atoms->nres += nres_extra;
 819     }
 820 }
 821
 822 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
 823 {
 824     atoms->nr   = natoms;
 825     atoms->nres = 0;
 826     snew(atoms->atomname, natoms);
 827     atoms->atomtype  = NULL;
 828     atoms->atomtypeB = NULL;
 829     snew(atoms->resinfo, natoms);
 830     snew(atoms->atom, natoms);
 831     if (bPdbinfo)
 832     {
 833         snew(atoms->pdbinfo, natoms);
 834     }
 835     else
 836     {
 837         atoms->pdbinfo = NULL;
 838     }
 839 }
 840
 841 t_atoms *copy_t_atoms(t_atoms *src)
 842 {
 843     t_atoms *dst;
 844     int      i;
 845
 846     snew(dst, 1);
 847     init_t_atoms(dst, src->nr, (NULL != src->pdbinfo));
 848     dst->nr = src->nr;
 849     if (NULL != src->atomname)
 850     {
 851         snew(dst->atomname, src->nr);
 852     }
 853     if (NULL != src->atomtype)
 854     {
 855         snew(dst->atomtype, src->nr);
 856     }
 857     if (NULL != src->atomtypeB)
 858     {
 859         snew(dst->atomtypeB, src->nr);
 860     }
 861     for (i = 0; (i < src->nr); i++)
 862     {
 863         dst->atom[i] = src->atom[i];
 864         if (NULL != src->pdbinfo)
 865         {
 866             dst->pdbinfo[i] = src->pdbinfo[i];
 867         }
 868         if (NULL != src->atomname)
 869         {
 870             dst->atomname[i]  = src->atomname[i];
 871         }
 872         if (NULL != src->atomtype)
 873         {
 874             dst->atomtype[i] = src->atomtype[i];
 875         }
 876         if (NULL != src->atomtypeB)
 877         {
 878             dst->atomtypeB[i] = src->atomtypeB[i];
 879         }
 880     }
 881     dst->nres = src->nres;
 882     for (i = 0; (i < src->nres); i++)
 883     {
 884         dst->resinfo[i] = src->resinfo[i];
 885     }
 886     return dst;
 887 }
 888
 889 void t_atoms_set_resinfo(t_atoms *atoms, int atom_ind, t_symtab *symtab,
 890                          const char *resname, int resnr, unsigned char ic,
 891                          int chainnum, char chainid)
 892 {
 893     t_resinfo *ri;
 894
 895     ri           = &atoms->resinfo[atoms->atom[atom_ind].resind];
 896     ri->name     = put_symtab(symtab, resname);
 897     ri->rtp      = NULL;
 898     ri->nr       = resnr;
 899     ri->ic       = ic;
 900     ri->chainnum = chainnum;
 901     ri->chainid  = chainid;
 902 }
 903
 904 void free_t_atoms(t_atoms *atoms, gmx_bool bFreeNames)
 905 {
 906     int i;
 907
 908     if (bFreeNames && atoms->atomname != NULL)
 909     {
 910         for (i = 0; i < atoms->nr; i++)
 911         {
 912             if (atoms->atomname[i] != NULL)
 913             {
 914                 sfree(*atoms->atomname[i]);
 915                 *atoms->atomname[i] = NULL;
 916             }
 917         }
 918     }
 919     if (bFreeNames && atoms->resinfo != NULL)
 920     {
 921         for (i = 0; i < atoms->nres; i++)
 922         {
 923             if (atoms->resinfo[i].name != NULL)
 924             {
 925                 sfree(*atoms->resinfo[i].name);
 926                 *atoms->resinfo[i].name = NULL;
 927             }
 928         }
 929     }
 930     if (bFreeNames && atoms->atomtype != NULL)
 931     {
 932         for (i = 0; i < atoms->nr; i++)
 933         {
 934             if (atoms->atomtype[i] != NULL)
 935             {
 936                 sfree(*atoms->atomtype[i]);
 937                 *atoms->atomtype[i] = NULL;
 938             }
 939         }
 940     }
 941     if (bFreeNames && atoms->atomtypeB != NULL)
 942     {
 943         for (i = 0; i < atoms->nr; i++)
 944         {
 945             if (atoms->atomtypeB[i] != NULL)
 946             {
 947                 sfree(*atoms->atomtypeB[i]);
 948                 *atoms->atomtypeB[i] = NULL;
 949             }
 950         }
 951     }
 952     sfree(atoms->atomname);
 953     sfree(atoms->atomtype);
 954     sfree(atoms->atomtypeB);
 955     sfree(atoms->resinfo);
 956     sfree(atoms->atom);
 957     sfree(atoms->pdbinfo);
 958     atoms->nr        = 0;
 959     atoms->nres      = 0;
 960     atoms->atomname  = NULL;
 961     atoms->atomtype  = NULL;
 962     atoms->atomtypeB = NULL;
 963     atoms->resinfo   = NULL;
 964     atoms->atom      = NULL;
 965     atoms->pdbinfo   = NULL;
 966 }
 967
 968 real max_cutoff(real cutoff1, real cutoff2)
 969 {
 970     if (cutoff1 == 0 || cutoff2 == 0)
 971     {
 972         return 0;
 973     }
 974     else
 975     {
 976         return max(cutoff1, cutoff2);
 977     }
 978 }
 979
 980 void init_df_history(df_history_t *dfhist, int nlambda)
 981 {
 982     int i;
 983
 984     dfhist->nlambda  = nlambda;
 985     dfhist->bEquil   = 0;
 986     dfhist->wl_delta = 0;
 987
 988     if (nlambda > 0)
 989     {
 990         snew(dfhist->sum_weights, dfhist->nlambda);
 991         snew(dfhist->sum_dg, dfhist->nlambda);
 992         snew(dfhist->sum_minvar, dfhist->nlambda);
 993         snew(dfhist->sum_variance, dfhist->nlambda);
 994         snew(dfhist->n_at_lam, dfhist->nlambda);
 995         snew(dfhist->wl_histo, dfhist->nlambda);
 996
 997         /* allocate transition matrices here */
 998         snew(dfhist->Tij, dfhist->nlambda);
 999         snew(dfhist->Tij_empirical, dfhist->nlambda);
1000
1001         /* allocate accumulators for various transition matrix
1002            free energy methods here */
1003         snew(dfhist->accum_p, dfhist->nlambda);
1004         snew(dfhist->accum_m, dfhist->nlambda);
1005         snew(dfhist->accum_p2, dfhist->nlambda);
1006         snew(dfhist->accum_m2, dfhist->nlambda);
1007
1008         for (i = 0; i < dfhist->nlambda; i++)
1009         {
1010             snew(dfhist->Tij[i], dfhist->nlambda);
1011             snew(dfhist->Tij_empirical[i], dfhist->nlambda);
1012             snew((dfhist->accum_p)[i], dfhist->nlambda);
1013             snew((dfhist->accum_m)[i], dfhist->nlambda);
1014             snew((dfhist->accum_p2)[i], dfhist->nlambda);
1015             snew((dfhist->accum_m2)[i], dfhist->nlambda);
1016         }
1017     }
1018 }
1019
1020 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
1021 {
1022     int i, j;
1023
1024     /* Currently, there should not be any difference in nlambda between the two,
1025        but this is included for completeness for potential later functionality */
1026     df_dest->nlambda  = df_source->nlambda;
1027     df_dest->bEquil   = df_source->bEquil;
1028     df_dest->wl_delta = df_source->wl_delta;
1029
1030     for (i = 0; i < df_dest->nlambda; i++)
1031     {
1032         df_dest->sum_weights[i]  = df_source->sum_weights[i];
1033         df_dest->sum_dg[i]       = df_source->sum_dg[i];
1034         df_dest->sum_minvar[i]   = df_source->sum_minvar[i];
1035         df_dest->sum_variance[i] = df_source->sum_variance[i];
1036         df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
1037         df_dest->wl_histo[i]     = df_source->wl_histo[i];
1038     }
1039
1040     for (i = 0; i < df_dest->nlambda; i++)
1041     {
1042         for (j = 0; j < df_dest->nlambda; j++)
1043         {
1044             df_dest->accum_p[i][j]        = df_source->accum_p[i][j];
1045             df_dest->accum_m[i][j]        = df_source->accum_m[i][j];
1046             df_dest->accum_p2[i][j]       = df_source->accum_p2[i][j];
1047             df_dest->accum_m2[i][j]       = df_source->accum_m2[i][j];
1048             df_dest->Tij[i][j]            = df_source->Tij[i][j];
1049             df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
1050         }
1051     }
1052 }
1053
1054 void done_df_history(df_history_t *dfhist)
1055 {
1056     int i;
1057
1058     if (dfhist->nlambda > 0)
1059     {
1060         sfree(dfhist->n_at_lam);
1061         sfree(dfhist->wl_histo);
1062         sfree(dfhist->sum_weights);
1063         sfree(dfhist->sum_dg);
1064         sfree(dfhist->sum_minvar);
1065         sfree(dfhist->sum_variance);
1066
1067         for (i = 0; i < dfhist->nlambda; i++)
1068         {
1069             sfree(dfhist->Tij[i]);
1070             sfree(dfhist->Tij_empirical[i]);
1071             sfree(dfhist->accum_p[i]);
1072             sfree(dfhist->accum_m[i]);
1073             sfree(dfhist->accum_p2[i]);
1074             sfree(dfhist->accum_m2[i]);
1075         }
1076     }
1077     dfhist->bEquil   = 0;
1078     dfhist->nlambda  = 0;
1079     dfhist->wl_delta = 0;
1080 }