src/gmxlib/topsort.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-  */
   2
   3 /*
   4  *
   5  *                This source code is part of
   6  *
   7  *                 G   R   O   M   A   C   S
   8  *
   9  *          GROningen MAchine for Chemical Simulations
  10  *
  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-2008, 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  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43
  44 static bool ip_pert(int ftype,t_iparams *ip)
  45 {
  46     bool bPert;
  47     int  i;
  48
  49     if (NRFPB(ftype) == 0)
  50     {
  51         return FALSE;
  52     }
  53
  54     switch (ftype)
  55     {
  56     case F_BONDS:
  57     case F_G96BONDS:
  58     case F_HARMONIC:
  59     case F_ANGLES:
  60     case F_G96ANGLES:
  61     case F_IDIHS:
  62         bPert = (ip->harmonic.rA  != ip->harmonic.rB ||
  63                  ip->harmonic.krA != ip->harmonic.krB);
  64         break;
  65     case F_PDIHS:
  66     case F_ANGRES:
  67     case F_ANGRESZ:
  68         bPert = (ip->pdihs.phiA != ip->pdihs.phiB ||
  69                  ip->pdihs.cpA  != ip->pdihs.cpB);
  70         break;
  71     case F_RBDIHS:
  72         bPert = FALSE;
  73         for(i=0; i<NR_RBDIHS; i++)
  74         {
  75             if (ip->rbdihs.rbcA[i] != ip->rbdihs.rbcB[i])
  76             {
  77                 bPert = TRUE;
  78             }
  79         }
  80         break;
  81     case F_TABBONDS:
  82     case F_TABBONDSNC:
  83     case F_TABANGLES:
  84     case F_TABDIHS:
  85         bPert = (ip->tab.kA != ip->tab.kB);
  86         break;
  87     case F_POSRES:
  88         bPert = FALSE;
  89         for(i=0; i<DIM; i++)
  90         {
  91             if (ip->posres.pos0A[i] != ip->posres.pos0B[i] ||
  92                 ip->posres.fcA[i]   != ip->posres.fcB[i])
  93             {
  94                 bPert = TRUE;
  95             }
  96         }
  97         break;
  98     case F_LJ14:
  99         bPert = (ip->lj14.c6A  != ip->lj14.c6B ||
 100                  ip->lj14.c12A != ip->lj14.c12B);
 101         break;
 102     default:
 103         bPert = FALSE;
 104         gmx_fatal(FARGS,"Function type %s not implemented in ip_pert",
 105                   interaction_function[ftype].longname);
 106     }
 107
 108     return bPert;
 109 }
 110
 111 bool gmx_mtop_bondeds_free_energy(const gmx_mtop_t *mtop)
 112 {
 113     const gmx_ffparams_t *ffparams;
 114     int  i,ftype;
 115     bool bPert;
 116
 117     ffparams = &mtop->ffparams;
 118
 119     /* Loop over all the function types and compare the A/B parameters */
 120     bPert = FALSE;
 121     for(i=0; i<ffparams->ntypes; i++)
 122     {
 123         ftype = ffparams->functype[i];
 124         if (interaction_function[ftype].flags & IF_BOND)
 125         {
 126             if (ip_pert(ftype,&ffparams->iparams[i]))
 127             {
 128                 bPert = TRUE;
 129             }
 130         }
 131     }
 132
 133     return (bPert ? ilsortFE_UNSORTED : ilsortNO_FE);
 134 }
 135
 136 void gmx_sort_ilist_fe(t_idef *idef)
 137 {
 138     int  ftype,nral,i,ic,ib,a;
 139     t_iparams *iparams;
 140     t_ilist *ilist;
 141     t_iatom *iatoms;
 142     bool bPert;
 143     t_iatom *iabuf;
 144     int  iabuf_nalloc;
 145
 146     iabuf_nalloc = 0;
 147     iabuf        = NULL;
 148
 149     iparams = idef->iparams;
 150
 151     for(ftype=0; ftype<F_NRE; ftype++)
 152     {
 153         if (interaction_function[ftype].flags & IF_BOND)
 154         {
 155             ilist = &idef->il[ftype];
 156             iatoms = ilist->iatoms;
 157             nral  = NRAL(ftype);
 158             ic = 0;
 159             ib = 0;
 160             i  = 0;
 161             while (i < ilist->nr)
 162             {
 163                 /* The first element of ia gives the type */
 164                 if (ip_pert(ftype,&iparams[iatoms[i]]))
 165                 {
 166                     /* Copy to the perturbed buffer */
 167                     if (ib + 1 + nral > iabuf_nalloc)
 168                     {
 169                         iabuf_nalloc = over_alloc_large(ib+1+nral);
 170                         srenew(iabuf,iabuf_nalloc);
 171                     }
 172                     for(a=0; a<1+nral; a++)
 173                     {
 174                         iabuf[ib++] = iatoms[i++];
 175                     }
 176                 }
 177                 else
 178                 {
 179                     /* Copy in place */
 180                     for(a=0; a<1+nral; a++)
 181                     {
 182                         iatoms[ic++] = iatoms[i++];
 183                     }
 184                 }
 185             }
 186             /* Now we now the number of non-perturbed interactions */
 187             ilist->nr_nonperturbed = ic;
 188
 189             /* Copy the buffer with perturbed interactions to the ilist */
 190             for(a=0; a<ib; a++)
 191             {
 192                 iatoms[ic++] = iabuf[a];
 193             }
 194
 195             if (debug)
 196             {
 197                 fprintf(debug,"%s non-pert %d pert %d\n",
 198                         interaction_function[ftype].longname,
 199                         ilist->nr_nonperturbed,
 200                         ilist->nr-ilist->nr_nonperturbed);
 201             }
 202         }
 203     }
 204
 205     sfree(iabuf);
 206
 207     idef->ilsort = ilsortFE_SORTED;
 208 }