src/gromacs/gmxlib/typedefs.c

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