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