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_pull_group(t_pull_group *pgrp)
 391 {
 392     if (pgrp->nat > 0)
 393     {
 394         sfree(pgrp->ind);
 395         sfree(pgrp->ind_loc);
 396         sfree(pgrp->weight);
 397         sfree(pgrp->weight_loc);
 398     }
 399 }
 400
 401 static void done_pull(t_pull *pull)
 402 {
 403     int i;
 404
 405     for (i = 0; i < pull->ngroup+1; i++)
 406     {
 407         done_pull_group(pull->group);
 408         done_pull_group(pull->dyna);
 409     }
 410 }
 411
 412 void done_inputrec(t_inputrec *ir)
 413 {
 414     int m;
 415
 416     for (m = 0; (m < DIM); m++)
 417     {
 418         if (ir->ex[m].a)
 419         {
 420             sfree(ir->ex[m].a);
 421         }
 422         if (ir->ex[m].phi)
 423         {
 424             sfree(ir->ex[m].phi);
 425         }
 426         if (ir->et[m].a)
 427         {
 428             sfree(ir->et[m].a);
 429         }
 430         if (ir->et[m].phi)
 431         {
 432             sfree(ir->et[m].phi);
 433         }
 434     }
 435
 436     sfree(ir->opts.nrdf);
 437     sfree(ir->opts.ref_t);
 438     sfree(ir->opts.annealing);
 439     sfree(ir->opts.anneal_npoints);
 440     sfree(ir->opts.anneal_time);
 441     sfree(ir->opts.anneal_temp);
 442     sfree(ir->opts.tau_t);
 443     sfree(ir->opts.acc);
 444     sfree(ir->opts.nFreeze);
 445     sfree(ir->opts.QMmethod);
 446     sfree(ir->opts.QMbasis);
 447     sfree(ir->opts.QMcharge);
 448     sfree(ir->opts.QMmult);
 449     sfree(ir->opts.bSH);
 450     sfree(ir->opts.CASorbitals);
 451     sfree(ir->opts.CASelectrons);
 452     sfree(ir->opts.SAon);
 453     sfree(ir->opts.SAoff);
 454     sfree(ir->opts.SAsteps);
 455     sfree(ir->opts.bOPT);
 456     sfree(ir->opts.bTS);
 457
 458     if (ir->pull)
 459     {
 460         done_pull(ir->pull);
 461         sfree(ir->pull);
 462     }
 463 }
 464
 465 static void zero_history(history_t *hist)
 466 {
 467     hist->disre_initf = 0;
 468     hist->ndisrepairs = 0;
 469     hist->disre_rm3tav = NULL;
 470     hist->orire_initf = 0;
 471     hist->norire_Dtav = 0;
 472     hist->orire_Dtav = NULL;
 473 }
 474
 475 static void zero_ekinstate(ekinstate_t *eks)
 476 {
 477     eks->ekin_n         = 0;
 478     eks->ekinh          = NULL;
 479     eks->ekinf          = NULL;
 480     eks->ekinh_old      = NULL;
 481     eks->ekinscalef_nhc = NULL;
 482     eks->ekinscaleh_nhc = NULL;
 483     eks->vscale_nhc     = NULL;
 484     eks->dekindl        = 0;
 485     eks->mvcos          = 0;
 486 }
 487
 488 void init_energyhistory(energyhistory_t * enerhist)
 489 {
 490     enerhist->nener = 0;
 491
 492     enerhist->ener_ave     = NULL;
 493     enerhist->ener_sum     = NULL;
 494     enerhist->ener_sum_sim = NULL;
 495     enerhist->dht          = NULL;
 496
 497     enerhist->nsteps     = 0;
 498     enerhist->nsum       = 0;
 499     enerhist->nsteps_sim = 0;
 500     enerhist->nsum_sim   = 0;
 501
 502     enerhist->dht = NULL;
 503 }
 504
 505 static void done_delta_h_history(delta_h_history_t *dht)
 506 {
 507     int i;
 508
 509     for (i = 0; i < dht->nndh; i++)
 510     {
 511         sfree(dht->dh[i]);
 512     }
 513     sfree(dht->dh);
 514     sfree(dht->ndh);
 515 }
 516
 517 void done_energyhistory(energyhistory_t * enerhist)
 518 {
 519     sfree(enerhist->ener_ave);
 520     sfree(enerhist->ener_sum);
 521     sfree(enerhist->ener_sum_sim);
 522
 523     if (enerhist->dht != NULL)
 524     {
 525         done_delta_h_history(enerhist->dht);
 526         sfree(enerhist->dht);
 527     }
 528 }
 529
 530 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
 531 {
 532     int i, j;
 533
 534     state->ngtc          = ngtc;
 535     state->nnhpres       = nnhpres;
 536     state->nhchainlength = nhchainlength;
 537     if (state->ngtc > 0)
 538     {
 539         snew(state->nosehoover_xi, state->nhchainlength*state->ngtc);
 540         snew(state->nosehoover_vxi, state->nhchainlength*state->ngtc);
 541         snew(state->therm_integral, state->ngtc);
 542         for (i = 0; i < state->ngtc; i++)
 543         {
 544             for (j = 0; j < state->nhchainlength; j++)
 545             {
 546                 state->nosehoover_xi[i*state->nhchainlength + j]   = 0.0;
 547                 state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
 548             }
 549         }
 550         for (i = 0; i < state->ngtc; i++)
 551         {
 552             state->therm_integral[i]  = 0.0;
 553         }
 554     }
 555     else
 556     {
 557         state->nosehoover_xi  = NULL;
 558         state->nosehoover_vxi = NULL;
 559         state->therm_integral = NULL;
 560     }
 561
 562     if (state->nnhpres > 0)
 563     {
 564         snew(state->nhpres_xi, state->nhchainlength*nnhpres);
 565         snew(state->nhpres_vxi, state->nhchainlength*nnhpres);
 566         for (i = 0; i < nnhpres; i++)
 567         {
 568             for (j = 0; j < state->nhchainlength; j++)
 569             {
 570                 state->nhpres_xi[i*nhchainlength + j]   = 0.0;
 571                 state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
 572             }
 573         }
 574     }
 575     else
 576     {
 577         state->nhpres_xi  = NULL;
 578         state->nhpres_vxi = NULL;
 579     }
 580 }
 581
 582
 583 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
 584 {
 585     int i;
 586
 587     state->natoms = natoms;
 588     state->nrng   = 0;
 589     state->flags  = 0;
 590     state->lambda = 0;
 591     snew(state->lambda, efptNR);
 592     for (i = 0; i < efptNR; i++)
 593     {
 594         state->lambda[i] = 0;
 595     }
 596     state->veta   = 0;
 597     clear_mat(state->box);
 598     clear_mat(state->box_rel);
 599     clear_mat(state->boxv);
 600     clear_mat(state->pres_prev);
 601     clear_mat(state->svir_prev);
 602     clear_mat(state->fvir_prev);
 603     init_gtc_state(state, ngtc, nnhpres, nhchainlength);
 604     state->nalloc = state->natoms;
 605     if (state->nalloc > 0)
 606     {
 607         snew(state->x, state->nalloc);
 608         snew(state->v, state->nalloc);
 609     }
 610     else
 611     {
 612         state->x = NULL;
 613         state->v = NULL;
 614     }
 615     state->sd_X = NULL;
 616     state->cg_p = NULL;
 617     zero_history(&state->hist);
 618     zero_ekinstate(&state->ekinstate);
 619     init_energyhistory(&state->enerhist);
 620     init_df_history(&state->dfhist,nlambda);
 621     state->ddp_count       = 0;
 622     state->ddp_count_cg_gl = 0;
 623     state->cg_gl           = NULL;
 624     state->cg_gl_nalloc    = 0;
 625 }
 626
 627 void done_state(t_state *state)
 628 {
 629     if (state->x)
 630     {
 631         sfree(state->x);
 632     }
 633     if (state->v)
 634     {
 635         sfree(state->v);
 636     }
 637     if (state->sd_X)
 638     {
 639         sfree(state->sd_X);
 640     }
 641     if (state->cg_p)
 642     {
 643         sfree(state->cg_p);
 644     }
 645     state->nalloc = 0;
 646     if (state->cg_gl)
 647     {
 648         sfree(state->cg_gl);
 649     }
 650     state->cg_gl_nalloc = 0;
 651     if (state->lambda)
 652     {
 653         sfree(state->lambda);
 654     }
 655     if (state->ngtc > 0)
 656     {
 657         sfree(state->nosehoover_xi);
 658         sfree(state->nosehoover_vxi);
 659         sfree(state->therm_integral);
 660     }
 661 }
 662
 663 static void do_box_rel(t_inputrec *ir, matrix box_rel, matrix b, gmx_bool bInit)
 664 {
 665     int d, d2;
 666
 667     for (d = YY; d <= ZZ; d++)
 668     {
 669         for (d2 = XX; d2 <= (ir->epct == epctSEMIISOTROPIC ? YY : ZZ); d2++)
 670         {
 671             /* We need to check if this box component is deformed
 672              * or if deformation of another component might cause
 673              * changes in this component due to box corrections.
 674              */
 675             if (ir->deform[d][d2] == 0 &&
 676                 !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
 677                   (b[YY][d2] != 0 || ir->deform[YY][d2] != 0)))
 678             {
 679                 if (bInit)
 680                 {
 681                     box_rel[d][d2] = b[d][d2]/b[XX][XX];
 682                 }
 683                 else
 684                 {
 685                     b[d][d2] = b[XX][XX]*box_rel[d][d2];
 686                 }
 687             }
 688         }
 689     }
 690 }
 691
 692 void set_box_rel(t_inputrec *ir, t_state *state)
 693 {
 694     /* Make sure the box obeys the restrictions before we fix the ratios */
 695     correct_box(NULL, 0, state->box, NULL);
 696
 697     clear_mat(state->box_rel);
 698
 699     if (PRESERVE_SHAPE(*ir))
 700     {
 701         do_box_rel(ir, state->box_rel, state->box, TRUE);
 702     }
 703 }
 704
 705 void preserve_box_shape(t_inputrec *ir, matrix box_rel, matrix b)
 706 {
 707     if (PRESERVE_SHAPE(*ir))
 708     {
 709         do_box_rel(ir, box_rel, b, FALSE);
 710     }
 711 }
 712
 713 void add_t_atoms(t_atoms *atoms, int natom_extra, int nres_extra)
 714 {
 715     int i;
 716
 717     if (natom_extra > 0)
 718     {
 719         srenew(atoms->atomname, atoms->nr+natom_extra);
 720         srenew(atoms->atom, atoms->nr+natom_extra);
 721         if (NULL != atoms->pdbinfo)
 722         {
 723             srenew(atoms->pdbinfo, atoms->nr+natom_extra);
 724         }
 725         if (NULL != atoms->atomtype)
 726         {
 727             srenew(atoms->atomtype, atoms->nr+natom_extra);
 728         }
 729         if (NULL != atoms->atomtypeB)
 730         {
 731             srenew(atoms->atomtypeB, atoms->nr+natom_extra);
 732         }
 733         for (i = atoms->nr; (i < atoms->nr+natom_extra); i++)
 734         {
 735             atoms->atomname[i] = NULL;
 736             memset(&atoms->atom[i], 0, sizeof(atoms->atom[i]));
 737             if (NULL != atoms->pdbinfo)
 738             {
 739                 memset(&atoms->pdbinfo[i], 0, sizeof(atoms->pdbinfo[i]));
 740             }
 741             if (NULL != atoms->atomtype)
 742             {
 743                 atoms->atomtype[i] = NULL;
 744             }
 745             if (NULL != atoms->atomtypeB)
 746             {
 747                 atoms->atomtypeB[i] = NULL;
 748             }
 749         }
 750         atoms->nr += natom_extra;
 751     }
 752     if (nres_extra > 0)
 753     {
 754         srenew(atoms->resinfo, atoms->nres+nres_extra);
 755         for (i = atoms->nres; (i < atoms->nres+nres_extra); i++)
 756         {
 757             memset(&atoms->resinfo[i], 0, sizeof(atoms->resinfo[i]));
 758         }
 759         atoms->nres += nres_extra;
 760     }
 761 }
 762
 763 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
 764 {
 765     atoms->nr   = natoms;
 766     atoms->nres = 0;
 767     snew(atoms->atomname, natoms);
 768     atoms->atomtype  = NULL;
 769     atoms->atomtypeB = NULL;
 770     snew(atoms->resinfo, natoms);
 771     snew(atoms->atom, natoms);
 772     if (bPdbinfo)
 773     {
 774         snew(atoms->pdbinfo, natoms);
 775     }
 776     else
 777     {
 778         atoms->pdbinfo = NULL;
 779     }
 780 }
 781
 782 t_atoms *copy_t_atoms(t_atoms *src)
 783 {
 784     t_atoms *dst;
 785     int      i;
 786
 787     snew(dst, 1);
 788     init_t_atoms(dst, src->nr, (NULL != src->pdbinfo));
 789     dst->nr = src->nr;
 790     if (NULL != src->atomname)
 791     {
 792         snew(dst->atomname, src->nr);
 793     }
 794     if (NULL != src->atomtype)
 795     {
 796         snew(dst->atomtype, src->nr);
 797     }
 798     if (NULL != src->atomtypeB)
 799     {
 800         snew(dst->atomtypeB, src->nr);
 801     }
 802     for (i = 0; (i < src->nr); i++)
 803     {
 804         dst->atom[i] = src->atom[i];
 805         if (NULL != src->pdbinfo)
 806         {
 807             dst->pdbinfo[i] = src->pdbinfo[i];
 808         }
 809         if (NULL != src->atomname)
 810         {
 811             dst->atomname[i]  = src->atomname[i];
 812         }
 813         if (NULL != src->atomtype)
 814         {
 815             dst->atomtype[i] = src->atomtype[i];
 816         }
 817         if (NULL != src->atomtypeB)
 818         {
 819             dst->atomtypeB[i] = src->atomtypeB[i];
 820         }
 821     }
 822     dst->nres = src->nres;
 823     for (i = 0; (i < src->nres); i++)
 824     {
 825         dst->resinfo[i] = src->resinfo[i];
 826     }
 827     return dst;
 828 }
 829
 830 void t_atoms_set_resinfo(t_atoms *atoms, int atom_ind, t_symtab *symtab,
 831                          const char *resname, int resnr, unsigned char ic,
 832                          int chainnum, char chainid)
 833 {
 834     t_resinfo *ri;
 835
 836     ri           = &atoms->resinfo[atoms->atom[atom_ind].resind];
 837     ri->name     = put_symtab(symtab, resname);
 838     ri->rtp      = NULL;
 839     ri->nr       = resnr;
 840     ri->ic       = ic;
 841     ri->chainnum = chainnum;
 842     ri->chainid  = chainid;
 843 }
 844
 845 void free_t_atoms(t_atoms *atoms, gmx_bool bFreeNames)
 846 {
 847     int i;
 848
 849     if (bFreeNames && atoms->atomname != NULL)
 850     {
 851         for (i = 0; i < atoms->nr; i++)
 852         {
 853             if (atoms->atomname[i] != NULL)
 854             {
 855                 sfree(*atoms->atomname[i]);
 856                 *atoms->atomname[i] = NULL;
 857             }
 858         }
 859     }
 860     if (bFreeNames && atoms->resinfo != NULL)
 861     {
 862         for (i = 0; i < atoms->nres; i++)
 863         {
 864             if (atoms->resinfo[i].name != NULL)
 865             {
 866                 sfree(*atoms->resinfo[i].name);
 867                 *atoms->resinfo[i].name = NULL;
 868             }
 869         }
 870     }
 871     if (bFreeNames && atoms->atomtype != NULL)
 872     {
 873         for (i = 0; i < atoms->nr; i++)
 874         {
 875             if (atoms->atomtype[i] != NULL)
 876             {
 877                 sfree(*atoms->atomtype[i]);
 878                 *atoms->atomtype[i] = NULL;
 879             }
 880         }
 881     }
 882     if (bFreeNames && atoms->atomtypeB != NULL)
 883     {
 884         for (i = 0; i < atoms->nr; i++)
 885         {
 886             if (atoms->atomtypeB[i] != NULL)
 887             {
 888                 sfree(*atoms->atomtypeB[i]);
 889                 *atoms->atomtypeB[i] = NULL;
 890             }
 891         }
 892     }
 893     sfree(atoms->atomname);
 894     sfree(atoms->atomtype);
 895     sfree(atoms->atomtypeB);
 896     sfree(atoms->resinfo);
 897     sfree(atoms->atom);
 898     sfree(atoms->pdbinfo);
 899     atoms->nr        = 0;
 900     atoms->nres      = 0;
 901     atoms->atomname  = NULL;
 902     atoms->atomtype  = NULL;
 903     atoms->atomtypeB = NULL;
 904     atoms->resinfo   = NULL;
 905     atoms->atom      = NULL;
 906     atoms->pdbinfo   = NULL;
 907 }
 908
 909 real max_cutoff(real cutoff1, real cutoff2)
 910 {
 911     if (cutoff1 == 0 || cutoff2 == 0)
 912     {
 913         return 0;
 914     }
 915     else
 916     {
 917         return max(cutoff1, cutoff2);
 918     }
 919 }
 920
 921 void init_df_history(df_history_t *dfhist, int nlambda)
 922 {
 923     int i;
 924
 925     dfhist->nlambda  = nlambda;
 926     dfhist->bEquil   = 0;
 927     dfhist->wl_delta = 0;
 928
 929     if (nlambda > 0)
 930     {
 931         snew(dfhist->sum_weights, dfhist->nlambda);
 932         snew(dfhist->sum_dg, dfhist->nlambda);
 933         snew(dfhist->sum_minvar, dfhist->nlambda);
 934         snew(dfhist->sum_variance, dfhist->nlambda);
 935         snew(dfhist->n_at_lam, dfhist->nlambda);
 936         snew(dfhist->wl_histo, dfhist->nlambda);
 937
 938         /* allocate transition matrices here */
 939         snew(dfhist->Tij, dfhist->nlambda);
 940         snew(dfhist->Tij_empirical, dfhist->nlambda);
 941
 942         /* allocate accumulators for various transition matrix
 943            free energy methods here */
 944         snew(dfhist->accum_p, dfhist->nlambda);
 945         snew(dfhist->accum_m, dfhist->nlambda);
 946         snew(dfhist->accum_p2, dfhist->nlambda);
 947         snew(dfhist->accum_m2, dfhist->nlambda);
 948
 949         for (i = 0; i < dfhist->nlambda; i++)
 950         {
 951             snew(dfhist->Tij[i], dfhist->nlambda);
 952             snew(dfhist->Tij_empirical[i], dfhist->nlambda);
 953             snew((dfhist->accum_p)[i], dfhist->nlambda);
 954             snew((dfhist->accum_m)[i], dfhist->nlambda);
 955             snew((dfhist->accum_p2)[i], dfhist->nlambda);
 956             snew((dfhist->accum_m2)[i], dfhist->nlambda);
 957         }
 958     }
 959 }
 960
 961 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
 962 {
 963     int i, j;
 964
 965     /* Currently, there should not be any difference in nlambda between the two,
 966        but this is included for completeness for potential later functionality */
 967     df_dest->nlambda = df_source->nlambda;
 968     df_dest->bEquil  = df_source->bEquil;
 969     df_dest->wl_delta = df_source->wl_delta;
 970
 971     for (i = 0; i < df_dest->nlambda; i++)
 972     {
 973         df_dest->sum_weights[i]  = df_source->sum_weights[i];
 974         df_dest->sum_dg[i]       = df_source->sum_dg[i];
 975         df_dest->sum_minvar[i]   = df_source->sum_minvar[i];
 976         df_dest->sum_variance[i] = df_source->sum_variance[i];
 977         df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
 978         df_dest->wl_histo[i]     = df_source->wl_histo[i];
 979     }
 980
 981     for (i = 0; i < df_dest->nlambda; i++)
 982     {
 983         for (j = 0; j < df_dest->nlambda; j++)
 984         {
 985             df_dest->accum_p[i][j]      = df_source->accum_p[i][j];
 986             df_dest->accum_m[i][j]      = df_source->accum_m[i][j];
 987             df_dest->accum_p2[i][j]     = df_source->accum_p2[i][j];
 988             df_dest->accum_m2[i][j]     = df_source->accum_m2[i][j];
 989             df_dest->Tij[i][j]            = df_source->Tij[i][j];
 990             df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
 991         }
 992     }
 993 }
 994
 995 void done_df_history(df_history_t *dfhist)
 996 {
 997     int i;
 998
 999     if (dfhist->nlambda > 0)
1000     {
1001         sfree(dfhist->n_at_lam);
1002         sfree(dfhist->wl_histo);
1003         sfree(dfhist->sum_weights);
1004         sfree(dfhist->sum_dg);
1005         sfree(dfhist->sum_minvar);
1006         sfree(dfhist->sum_variance);
1007
1008         for (i=0;i<dfhist->nlambda;i++)
1009         {
1010             sfree(dfhist->Tij[i]);
1011             sfree(dfhist->Tij_empirical[i]);
1012             sfree(dfhist->accum_p[i]);
1013             sfree(dfhist->accum_m[i]);
1014             sfree(dfhist->accum_p2[i]);
1015             sfree(dfhist->accum_m2[i]);
1016         }
1017     }
1018     dfhist->bEquil   = 0;
1019     dfhist->nlambda  = 0;
1020     dfhist->wl_delta = 0;
1021 }