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 "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 static void init_swapstate(swapstate_t *swapstate)
 490 {
 491     int ii, ic;
 492
 493     swapstate->eSwapCoords = 0;
 494     swapstate->nAverage    = 0;
 495
 496     /* Ion/water position swapping */
 497     for (ic = 0; ic < eCompNR; ic++)
 498     {
 499         for (ii = 0; ii < eIonNR; ii++)
 500         {
 501             swapstate->nat_req[ic][ii]        = 0;
 502             swapstate->nat_req_p[ic][ii]      = NULL;
 503             swapstate->inflow_netto[ic][ii]   = 0;
 504             swapstate->inflow_netto_p[ic][ii] = NULL;
 505             swapstate->nat_past[ic][ii]       = NULL;
 506             swapstate->nat_past_p[ic][ii]     = NULL;
 507             swapstate->fluxfromAtoB[ic][ii]   = 0;
 508             swapstate->fluxfromAtoB_p[ic][ii] = NULL;
 509         }
 510     }
 511     swapstate->fluxleak               = NULL;
 512     swapstate->nions                  = 0;
 513     swapstate->comp_from              = NULL;
 514     swapstate->channel_label          = NULL;
 515     swapstate->bFromCpt               = 0;
 516     swapstate->nat[eChan0]            = 0;
 517     swapstate->nat[eChan1]            = 0;
 518     swapstate->xc_old_whole[eChan0]   = NULL;
 519     swapstate->xc_old_whole[eChan1]   = NULL;
 520     swapstate->xc_old_whole_p[eChan0] = NULL;
 521     swapstate->xc_old_whole_p[eChan1] = NULL;
 522 }
 523
 524 void init_energyhistory(energyhistory_t * enerhist)
 525 {
 526     enerhist->nener = 0;
 527
 528     enerhist->ener_ave     = NULL;
 529     enerhist->ener_sum     = NULL;
 530     enerhist->ener_sum_sim = NULL;
 531     enerhist->dht          = NULL;
 532
 533     enerhist->nsteps     = 0;
 534     enerhist->nsum       = 0;
 535     enerhist->nsteps_sim = 0;
 536     enerhist->nsum_sim   = 0;
 537
 538     enerhist->dht = NULL;
 539 }
 540
 541 static void done_delta_h_history(delta_h_history_t *dht)
 542 {
 543     int i;
 544
 545     for (i = 0; i < dht->nndh; i++)
 546     {
 547         sfree(dht->dh[i]);
 548     }
 549     sfree(dht->dh);
 550     sfree(dht->ndh);
 551 }
 552
 553 void done_energyhistory(energyhistory_t * enerhist)
 554 {
 555     sfree(enerhist->ener_ave);
 556     sfree(enerhist->ener_sum);
 557     sfree(enerhist->ener_sum_sim);
 558
 559     if (enerhist->dht != NULL)
 560     {
 561         done_delta_h_history(enerhist->dht);
 562         sfree(enerhist->dht);
 563     }
 564 }
 565
 566 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
 567 {
 568     int i, j;
 569
 570     state->ngtc          = ngtc;
 571     state->nnhpres       = nnhpres;
 572     state->nhchainlength = nhchainlength;
 573     if (state->ngtc > 0)
 574     {
 575         snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
 576         snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
 577         snew(state->therm_integral, state->ngtc);
 578         for (i = 0; i < state->ngtc; i++)
 579         {
 580             for (j = 0; j < state->nhchainlength; j++)
 581             {
 582                 state->nosehoover_xi[i*state->nhchainlength + j]   = 0.0;
 583                 state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
 584             }
 585         }
 586         for (i = 0; i < state->ngtc; i++)
 587         {
 588             state->therm_integral[i]  = 0.0;
 589         }
 590     }
 591     else
 592     {
 593         state->nosehoover_xi  = NULL;
 594         state->nosehoover_vxi = NULL;
 595         state->therm_integral = NULL;
 596     }
 597
 598     if (state->nnhpres > 0)
 599     {
 600         snew(state->nhpres_xi, state->nhchainlength*nnhpres);
 601         snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
 602         for (i = 0; i < nnhpres; i++)
 603         {
 604             for (j = 0; j < state->nhchainlength; j++)
 605             {
 606                 state->nhpres_xi[i*nhchainlength + j]   = 0.0;
 607                 state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
 608             }
 609         }
 610     }
 611     else
 612     {
 613         state->nhpres_xi  = NULL;
 614         state->nhpres_vxi = NULL;
 615     }
 616 }
 617
 618
 619 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
 620 {
 621     int i;
 622
 623     state->natoms = natoms;
 624     state->nrng   = 0;
 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 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
 701 {
 702     int d, d2;
 703
 704     for (d = YY; d <= ZZ; d++)
 705     {
 706         for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
 707         {
 708             /* We need to check if this box component is deformed
 709              * or if deformation of another component might cause
 710              * changes in this component due to box corrections.
 711              */
 712             if (ir->deform[d][d2] == 0 &&
 713                 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
 714                   (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
 715             {
 716                 if (bInit)
 717                 {
 718                     box_rel[d][d2] = b[d][d2]/b[XX][XX];
 719                 }
 720                 else
 721                 {
 722                     b[d][d2] = b[XX][XX]*box_rel[d][d2];
 723                 }
 724             }
 725         }
 726     }
 727 }
 728
 729 void set_box_rel(t_inputrec *ir, t_state *state)
 730 {
 731     /* Make sure the box obeys the restrictions before we fix the ratios */
 732     correct_box(NULL, 0, state->box, NULL);
 733
 734     clear_mat(state->box_rel);
 735
 736     if (PRESERVE_SHAPE(*ir))
 737     {
 738         do_box_rel(ir, state->box_rel, state->box, TRUE);
 739     }
 740 }
 741
 742 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
 743 {
 744     if (PRESERVE_SHAPE(*ir))
 745     {
 746         do_box_rel(ir, box_rel, b, FALSE);
 747     }
 748 }
 749
 750 void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
 751 {
 752     int i;
 753
 754     if (natom_extra > 0)
 755     {
 756         srenew(atoms->atomname, atoms->nr+natom_extra);
 757         srenew(atoms->atom, atoms->nr+natom_extra);
 758         if (NULL != atoms->pdbinfo)
 759         {
 760             srenew(atoms->pdbinfo, atoms->nr+natom_extra);
 761         }
 762         if (NULL != atoms->atomtype)
 763         {
 764             srenew(atoms->atomtype, atoms->nr+natom_extra);
 765         }
 766         if (NULL != atoms->atomtypeB)
 767         {
 768             srenew(atoms->atomtypeB, atoms->nr+natom_extra);
 769         }
 770         for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
 771         {
 772             atoms->atomname[i] = NULL;
 773             memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
 774             if (NULL != atoms->pdbinfo)
 775             {
 776                 memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
 777             }
 778             if (NULL != atoms->atomtype)
 779             {
 780                 atoms->atomtype[i] = NULL;
 781             }
 782             if (NULL != atoms->atomtypeB)
 783             {
 784                 atoms->atomtypeB[i] = NULL;
 785             }
 786         }
 787         atoms->nr += natom_extra;
 788     }
 789     if (nres_extra > 0)
 790     {
 791         srenew(atoms->resinfo, atoms->nres+nres_extra);
 792         for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
 793         {
 794             memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
 795         }
 796         atoms->nres += nres_extra;
 797     }
 798 }
 799
 800 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
 801 {
 802     atoms->nr   = natoms;
 803     atoms->nres = 0;
 804     snew(atoms->atomname, natoms);
 805     atoms->atomtype  = NULL;
 806     atoms->atomtypeB = NULL;
 807     snew(atoms->resinfo, natoms);
 808     snew(atoms->atom, natoms);
 809     if (bPdbinfo)
 810     {
 811         snew(atoms->pdbinfo, natoms);
 812     }
 813     else
 814     {
 815         atoms->pdbinfo = NULL;
 816     }
 817 }
 818
 819 t_atoms *copy_t_atoms(t_atoms *src)
 820 {
 821     t_atoms *dst;
 822     int      i;
 823
 824     snew(dst, 1);
 825     init_t_atoms(dst, src->nr, (NULL != src->pdbinfo));
 826     dst->nr = src->nr;
 827     if (NULL != src->atomname)
 828     {
 829         snew(dst->atomname, src->nr);
 830     }
 831     if (NULL != src->atomtype)
 832     {
 833         snew(dst->atomtype, src->nr);
 834     }
 835     if (NULL != src->atomtypeB)
 836     {
 837         snew(dst->atomtypeB, src->nr);
 838     }
 839     for (i = 0; (i < src->nr); i++)
 840     {
 841         dst->atom[i] = src->atom[i];
 842         if (NULL != src->pdbinfo)
 843         {
 844             dst->pdbinfo[i] = src->pdbinfo[i];
 845         }
 846         if (NULL != src->atomname)
 847         {
 848             dst->atomname[i]  = src->atomname[i];
 849         }
 850         if (NULL != src->atomtype)
 851         {
 852             dst->atomtype[i] = src->atomtype[i];
 853         }
 854         if (NULL != src->atomtypeB)
 855         {
 856             dst->atomtypeB[i] = src->atomtypeB[i];
 857         }
 858     }
 859     dst->nres = src->nres;
 860     for (i = 0; (i < src->nres); i++)
 861     {
 862         dst->resinfo[i] = src->resinfo[i];
 863     }
 864     return dst;
 865 }
 866
 867 void t_atoms_set_resinfo(t_atoms *atoms, int atom_ind, t_symtab *symtab,
 868                          const char *resname, int resnr, unsigned char ic,
 869                          int chainnum, char chainid)
 870 {
 871     t_resinfo *ri;
 872
 873     ri           = &atoms->resinfo[atoms->atom[atom_ind].resind];
 874     ri->name     = put_symtab(symtab, resname);
 875     ri->rtp      = NULL;
 876     ri->nr       = resnr;
 877     ri->ic       = ic;
 878     ri->chainnum = chainnum;
 879     ri->chainid  = chainid;
 880 }
 881
 882 void free_t_atoms(t_atoms *atoms, gmx_bool bFreeNames)
 883 {
 884     int i;
 885
 886     if (bFreeNames && atoms->atomname != NULL)
 887     {
 888         for (i = 0; i < atoms->nr; i++)
 889         {
 890             if (atoms->atomname[i] != NULL)
 891             {
 892                 sfree(*atoms->atomname[i]);
 893                 *atoms->atomname[i] = NULL;
 894             }
 895         }
 896     }
 897     if (bFreeNames && atoms->resinfo != NULL)
 898     {
 899         for (i = 0; i < atoms->nres; i++)
 900         {
 901             if (atoms->resinfo[i].name != NULL)
 902             {
 903                 sfree(*atoms->resinfo[i].name);
 904                 *atoms->resinfo[i].name = NULL;
 905             }
 906         }
 907     }
 908     if (bFreeNames && atoms->atomtype != NULL)
 909     {
 910         for (i = 0; i < atoms->nr; i++)
 911         {
 912             if (atoms->atomtype[i] != NULL)
 913             {
 914                 sfree(*atoms->atomtype[i]);
 915                 *atoms->atomtype[i] = NULL;
 916             }
 917         }
 918     }
 919     if (bFreeNames && atoms->atomtypeB != NULL)
 920     {
 921         for (i = 0; i < atoms->nr; i++)
 922         {
 923             if (atoms->atomtypeB[i] != NULL)
 924             {
 925                 sfree(*atoms->atomtypeB[i]);
 926                 *atoms->atomtypeB[i] = NULL;
 927             }
 928         }
 929     }
 930     sfree(atoms->atomname);
 931     sfree(atoms->atomtype);
 932     sfree(atoms->atomtypeB);
 933     sfree(atoms->resinfo);
 934     sfree(atoms->atom);
 935     sfree(atoms->pdbinfo);
 936     atoms->nr        = 0;
 937     atoms->nres      = 0;
 938     atoms->atomname  = NULL;
 939     atoms->atomtype  = NULL;
 940     atoms->atomtypeB = NULL;
 941     atoms->resinfo   = NULL;
 942     atoms->atom      = NULL;
 943     atoms->pdbinfo   = NULL;
 944 }
 945
 946 real max_cutoff(real cutoff1, real cutoff2)
 947 {
 948     if (cutoff1 == 0 || cutoff2 == 0)
 949     {
 950         return 0;
 951     }
 952     else
 953     {
 954         return max(cutoff1, cutoff2);
 955     }
 956 }
 957
 958 void init_df_history(df_history_t *dfhist, int nlambda)
 959 {
 960     int i;
 961
 962     dfhist->nlambda  = nlambda;
 963     dfhist->bEquil   = 0;
 964     dfhist->wl_delta = 0;
 965
 966     if (nlambda > 0)
 967     {
 968         snew(dfhist->sum_weights, dfhist->nlambda);
 969         snew(dfhist->sum_dg, dfhist->nlambda);
 970         snew(dfhist->sum_minvar, dfhist->nlambda);
 971         snew(dfhist->sum_variance, dfhist->nlambda);
 972         snew(dfhist->n_at_lam, dfhist->nlambda);
 973         snew(dfhist->wl_histo, dfhist->nlambda);
 974
 975         /* allocate transition matrices here */
 976         snew(dfhist->Tij, dfhist->nlambda);
 977         snew(dfhist->Tij_empirical, dfhist->nlambda);
 978
 979         /* allocate accumulators for various transition matrix
 980            free energy methods here */
 981         snew(dfhist->accum_p, dfhist->nlambda);
 982         snew(dfhist->accum_m, dfhist->nlambda);
 983         snew(dfhist->accum_p2, dfhist->nlambda);
 984         snew(dfhist->accum_m2, dfhist->nlambda);
 985
 986         for (i = 0; i < dfhist->nlambda; i++)
 987         {
 988             snew(dfhist->Tij[i], dfhist->nlambda);
 989             snew(dfhist->Tij_empirical[i], dfhist->nlambda);
 990             snew((dfhist->accum_p)[i], dfhist->nlambda);
 991             snew((dfhist->accum_m)[i], dfhist->nlambda);
 992             snew((dfhist->accum_p2)[i], dfhist->nlambda);
 993             snew((dfhist->accum_m2)[i], dfhist->nlambda);
 994         }
 995     }
 996 }
 997
 998 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
 999 {
1000     int i, j;
1001
1002     /* Currently, there should not be any difference in nlambda between the two,
1003        but this is included for completeness for potential later functionality */
1004     df_dest->nlambda  = df_source->nlambda;
1005     df_dest->bEquil   = df_source->bEquil;
1006     df_dest->wl_delta = df_source->wl_delta;
1007
1008     for (i = 0; i < df_dest->nlambda; i++)
1009     {
1010         df_dest->sum_weights[i]  = df_source->sum_weights[i];
1011         df_dest->sum_dg[i]       = df_source->sum_dg[i];
1012         df_dest->sum_minvar[i]   = df_source->sum_minvar[i];
1013         df_dest->sum_variance[i] = df_source->sum_variance[i];
1014         df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
1015         df_dest->wl_histo[i]     = df_source->wl_histo[i];
1016     }
1017
1018     for (i = 0; i < df_dest->nlambda; i++)
1019     {
1020         for (j = 0; j < df_dest->nlambda; j++)
1021         {
1022             df_dest->accum_p[i][j]        = df_source->accum_p[i][j];
1023             df_dest->accum_m[i][j]        = df_source->accum_m[i][j];
1024             df_dest->accum_p2[i][j]       = df_source->accum_p2[i][j];
1025             df_dest->accum_m2[i][j]       = df_source->accum_m2[i][j];
1026             df_dest->Tij[i][j]            = df_source->Tij[i][j];
1027             df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
1028         }
1029     }
1030 }
1031
1032 void done_df_history(df_history_t *dfhist)
1033 {
1034     int i;
1035
1036     if (dfhist->nlambda > 0)
1037     {
1038         sfree(dfhist->n_at_lam);
1039         sfree(dfhist->wl_histo);
1040         sfree(dfhist->sum_weights);
1041         sfree(dfhist->sum_dg);
1042         sfree(dfhist->sum_minvar);
1043         sfree(dfhist->sum_variance);
1044
1045         for (i = 0; i < dfhist->nlambda; i++)
1046         {
1047             sfree(dfhist->Tij[i]);
1048             sfree(dfhist->Tij_empirical[i]);
1049             sfree(dfhist->accum_p[i]);
1050             sfree(dfhist->accum_m[i]);
1051             sfree(dfhist->accum_p2[i]);
1052             sfree(dfhist->accum_m2[i]);
1053         }
1054     }
1055     dfhist->bEquil   = 0;
1056     dfhist->nlambda  = 0;
1057     dfhist->wl_delta = 0;
1058 }