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