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