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_THREADS
  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 bool bOverAllocDD=FALSE;
  57 #ifdef GMX_THREADS
  58 static tMPI_Thread_mutex_t over_alloc_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
  59 #endif
  60
  61
  62 void set_over_alloc_dd(bool set)
  63 {
  64 #ifdef GMX_THREADS
  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_THREADS
  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 }
 198
 199 void stupid_fill_block(t_block *grp,int natom,bool bOneIndexGroup)
 200 {
 201   int i;
 202
 203   if (bOneIndexGroup) {
 204     grp->nalloc_index = 2;
 205     snew(grp->index,grp->nalloc_index);
 206     grp->index[0]=0;
 207     grp->index[1]=natom;
 208     grp->nr=1;
 209   }
 210   else {
 211     grp->nalloc_index = natom+1;
 212     snew(grp->index,grp->nalloc_index);
 213     snew(grp->index,natom+1);
 214     for(i=0; (i<=natom); i++)
 215       grp->index[i]=i;
 216     grp->nr=natom;
 217   }
 218 }
 219
 220 void stupid_fill_blocka(t_blocka *grp,int natom)
 221 {
 222   int i;
 223
 224   grp->nalloc_a = natom;
 225   snew(grp->a,grp->nalloc_a);
 226   for(i=0; (i<natom); i++)
 227     grp->a[i]=i;
 228   grp->nra=natom;
 229
 230   grp->nalloc_index = natom + 1;
 231   snew(grp->index,grp->nalloc_index);
 232   for(i=0; (i<=natom); i++)
 233     grp->index[i]=i;
 234   grp->nr=natom;
 235 }
 236
 237 void copy_blocka(const t_blocka *src,t_blocka *dest)
 238 {
 239   int i;
 240
 241   dest->nr = src->nr;
 242   dest->nalloc_index = dest->nr + 1;
 243   snew(dest->index,dest->nalloc_index);
 244   for(i=0; i<dest->nr+1; i++) {
 245     dest->index[i] = src->index[i];
 246   }
 247   dest->nra = src->nra;
 248   dest->nalloc_a = dest->nra + 1;
 249   snew(dest->a,dest->nalloc_a);
 250   for(i=0; i<dest->nra+1; i++) {
 251     dest->a[i] = src->a[i];
 252   }
 253 }
 254
 255 void done_block(t_block *block)
 256 {
 257   block->nr    = 0;
 258   sfree(block->index);
 259   block->nalloc_index = 0;
 260 }
 261
 262 void done_blocka(t_blocka *block)
 263 {
 264   block->nr    = 0;
 265   block->nra   = 0;
 266   sfree(block->index);
 267   if (block->a)
 268     sfree(block->a);
 269   block->nalloc_index = 0;
 270   block->nalloc_a = 0;
 271 }
 272
 273 void done_atom (t_atoms *at)
 274 {
 275   at->nr       = 0;
 276   at->nres     = 0;
 277   sfree(at->atom);
 278   sfree(at->resinfo);
 279   sfree(at->atomname);
 280   sfree(at->atomtype);
 281   sfree(at->atomtypeB);
 282 }
 283
 284 void done_atomtypes(t_atomtypes *atype)
 285 {
 286   atype->nr = 0;
 287   sfree(atype->radius);
 288   sfree(atype->vol);
 289   sfree(atype->surftens);
 290   sfree(atype->atomnumber);
 291   sfree(atype->gb_radius);
 292   sfree(atype->S_hct);
 293 }
 294
 295 void done_moltype(gmx_moltype_t *molt)
 296 {
 297   int f;
 298
 299   done_atom(&molt->atoms);
 300   done_block(&molt->cgs);
 301   done_blocka(&molt->excls);
 302
 303   for(f=0; f<F_NRE; f++) {
 304     sfree(molt->ilist[f].iatoms);
 305     molt->ilist[f].nalloc = 0;
 306   }
 307 }
 308
 309 void done_molblock(gmx_molblock_t *molb)
 310 {
 311   if (molb->nposres_xA > 0) {
 312     molb->nposres_xA = 0;
 313     free(molb->posres_xA);
 314   }
 315   if (molb->nposres_xB > 0) {
 316     molb->nposres_xB = 0;
 317     free(molb->posres_xB);
 318   }
 319 }
 320
 321 void done_mtop(gmx_mtop_t *mtop,bool bDoneSymtab)
 322 {
 323   int i;
 324
 325   if (bDoneSymtab) {
 326     done_symtab(&mtop->symtab);
 327   }
 328
 329   sfree(mtop->ffparams.functype);
 330   sfree(mtop->ffparams.iparams);
 331
 332   for(i=0; i<mtop->nmoltype; i++) {
 333     done_moltype(&mtop->moltype[i]);
 334   }
 335   sfree(mtop->moltype);
 336   for(i=0; i<mtop->nmolblock; i++) {
 337     done_molblock(&mtop->molblock[i]);
 338   }
 339   sfree(mtop->molblock);
 340   done_block(&mtop->mols);
 341 }
 342
 343 void done_top(t_topology *top)
 344 {
 345   int f;
 346
 347   sfree(top->idef.functype);
 348   sfree(top->idef.iparams);
 349   for (f = 0; f < F_NRE; ++f)
 350   {
 351       sfree(top->idef.il[f].iatoms);
 352       top->idef.il[f].iatoms = NULL;
 353       top->idef.il[f].nalloc = 0;
 354   }
 355
 356   done_atom (&(top->atoms));
 357
 358   /* For GB */
 359   done_atomtypes(&(top->atomtypes));
 360
 361   done_symtab(&(top->symtab));
 362   done_block(&(top->cgs));
 363   done_block(&(top->mols));
 364   done_blocka(&(top->excls));
 365 }
 366
 367 static void done_pullgrp(t_pullgrp *pgrp)
 368 {
 369   sfree(pgrp->ind);
 370   sfree(pgrp->ind_loc);
 371   sfree(pgrp->weight);
 372   sfree(pgrp->weight_loc);
 373 }
 374
 375 static void done_pull(t_pull *pull)
 376 {
 377   int i;
 378
 379   for(i=0; i<pull->ngrp+1; i++) {
 380     done_pullgrp(pull->grp);
 381     done_pullgrp(pull->dyna);
 382   }
 383 }
 384
 385 void done_inputrec(t_inputrec *ir)
 386 {
 387   int m;
 388
 389   for(m=0; (m<DIM); m++) {
 390     if (ir->ex[m].a)   sfree(ir->ex[m].a);
 391     if (ir->ex[m].phi) sfree(ir->ex[m].phi);
 392     if (ir->et[m].a)   sfree(ir->et[m].a);
 393     if (ir->et[m].phi) sfree(ir->et[m].phi);
 394   }
 395
 396   sfree(ir->opts.nrdf);
 397   sfree(ir->opts.ref_t);
 398   sfree(ir->opts.annealing);
 399   sfree(ir->opts.anneal_npoints);
 400   sfree(ir->opts.anneal_time);
 401   sfree(ir->opts.anneal_temp);
 402   sfree(ir->opts.tau_t);
 403   sfree(ir->opts.acc);
 404   sfree(ir->opts.nFreeze);
 405   sfree(ir->opts.QMmethod);
 406   sfree(ir->opts.QMbasis);
 407   sfree(ir->opts.QMcharge);
 408   sfree(ir->opts.QMmult);
 409   sfree(ir->opts.bSH);
 410   sfree(ir->opts.CASorbitals);
 411   sfree(ir->opts.CASelectrons);
 412   sfree(ir->opts.SAon);
 413   sfree(ir->opts.SAoff);
 414   sfree(ir->opts.SAsteps);
 415   sfree(ir->opts.bOPT);
 416   sfree(ir->opts.bTS);
 417
 418   if (ir->pull) {
 419     done_pull(ir->pull);
 420     sfree(ir->pull);
 421   }
 422 }
 423
 424 static void init_ekinstate(ekinstate_t *eks)
 425 {
 426   eks->ekin_n         = 0;
 427   eks->ekinh          = NULL;
 428   eks->ekinf          = NULL;
 429   eks->ekinh_old      = NULL;
 430   eks->ekinscalef_nhc = NULL;
 431   eks->ekinscaleh_nhc = NULL;
 432   eks->vscale_nhc     = NULL;
 433   eks->dekindl        = 0;
 434   eks->mvcos          = 0;
 435 }
 436
 437 static void init_energyhistory(energyhistory_t *enh)
 438 {
 439   enh->ener_ave     = NULL;
 440   enh->ener_sum     = NULL;
 441   enh->ener_sum_sim = NULL;
 442   enh->nener        = 0;
 443 }
 444
 445 void init_gtc_state(t_state *state, int ngtc, int nnhpres, int nhchainlength)
 446 {
 447     int i,j;
 448
 449     state->ngtc = ngtc;
 450     state->nnhpres = nnhpres;
 451     state->nhchainlength = nhchainlength;
 452     if (state->ngtc > 0)
 453     {
 454         snew(state->nosehoover_xi,state->nhchainlength*state->ngtc);
 455         snew(state->nosehoover_vxi,state->nhchainlength*state->ngtc);
 456         snew(state->therm_integral,state->ngtc);
 457         for(i=0; i<state->ngtc; i++)
 458         {
 459             for (j=0;j<state->nhchainlength;j++)
 460  {
 461                 state->nosehoover_xi[i*state->nhchainlength + j]  = 0.0;
 462                 state->nosehoover_vxi[i*state->nhchainlength + j]  = 0.0;
 463             }
 464         }
 465         for(i=0; i<state->ngtc; i++) {
 466             state->therm_integral[i]  = 0.0;
 467         }
 468     }
 469     else
 470     {
 471         state->nosehoover_xi  = NULL;
 472         state->nosehoover_vxi = NULL;
 473         state->therm_integral = NULL;
 474     }
 475
 476     if (state->nnhpres > 0)
 477     {
 478         snew(state->nhpres_xi,state->nhchainlength*nnhpres);
 479         snew(state->nhpres_vxi,state->nhchainlength*nnhpres);
 480         for(i=0; i<nnhpres; i++)
 481         {
 482             for (j=0;j<state->nhchainlength;j++)
 483             {
 484                 state->nhpres_xi[i*nhchainlength + j]  = 0.0;
 485                 state->nhpres_vxi[i*nhchainlength + j]  = 0.0;
 486             }
 487         }
 488     }
 489     else
 490     {
 491         state->nhpres_xi  = NULL;
 492         state->nhpres_vxi = NULL;
 493     }
 494 }
 495
 496
 497 void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength)
 498 {
 499   int i;
 500
 501   state->natoms = natoms;
 502   state->nrng   = 0;
 503   state->flags  = 0;
 504   state->lambda = 0;
 505   state->veta   = 0;
 506   clear_mat(state->box);
 507   clear_mat(state->box_rel);
 508   clear_mat(state->boxv);
 509   clear_mat(state->pres_prev);
 510   clear_mat(state->svir_prev);
 511   clear_mat(state->fvir_prev);
 512   init_gtc_state(state,ngtc,nnhpres,nhchainlength);
 513   state->nalloc = state->natoms;
 514   if (state->nalloc > 0) {
 515     snew(state->x,state->nalloc);
 516     snew(state->v,state->nalloc);
 517   } else {
 518     state->x = NULL;
 519     state->v = NULL;
 520   }
 521   state->sd_X = NULL;
 522   state->cg_p = NULL;
 523
 524   init_ekinstate(&state->ekinstate);
 525
 526   init_energyhistory(&state->enerhist);
 527
 528   state->ddp_count = 0;
 529   state->ddp_count_cg_gl = 0;
 530   state->cg_gl = NULL;
 531   state->cg_gl_nalloc = 0;
 532 }
 533
 534 void done_state(t_state *state)
 535 {
 536   if (state->nosehoover_xi) sfree(state->nosehoover_xi);
 537   if (state->x) sfree(state->x);
 538   if (state->v) sfree(state->v);
 539   if (state->sd_X) sfree(state->sd_X);
 540   if (state->cg_p) sfree(state->cg_p);
 541   state->nalloc = 0;
 542   if (state->cg_gl) sfree(state->cg_gl);
 543   state->cg_gl_nalloc = 0;
 544 }
 545
 546 static void do_box_rel(t_inputrec *ir,matrix box_rel,matrix b,bool bInit)
 547 {
 548   int d,d2;
 549
 550   for(d=YY; d<=ZZ; d++) {
 551     for(d2=XX; d2<=(ir->epct==epctSEMIISOTROPIC ? YY : ZZ); d2++) {
 552       /* We need to check if this box component is deformed
 553        * or if deformation of another component might cause
 554        * changes in this component due to box corrections.
 555        */
 556       if (ir->deform[d][d2] == 0 &&
 557           !(d == ZZ && d2 == XX && ir->deform[d][YY] != 0 &&
 558             (b[YY][d2] != 0 || ir->deform[YY][d2] != 0))) {
 559         if (bInit) {
 560           box_rel[d][d2] = b[d][d2]/b[XX][XX];
 561         } else {
 562           b[d][d2] = b[XX][XX]*box_rel[d][d2];
 563         }
 564       }
 565     }
 566   }
 567 }
 568
 569 void set_box_rel(t_inputrec *ir,t_state *state)
 570 {
 571   /* Make sure the box obeys the restrictions before we fix the ratios */
 572   correct_box(NULL,0,state->box,NULL);
 573
 574   clear_mat(state->box_rel);
 575
 576   if (PRESERVE_SHAPE(*ir))
 577     do_box_rel(ir,state->box_rel,state->box,TRUE);
 578 }
 579
 580 void preserve_box_shape(t_inputrec *ir,matrix box_rel,matrix b)
 581 {
 582   if (PRESERVE_SHAPE(*ir))
 583     do_box_rel(ir,box_rel,b,FALSE);
 584 }
 585
 586 void add_t_atoms(t_atoms *atoms,int natom_extra,int nres_extra)
 587 {
 588     int i;
 589
 590     if (natom_extra > 0)
 591     {
 592         srenew(atoms->atomname,atoms->nr+natom_extra);
 593         srenew(atoms->atom,atoms->nr+natom_extra);
 594         if (NULL != atoms->pdbinfo)
 595             srenew(atoms->pdbinfo,atoms->nr+natom_extra);
 596         if (NULL != atoms->atomtype)
 597             srenew(atoms->atomtype,atoms->nr+natom_extra);
 598         if (NULL != atoms->atomtypeB)
 599             srenew(atoms->atomtypeB,atoms->nr+natom_extra);
 600         for(i=atoms->nr; (i<atoms->nr+natom_extra); i++) {
 601             atoms->atomname[i] = NULL;
 602             memset(&atoms->atom[i],0,sizeof(atoms->atom[i]));
 603             if (NULL != atoms->pdbinfo)
 604                 memset(&atoms->pdbinfo[i],0,sizeof(atoms->pdbinfo[i]));
 605             if (NULL != atoms->atomtype)
 606                 atoms->atomtype[i] = NULL;
 607             if (NULL != atoms->atomtypeB)
 608                 atoms->atomtypeB[i] = NULL;
 609         }
 610         atoms->nr += natom_extra;
 611     }
 612     if (nres_extra > 0)
 613     {
 614         srenew(atoms->resinfo,atoms->nres+nres_extra);
 615         for(i=atoms->nres; (i<atoms->nres+nres_extra); i++) {
 616             memset(&atoms->resinfo[i],0,sizeof(atoms->resinfo[i]));
 617         }
 618         atoms->nres += nres_extra;
 619     }
 620 }
 621
 622 void init_t_atoms(t_atoms *atoms, int natoms, bool bPdbinfo)
 623 {
 624   atoms->nr=natoms;
 625   atoms->nres=0;
 626   snew(atoms->atomname,natoms);
 627   atoms->atomtype=NULL;
 628   atoms->atomtypeB=NULL;
 629   snew(atoms->resinfo,natoms);
 630   snew(atoms->atom,natoms);
 631   if (bPdbinfo)
 632     snew(atoms->pdbinfo,natoms);
 633   else
 634     atoms->pdbinfo=NULL;
 635 }
 636
 637 t_atoms *copy_t_atoms(t_atoms *src)
 638 {
 639   t_atoms *dst;
 640   int i;
 641
 642   snew(dst,1);
 643   init_t_atoms(dst,src->nr,(NULL != src->pdbinfo));
 644   dst->nr = src->nr;
 645   if (NULL != src->atomname)
 646       snew(dst->atomname,src->nr);
 647   if (NULL != src->atomtype)
 648       snew(dst->atomtype,src->nr);
 649   if (NULL != src->atomtypeB)
 650       snew(dst->atomtypeB,src->nr);
 651   for(i=0; (i<src->nr); i++) {
 652     dst->atom[i] = src->atom[i];
 653     if (NULL != src->pdbinfo)
 654       dst->pdbinfo[i] = src->pdbinfo[i];
 655     if (NULL != src->atomname)
 656         dst->atomname[i]  = src->atomname[i];
 657     if (NULL != src->atomtype)
 658         dst->atomtype[i] = src->atomtype[i];
 659     if (NULL != src->atomtypeB)
 660         dst->atomtypeB[i] = src->atomtypeB[i];
 661   }
 662   dst->nres = src->nres;
 663   for(i=0; (i<src->nres); i++) {
 664     dst->resinfo[i] = src->resinfo[i];
 665   }
 666   return dst;
 667 }
 668
 669 void t_atoms_set_resinfo(t_atoms *atoms,int atom_ind,t_symtab *symtab,
 670                          const char *resname,int resnr,unsigned char ic,
 671                          int chainnum, char chainid)
 672 {
 673   t_resinfo *ri;
 674
 675   ri = &atoms->resinfo[atoms->atom[atom_ind].resind];
 676   ri->name  = put_symtab(symtab,resname);
 677   ri->rtp   = NULL;
 678   ri->nr    = resnr;
 679   ri->ic    = ic;
 680   ri->chainnum = chainnum;
 681   ri->chainid = chainid;
 682 }
 683
 684 void free_t_atoms(t_atoms *atoms,bool bFreeNames)
 685 {
 686   int i;
 687
 688   if (bFreeNames) {
 689     for(i=0; i<atoms->nr; i++) {
 690       sfree(*atoms->atomname[i]);
 691       *atoms->atomname[i]=NULL;
 692     }
 693     for(i=0; i<atoms->nres; i++) {
 694       sfree(*atoms->resinfo[i].name);
 695       *atoms->resinfo[i].name=NULL;
 696     }
 697   }
 698   sfree(atoms->atomname);
 699   /* Do we need to free atomtype and atomtypeB as well ? */
 700   sfree(atoms->resinfo);
 701   sfree(atoms->atom);
 702   if (atoms->pdbinfo)
 703     sfree(atoms->pdbinfo);
 704   atoms->nr=0;
 705   atoms->nres=0;
 706 }
 707
 708 real max_cutoff(real cutoff1,real cutoff2)
 709 {
 710     if (cutoff1 == 0 || cutoff2 == 0)
 711     {
 712         return 0;
 713     }
 714     else
 715     {
 716         return max(cutoff1,cutoff2);
 717     }
 718 }