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