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     return OVER_ALLOC_FAC*n + 100;
  80   else
  81     return n;
  82 }
  83
  84 int gmx_large_int_to_int(gmx_large_int_t step,const char *warn)
  85 {
  86   int i;
  87
  88   i = (int)step;
  89
  90   if (warn != NULL && (step < INT_MIN || step > INT_MAX)) {
  91     fprintf(stderr,"\nWARNING during %s:\n",warn);
  92     fprintf(stderr,"step value ");
  93     fprintf(stderr,gmx_large_int_pfmt,step);
  94     fprintf(stderr," does not fit in int, converted to %d\n\n",i);
  95   }
  96
  97   return i;
  98 }
  99
 100 char *gmx_step_str(gmx_large_int_t i,char *buf)
 101 {
 102   sprintf(buf,gmx_large_int_pfmt,i);
 103
 104   return buf;
 105 }
 106
 107 void init_block(t_block *block)
 108 {
 109   int i;
 110
 111   block->nr           = 0;
 112   block->nalloc_index = 1;
 113   snew(block->index,block->nalloc_index);
 114   block->index[0]     = 0;
 115 }
 116
 117 void init_blocka(t_blocka *block)
 118 {
 119   int i;
 120
 121   block->nr           = 0;
 122   block->nra          = 0;
 123   block->nalloc_index = 1;
 124   snew(block->index,block->nalloc_index);
 125   block->index[0]     = 0;
 126   block->nalloc_a     = 0;
 127   block->a            = NULL;
 128 }
 129
 130 void init_atom(t_atoms *at)
 131 {
 132   int i;
 133
 134   at->nr       = 0;
 135   at->nres     = 0;
 136   at->atom     = NULL;
 137   at->resinfo  = NULL;
 138   at->atomname = NULL;
 139   at->atomtype = NULL;
 140   at->atomtypeB= NULL;
 141   at->pdbinfo  = NULL;
 142 }
 143
 144 void init_atomtypes(t_atomtypes *at)
 145 {
 146   at->nr = 0;
 147   at->radius = NULL;
 148   at->vol = NULL;
 149   at->atomnumber = NULL;
 150   at->gb_radius = NULL;
 151   at->S_hct = NULL;
 152 }
 153
 154 void init_groups(gmx_groups_t *groups)
 155 {
 156   int g;
 157
 158   groups->ngrpname = 0;
 159   groups->grpname  = NULL;
 160   for(g=0; (g<egcNR); g++) {
 161     groups->grps[g].nm_ind = NULL;
 162     groups->ngrpnr[g] = 0;
 163     groups->grpnr[g]  = NULL;
 164   }
 165
 166 }
 167
 168 void init_mtop(gmx_mtop_t *mtop)
 169 {
 170   mtop->name = NULL;
 171   mtop->nmoltype = 0;
 172   mtop->moltype = NULL;
 173   mtop->nmolblock = 0;
 174   mtop->molblock = NULL;
 175   mtop->maxres_renum = 0;
 176   mtop->maxresnr = -1;
 177   init_groups(&mtop->groups);
 178   init_block(&mtop->mols);
 179   open_symtab(&mtop->symtab);
 180 }
 181
 182 void init_top (t_topology *top)
 183 {
 184   int i;
 185
 186   top->name = NULL;
 187   init_atom (&(top->atoms));
 188   init_atomtypes(&(top->atomtypes));
 189   init_block(&top->cgs);
 190   init_block(&top->mols);
 191   init_blocka(&top->excls);
 192   open_symtab(&top->symtab);
 193 }
 194
 195 void init_inputrec(t_inputrec *ir)
 196 {
 197     memset(ir,0,(size_t)sizeof(*ir));
 198     snew(ir->fepvals,1);
 199     snew(ir->expandedvals,1);
 200     snew(ir->simtempvals,1);
 201 }
 202
 203 void stupid_fill_block(t_block *grp,int natom,gmx_bool bOneIndexGroup)
 204 {
 205   int i;
 206
 207   if (bOneIndexGroup) {
 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     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       grp->index[i]=i;
 220     grp->nr=natom;
 221   }
 222 }
 223
 224 void stupid_fill_blocka(t_blocka *grp,int natom)
 225 {
 226   int i;
 227
 228   grp->nalloc_a = natom;
 229   snew(grp->a,grp->nalloc_a);
 230   for(i=0; (i<natom); i++)
 231     grp->a[i]=i;
 232   grp->nra=natom;
 233
 234   grp->nalloc_index = natom + 1;
 235   snew(grp->index,grp->nalloc_index);
 236   for(i=0; (i<=natom); i++)
 237     grp->index[i]=i;
 238   grp->nr=natom;
 239 }
 240
 241 void copy_blocka(const t_blocka *src,t_blocka *dest)
 242 {
 243   int i;
 244
 245   dest->nr = src->nr;
 246   dest->nalloc_index = dest->nr + 1;
 247   snew(dest->index,dest->nalloc_index);
 248   for(i=0; i<dest->nr+1; i++) {
 249     dest->index[i] = src->index[i];
 250   }
 251   dest->nra = src->nra;
 252   dest->nalloc_a = dest->nra + 1;
 253   snew(dest->a,dest->nalloc_a);
 254   for(i=0; i<dest->nra+1; i++) {
 255     dest->a[i] = src->a[i];
 256   }
 257 }
 258
 259 void done_block(t_block *block)
 260 {
 261   block->nr    = 0;
 262   sfree(block->index);
 263   block->nalloc_index = 0;
 264 }
 265
 266 void done_blocka(t_blocka *block)
 267 {
 268   block->nr    = 0;
 269   block->nra   = 0;
 270   sfree(block->index);
 271   if (block->a)
 272     sfree(block->a);
 273   block->nalloc_index = 0;
 274   block->nalloc_a = 0;
 275 }
 276
 277 void done_atom (t_atoms *at)
 278 {
 279   at->nr       = 0;
 280   at->nres     = 0;
 281   sfree(at->atom);
 282   sfree(at->resinfo);
 283   sfree(at->atomname);
 284   sfree(at->atomtype);
 285   sfree(at->atomtypeB);
 286   if (at->pdbinfo)
 287     sfree(at->pdbinfo);
 288 }
 289
 290 void done_atomtypes(t_atomtypes *atype)
 291 {
 292   atype->nr = 0;
 293   sfree(atype->radius);
 294   sfree(atype->vol);
 295   sfree(atype->surftens);
 296   sfree(atype->atomnumber);
 297   sfree(atype->gb_radius);
 298   sfree(atype->S_hct);
 299 }
 300
 301 void done_moltype(gmx_moltype_t *molt)
 302 {
 303   int f;
 304
 305   done_atom(&molt->atoms);
 306   done_block(&molt->cgs);
 307   done_blocka(&molt->excls);
 308
 309   for(f=0; f<F_NRE; f++) {
 310     sfree(molt->ilist[f].iatoms);
 311     molt->ilist[f].nalloc = 0;
 312   }
 313 }
 314
 315 void done_molblock(gmx_molblock_t *molb)
 316 {
 317   if (molb->nposres_xA > 0) {
 318     molb->nposres_xA = 0;
 319     free(molb->posres_xA);
 320   }
 321   if (molb->nposres_xB > 0) {
 322     molb->nposres_xB = 0;
 323     free(molb->posres_xB);
 324   }
 325 }
 326
 327 void done_mtop(gmx_mtop_t *mtop,gmx_bool bDoneSymtab)
 328 {
 329   int i;
 330
 331   if (bDoneSymtab) {
 332     done_symtab(&mtop->symtab);
 333   }
 334
 335   sfree(mtop->ffparams.functype);
 336   sfree(mtop->ffparams.iparams);
 337
 338   for(i=0; i<mtop->nmoltype; i++) {
 339     done_moltype(&mtop->moltype[i]);
 340   }
 341   sfree(mtop->moltype);
 342   for(i=0; i<mtop->nmolblock; i++) {
 343     done_molblock(&mtop->molblock[i]);
 344   }
 345   sfree(mtop->molblock);
 346   done_block(&mtop->mols);
 347 }
 348
 349 void done_top(t_topology *top)
 350 {
 351   int f;
 352
 353   sfree(top->idef.functype);
 354   sfree(top->idef.iparams);
 355   for (f = 0; f < F_NRE; ++f)
 356   {
 357       sfree(top->idef.il[f].iatoms);
 358       top->idef.il[f].iatoms = NULL;
 359       top->idef.il[f].nalloc = 0;
 360   }
 361
 362   done_atom (&(top->atoms));
 363
 364   /* For GB */
 365   done_atomtypes(&(top->atomtypes));
 366
 367   done_symtab(&(top->symtab));
 368   done_block(&(top->cgs));
 369   done_block(&(top->mols));
 370   done_blocka(&(top->excls));
 371 }
 372
 373 static void done_pullgrp(t_pullgrp *pgrp)
 374 {
 375   sfree(pgrp->ind);
 376   sfree(pgrp->ind_loc);
 377   sfree(pgrp->weight);
 378   sfree(pgrp->weight_loc);
 379 }
 380
 381 static void done_pull(t_pull *pull)
 382 {
 383   int i;
 384
 385   for(i=0; i<pull->ngrp+1; i++) {
 386     done_pullgrp(pull->grp);
 387     done_pullgrp(pull->dyna);
 388   }
 389 }
 390
 391 void done_inputrec(t_inputrec *ir)
 392 {
 393   int m;
 394
 395   for(m=0; (m<DIM); m++) {
 396     if (ir->ex[m].a)   sfree(ir->ex[m].a);
 397     if (ir->ex[m].phi) sfree(ir->ex[m].phi);
 398     if (ir->et[m].a)   sfree(ir->et[m].a);
 399     if (ir->et[m].phi) sfree(ir->et[m].phi);
 400   }
 401
 402   sfree(ir->opts.nrdf);
 403   sfree(ir->opts.ref_t);
 404   sfree(ir->opts.annealing);
 405   sfree(ir->opts.anneal_npoints);
 406   sfree(ir->opts.anneal_time);
 407   sfree(ir->opts.anneal_temp);
 408   sfree(ir->opts.tau_t);
 409   sfree(ir->opts.acc);
 410   sfree(ir->opts.nFreeze);
 411   sfree(ir->opts.QMmethod);
 412   sfree(ir->opts.QMbasis);
 413   sfree(ir->opts.QMcharge);
 414   sfree(ir->opts.QMmult);
 415   sfree(ir->opts.bSH);
 416   sfree(ir->opts.CASorbitals);
 417   sfree(ir->opts.CASelectrons);
 418   sfree(ir->opts.SAon);
 419   sfree(ir->opts.SAoff);
 420   sfree(ir->opts.SAsteps);
 421   sfree(ir->opts.bOPT);
 422   sfree(ir->opts.bTS);
 423
 424   if (ir->pull) {
 425     done_pull(ir->pull);
 426     sfree(ir->pull);
 427   }
 428 }
 429
 430 static void zero_ekinstate(ekinstate_t *eks)
 431 {
 432   eks->ekin_n         = 0;
 433   eks->ekinh          = NULL;
 434   eks->ekinf          = NULL;
 435   eks->ekinh_old      = NULL;
 436   eks->ekinscalef_nhc = NULL;
 437   eks->ekinscaleh_nhc = NULL;
 438   eks->vscale_nhc     = NULL;
 439   eks->dekindl        = 0;
 440   eks->mvcos          = 0;
 441 }
 442
 443 void init_energyhistory(energyhistory_t * enerhist)
 444 {
 445     enerhist->nener = 0;
 446
 447     enerhist->ener_ave     = NULL;
 448     enerhist->ener_sum     = NULL;
 449     enerhist->ener_sum_sim = NULL;
 450     enerhist->dht          = NULL;
 451
 452     enerhist->nsteps     = 0;
 453     enerhist->nsum       = 0;
 454     enerhist->nsteps_sim = 0;
 455     enerhist->nsum_sim   = 0;
 456
 457     enerhist->dht = NULL;
 458 }
 459
 460 static void done_delta_h_history(delta_h_history_t *dht)
 461 {
 462     int i;
 463
 464     for(i=0; i<dht->nndh; i++)
 465     {
 466         sfree(dht->dh[i]);
 467     }
 468     sfree(dht->dh);
 469     sfree(dht->ndh);
 470 }
 471
 472 void done_energyhistory(energyhistory_t * enerhist)
 473 {
 474     sfree(enerhist->ener_ave);
 475     sfree(enerhist->ener_sum);
 476     sfree(enerhist->ener_sum_sim);
 477
 478     if (enerhist->dht != NULL)
 479     {
 480         done_delta_h_history(enerhist->dht);
 481         sfree(enerhist->dht);
 482     }
 483 }
 484
 485 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
 486 {
 487     int i,j;
 488
 489     state->ngtc = ngtc;
 490     state->nnhpres = nnhpres;
 491     state->nhchainlength = nhchainlength;
 492     if (state->ngtc > 0)
 493     {
 494         snew(state->nosehoover_xi,state->nhchainlength*state->ngtc);
 495         snew(state->nosehoover_vxi,state->nhchainlength*state->ngtc);
 496         snew(state->therm_integral,state->ngtc);
 497         for(i=0; i<state->ngtc; i++)
 498         {
 499             for (j=0;j<state->nhchainlength;j++)
 500             {
 501                 state->nosehoover_xi[i*state->nhchainlength + j]  = 0.0;
 502                 state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
 503             }
 504         }
 505         for(i=0; i<state->ngtc; i++) {
 506             state->therm_integral[i]  = 0.0;
 507         }
 508     }
 509     else
 510     {
 511         state->nosehoover_xi  = NULL;
 512         state->nosehoover_vxi = NULL;
 513         state->therm_integral = NULL;
 514     }
 515
 516     if (state->nnhpres > 0)
 517     {
 518         snew(state->nhpres_xi,state->nhchainlength*nnhpres);
 519         snew(state->nhpres_vxi,state->nhchainlength*nnhpres);
 520         for(i=0; i<nnhpres; i++)
 521         {
 522             for (j=0;j<state->nhchainlength;j++)
 523             {
 524                 state->nhpres_xi[i*nhchainlength + j]  = 0.0;
 525                 state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
 526             }
 527         }
 528     }
 529     else
 530     {
 531         state->nhpres_xi  = NULL;
 532         state->nhpres_vxi = NULL;
 533     }
 534 }
 535
 536
 537 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
 538 {
 539   int i;
 540
 541   state->natoms = natoms;
 542   state->nrng   = 0;
 543   state->flags  = 0;
 544   state->lambda = 0;
 545   snew(state->lambda,efptNR);
 546   for (i=0;i<efptNR;i++)
 547   {
 548       state->lambda[i] = 0;
 549   }
 550   state->veta   = 0;
 551   clear_mat(state->box);
 552   clear_mat(state->box_rel);
 553   clear_mat(state->boxv);
 554   clear_mat(state->pres_prev);
 555   clear_mat(state->svir_prev);
 556   clear_mat(state->fvir_prev);
 557   init_gtc_state(state,ngtc,nnhpres,nhchainlength);
 558   state->nalloc = state->natoms;
 559   if (state->nalloc > 0) {
 560     snew(state->x,state->nalloc);
 561     snew(state->v,state->nalloc);
 562   } else {
 563     state->x = NULL;
 564     state->v = NULL;
 565   }
 566   state->sd_X = NULL;
 567   state->cg_p = NULL;
 568
 569   zero_ekinstate(&state->ekinstate);
 570
 571   init_energyhistory(&state->enerhist);
 572
 573   init_df_history(&state->dfhist,nlambda,0);
 574
 575   state->ddp_count = 0;
 576   state->ddp_count_cg_gl = 0;
 577   state->cg_gl = NULL;
 578   state->cg_gl_nalloc = 0;
 579 }
 580
 581 void done_state(t_state *state)
 582 {
 583   if (state->nosehoover_xi) sfree(state->nosehoover_xi);
 584   if (state->x) sfree(state->x);
 585   if (state->v) sfree(state->v);
 586   if (state->sd_X) sfree(state->sd_X);
 587   if (state->cg_p) sfree(state->cg_p);
 588   state->nalloc = 0;
 589   if (state->cg_gl) sfree(state->cg_gl);
 590   state->cg_gl_nalloc = 0;
 591 }
 592
 593 static void do_box_rel(t_inputrec *ir,matrix box_rel,matrix b,gmx_bool bInit)
 594 {
 595   int d,d2;
 596
 597   for(d=YY; d<=ZZ; d++) {
 598     for(d2=XX; d2<=(ir->epct==epctSEMIISOTROPIC ? YY : ZZ); d2++) {
 599       /* We need to check if this box component is deformed
 600        * or if deformation of another component might cause
 601        * changes in this component due to box corrections.
 602        */
 603       if (ir->deform[d][d2] == 0 &&
 604           !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
 605             (b[YY][d2] != 0 || ir->deform[YY][d2] != 0))) {
 606         if (bInit) {
 607           box_rel[d][d2] = b[d][d2]/b[XX][XX];
 608         } else {
 609           b[d][d2] = b[XX][XX]*box_rel[d][d2];
 610         }
 611       }
 612     }
 613   }
 614 }
 615
 616 void set_box_rel(t_inputrec *ir,t_state *state)
 617 {
 618   /* Make sure the box obeys the restrictions before we fix the ratios */
 619   correct_box(NULL,0,state->box,NULL);
 620
 621   clear_mat(state->box_rel);
 622
 623   if (PRESERVE_SHAPE(*ir))
 624     do_box_rel(ir,state->box_rel,state->box,TRUE);
 625 }
 626
 627 void preserve_box_shape(t_inputrec *ir,matrix box_rel,matrix b)
 628 {
 629   if (PRESERVE_SHAPE(*ir))
 630     do_box_rel(ir,box_rel,b,FALSE);
 631 }
 632
 633 void add_t_atoms(t_atoms *atoms,int natom_extra,int nres_extra)
 634 {
 635     int i;
 636
 637     if (natom_extra > 0)
 638     {
 639         srenew(atoms->atomname,atoms->nr+natom_extra);
 640         srenew(atoms->atom,atoms->nr+natom_extra);
 641         if (NULL != atoms->pdbinfo)
 642             srenew(atoms->pdbinfo,atoms->nr+natom_extra);
 643         if (NULL != atoms->atomtype)
 644             srenew(atoms->atomtype,atoms->nr+natom_extra);
 645         if (NULL != atoms->atomtypeB)
 646             srenew(atoms->atomtypeB,atoms->nr+natom_extra);
 647         for(i=atoms->nr; (i<atoms->nr+natom_extra); i++) {
 648             atoms->atomname[i] = NULL;
 649             memset(&atoms->atom[i],0,sizeof(atoms->atom[i]));
 650             if (NULL != atoms->pdbinfo)
 651                 memset(&atoms->pdbinfo[i],0,sizeof(atoms->pdbinfo[i]));
 652             if (NULL != atoms->atomtype)
 653                 atoms->atomtype[i] = NULL;
 654             if (NULL != atoms->atomtypeB)
 655                 atoms->atomtypeB[i] = NULL;
 656         }
 657         atoms->nr += natom_extra;
 658     }
 659     if (nres_extra > 0)
 660     {
 661         srenew(atoms->resinfo,atoms->nres+nres_extra);
 662         for(i=atoms->nres; (i<atoms->nres+nres_extra); i++) {
 663             memset(&atoms->resinfo[i],0,sizeof(atoms->resinfo[i]));
 664         }
 665         atoms->nres += nres_extra;
 666     }
 667 }
 668
 669 void init_t_atoms(t_atoms *atoms, int natoms, gmx_bool bPdbinfo)
 670 {
 671   atoms->nr=natoms;
 672   atoms->nres=0;
 673   snew(atoms->atomname,natoms);
 674   atoms->atomtype=NULL;
 675   atoms->atomtypeB=NULL;
 676   snew(atoms->resinfo,natoms);
 677   snew(atoms->atom,natoms);
 678   if (bPdbinfo)
 679     snew(atoms->pdbinfo,natoms);
 680   else
 681     atoms->pdbinfo=NULL;
 682 }
 683
 684 t_atoms *copy_t_atoms(t_atoms *src)
 685 {
 686   t_atoms *dst;
 687   int i;
 688
 689   snew(dst,1);
 690   init_t_atoms(dst,src->nr,(NULL != src->pdbinfo));
 691   dst->nr = src->nr;
 692   if (NULL != src->atomname)
 693       snew(dst->atomname,src->nr);
 694   if (NULL != src->atomtype)
 695       snew(dst->atomtype,src->nr);
 696   if (NULL != src->atomtypeB)
 697       snew(dst->atomtypeB,src->nr);
 698   for(i=0; (i<src->nr); i++) {
 699     dst->atom[i] = src->atom[i];
 700     if (NULL != src->pdbinfo)
 701       dst->pdbinfo[i] = src->pdbinfo[i];
 702     if (NULL != src->atomname)
 703         dst->atomname[i]  = src->atomname[i];
 704     if (NULL != src->atomtype)
 705         dst->atomtype[i] = src->atomtype[i];
 706     if (NULL != src->atomtypeB)
 707         dst->atomtypeB[i] = src->atomtypeB[i];
 708   }
 709   dst->nres = src->nres;
 710   for(i=0; (i<src->nres); i++) {
 711     dst->resinfo[i] = src->resinfo[i];
 712   }
 713   return dst;
 714 }
 715
 716 void t_atoms_set_resinfo(t_atoms *atoms,int atom_ind,t_symtab *symtab,
 717                          const char *resname,int resnr,unsigned char ic,
 718                          int chainnum, char chainid)
 719 {
 720   t_resinfo *ri;
 721
 722   ri = &atoms->resinfo[atoms->atom[atom_ind].resind];
 723   ri->name  = put_symtab(symtab,resname);
 724   ri->rtp   = NULL;
 725   ri->nr    = resnr;
 726   ri->ic    = ic;
 727   ri->chainnum = chainnum;
 728   ri->chainid = chainid;
 729 }
 730
 731 void free_t_atoms(t_atoms *atoms,gmx_bool bFreeNames)
 732 {
 733   int i;
 734
 735   if (bFreeNames) {
 736     for(i=0; i<atoms->nr; i++) {
 737       sfree(*atoms->atomname[i]);
 738       *atoms->atomname[i]=NULL;
 739     }
 740     for(i=0; i<atoms->nres; i++) {
 741       sfree(*atoms->resinfo[i].name);
 742       *atoms->resinfo[i].name=NULL;
 743     }
 744   }
 745   sfree(atoms->atomname);
 746   /* Do we need to free atomtype and atomtypeB as well ? */
 747   sfree(atoms->resinfo);
 748   sfree(atoms->atom);
 749   if (atoms->pdbinfo)
 750     sfree(atoms->pdbinfo);
 751   atoms->nr=0;
 752   atoms->nres=0;
 753   atoms->atomname = NULL;
 754   atoms->resinfo = NULL;
 755   atoms->atom = NULL;
 756   atoms->pdbinfo = NULL;
 757 }
 758
 759 real max_cutoff(real cutoff1,real cutoff2)
 760 {
 761     if (cutoff1 == 0 || cutoff2 == 0)
 762     {
 763         return 0;
 764     }
 765     else
 766     {
 767         return max(cutoff1,cutoff2);
 768     }
 769 }
 770
 771 extern void init_df_history(df_history_t *dfhist, int nlambda, real wl_delta)
 772 {
 773     int i;
 774
 775     dfhist->bEquil = 0;
 776     dfhist->nlambda = nlambda;
 777     dfhist->wl_delta = wl_delta;
 778     snew(dfhist->sum_weights,dfhist->nlambda);
 779     snew(dfhist->sum_dg,dfhist->nlambda);
 780     snew(dfhist->sum_minvar,dfhist->nlambda);
 781     snew(dfhist->sum_variance,dfhist->nlambda);
 782     snew(dfhist->n_at_lam,dfhist->nlambda);
 783     snew(dfhist->wl_histo,dfhist->nlambda);
 784
 785     /* allocate transition matrices here */
 786     snew(dfhist->Tij,dfhist->nlambda);
 787     snew(dfhist->Tij_empirical,dfhist->nlambda);
 788
 789     for (i=0;i<dfhist->nlambda;i++) {
 790         snew(dfhist->Tij[i],dfhist->nlambda);
 791         snew(dfhist->Tij_empirical[i],dfhist->nlambda);
 792     }
 793
 794     snew(dfhist->accum_p,dfhist->nlambda);
 795     snew(dfhist->accum_m,dfhist->nlambda);
 796     snew(dfhist->accum_p2,dfhist->nlambda);
 797     snew(dfhist->accum_m2,dfhist->nlambda);
 798
 799     for (i=0;i<dfhist->nlambda;i++) {
 800         snew((dfhist->accum_p)[i],dfhist->nlambda);
 801         snew((dfhist->accum_m)[i],dfhist->nlambda);
 802         snew((dfhist->accum_p2)[i],dfhist->nlambda);
 803         snew((dfhist->accum_m2)[i],dfhist->nlambda);
 804     }
 805 }
 806
 807 extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
 808 {
 809     int i,j;
 810
 811     init_df_history(df_dest,df_source->nlambda,df_source->wl_delta);
 812     df_dest->nlambda = df_source->nlambda;
 813     df_dest->bEquil = df_source->bEquil;
 814     for (i=0;i<df_dest->nlambda;i++)
 815     {
 816         df_dest->sum_weights[i]  = df_source->sum_weights[i];
 817         df_dest->sum_dg[i]       = df_source->sum_dg[i];
 818         df_dest->sum_minvar[i]   = df_source->sum_minvar[i];
 819         df_dest->sum_variance[i] = df_source->sum_variance[i];
 820         df_dest->n_at_lam[i]     = df_source->n_at_lam[i];
 821         df_dest->wl_histo[i]     = df_source->wl_histo[i];
 822         df_dest->accum_p[i]      = df_source->accum_p[i];
 823         df_dest->accum_m[i]      = df_source->accum_m[i];
 824         df_dest->accum_p2[i]     = df_source->accum_p2[i];
 825         df_dest->accum_m2[i]     = df_source->accum_m2[i];
 826     }
 827
 828     for (i=0;i<df_dest->nlambda;i++)
 829     {
 830         for (j=0;j<df_dest->nlambda;j++)
 831         {
 832             df_dest->Tij[i][j]  = df_source->Tij[i][j];
 833             df_dest->Tij_empirical[i][j]  = df_source->Tij_empirical[i][j];
 834         }
 835     }
 836 }