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  * This file is part of the GROMACS molecular simulation package.
   5  *
   6  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   7  * Copyright (c) 2001-2008, The GROMACS development team,
   8  * check out http://www.gromacs.org for more information.
   9  * Copyright (c) 2012,2013, by the GROMACS development team, led by
  10  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
  11  * others, as listed in the AUTHORS file in the top-level source
  12  * directory and at http://www.gromacs.org.
  13  *
  14  * GROMACS is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU Lesser General Public License
  16  * as published by the Free Software Foundation; either version 2.1
  17  * of the License, or (at your option) any later version.
  18  *
  19  * GROMACS is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  22  * Lesser General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Lesser General Public
  25  * License along with GROMACS; if not, see
  26  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  27  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  28  *
  29  * If you want to redistribute modifications to GROMACS, please
  30  * consider that scientific software is very special. Version
  31  * control is crucial - bugs must be traceable. We will be happy to
  32  * consider code for inclusion in the official distribution, but
  33  * derived work must not be called official GROMACS. Details are found
  34  * in the README & COPYING files - if they are missing, get the
  35  * official version at http://www.gromacs.org.
  36  *
  37  * To help us fund GROMACS development, we humbly ask that you cite
  38  * the research papers on the package. Check out http://www.gromacs.org.
  39  */
  40
  41 #ifdef HAVE_CONFIG_H
  42 #include <config.h>
  43 #endif
  44
  45 #include "typedefs.h"
  46 #include "topsort.h"
  47 #include "smalloc.h"
  48 #include "gmx_fatal.h"
  49
  50 static gmx_bool ip_pert(int ftype, const t_iparams *ip)
  51 {
  52     gmx_bool bPert;
  53     int      i;
  54
  55     if (NRFPB(ftype) == 0)
  56     {
  57         return FALSE;
  58     }
  59
  60     switch (ftype)
  61     {
  62         case F_BONDS:
  63         case F_G96BONDS:
  64         case F_HARMONIC:
  65         case F_ANGLES:
  66         case F_G96ANGLES:
  67         case F_IDIHS:
  68             bPert = (ip->harmonic.rA  != ip->harmonic.rB ||
  69                      ip->harmonic.krA != ip->harmonic.krB);
  70             break;
  71         case F_MORSE:
  72             bPert = (ip->morse.b0A  != ip->morse.b0B ||
  73                      ip->morse.cbA  != ip->morse.cbB ||
  74                      ip->morse.betaA  != ip->morse.betaB);
  75             break;
  76         case F_RESTRBONDS:
  77             bPert = (ip->restraint.lowA  != ip->restraint.lowB ||
  78                      ip->restraint.up1A  != ip->restraint.up1B ||
  79                      ip->restraint.up2A  != ip->restraint.up2B ||
  80                      ip->restraint.kA    != ip->restraint.kB);
  81             break;
  82         case F_UREY_BRADLEY:
  83             bPert = (ip->u_b.thetaA  != ip->u_b.thetaB  ||
  84                      ip->u_b.kthetaA != ip->u_b.kthetaB ||
  85                      ip->u_b.r13A    != ip->u_b.r13B    ||
  86                      ip->u_b.kUBA    != ip->u_b.kUBB);
  87             break;
  88         case F_PDIHS:
  89         case F_PIDIHS:
  90         case F_ANGRES:
  91         case F_ANGRESZ:
  92             bPert = (ip->pdihs.phiA != ip->pdihs.phiB ||
  93                      ip->pdihs.cpA  != ip->pdihs.cpB);
  94             break;
  95         case F_RBDIHS:
  96             bPert = FALSE;
  97             for (i = 0; i < NR_RBDIHS; i++)
  98             {
  99                 if (ip->rbdihs.rbcA[i] != ip->rbdihs.rbcB[i])
 100                 {
 101                     bPert = TRUE;
 102                 }
 103             }
 104             break;
 105         case F_TABBONDS:
 106         case F_TABBONDSNC:
 107         case F_TABANGLES:
 108         case F_TABDIHS:
 109             bPert = (ip->tab.kA != ip->tab.kB);
 110             break;
 111         case F_POSRES:
 112             bPert = FALSE;
 113             for (i = 0; i < DIM; i++)
 114             {
 115                 if (ip->posres.pos0A[i] != ip->posres.pos0B[i] ||
 116                     ip->posres.fcA[i]   != ip->posres.fcB[i])
 117                 {
 118                     bPert = TRUE;
 119                 }
 120             }
 121             break;
 122         case F_DIHRES:
 123             bPert = ((ip->dihres.phiA != ip->dihres.phiB) ||
 124                      (ip->dihres.dphiA != ip->dihres.dphiB) ||
 125                      (ip->dihres.kfacA != ip->dihres.kfacB));
 126             break;
 127         case F_LJ14:
 128             bPert = (ip->lj14.c6A  != ip->lj14.c6B ||
 129                      ip->lj14.c12A != ip->lj14.c12B);
 130             break;
 131         case F_CMAP:
 132             bPert = FALSE;
 133             break;
 134         default:
 135             bPert = FALSE;
 136             gmx_fatal(FARGS, "Function type %s not implemented in ip_pert",
 137                       interaction_function[ftype].longname);
 138     }
 139
 140     return bPert;
 141 }
 142
 143 static gmx_bool ip_q_pert(int ftype, const t_iatom *ia,
 144                           const t_iparams *ip, const real *qA, const real *qB)
 145 {
 146     /* 1-4 interactions do not have the charges stored in the iparams list,
 147      * so we need a separate check for those.
 148      */
 149     return (ip_pert(ftype, ip+ia[0]) ||
 150             (ftype == F_LJ14 && (qA[ia[1]] != qB[ia[1]] ||
 151                                  qA[ia[2]] != qB[ia[2]])));
 152 }
 153
 154 gmx_bool gmx_mtop_bondeds_free_energy(const gmx_mtop_t *mtop)
 155 {
 156     const gmx_ffparams_t *ffparams;
 157     int                   i, ftype;
 158     int                   mb;
 159     t_atom               *atom;
 160     t_ilist              *il;
 161     t_iatom              *ia;
 162     gmx_bool              bPert;
 163
 164     ffparams = &mtop->ffparams;
 165
 166     /* Loop over all the function types and compare the A/B parameters */
 167     bPert = FALSE;
 168     for (i = 0; i < ffparams->ntypes; i++)
 169     {
 170         ftype = ffparams->functype[i];
 171         if (interaction_function[ftype].flags & IF_BOND)
 172         {
 173             if (ip_pert(ftype, &ffparams->iparams[i]))
 174             {
 175                 bPert = TRUE;
 176             }
 177         }
 178     }
 179
 180     /* Check perturbed charges for 1-4 interactions */
 181     for (mb = 0; mb < mtop->nmolblock; mb++)
 182     {
 183         atom = mtop->moltype[mtop->molblock[mb].type].atoms.atom;
 184         il   = &mtop->moltype[mtop->molblock[mb].type].ilist[F_LJ14];
 185         ia   = il->iatoms;
 186         for (i = 0; i < il->nr; i += 3)
 187         {
 188             if (atom[ia[i+1]].q != atom[ia[i+1]].qB ||
 189                 atom[ia[i+2]].q != atom[ia[i+2]].qB)
 190             {
 191                 bPert = TRUE;
 192             }
 193         }
 194     }
 195
 196     return bPert;
 197 }
 198
 199 void gmx_sort_ilist_fe(t_idef *idef, const real *qA, const real *qB)
 200 {
 201     int        ftype, nral, i, ic, ib, a;
 202     t_iparams *iparams;
 203     t_ilist   *ilist;
 204     t_iatom   *iatoms;
 205     gmx_bool   bPert;
 206     t_iatom   *iabuf;
 207     int        iabuf_nalloc;
 208
 209     if (qB == NULL)
 210     {
 211         qB = qA;
 212     }
 213
 214     iabuf_nalloc = 0;
 215     iabuf        = NULL;
 216
 217     iparams = idef->iparams;
 218
 219     for (ftype = 0; ftype < F_NRE; ftype++)
 220     {
 221         if (interaction_function[ftype].flags & IF_BOND)
 222         {
 223             ilist  = &idef->il[ftype];
 224             iatoms = ilist->iatoms;
 225             nral   = NRAL(ftype);
 226             ic     = 0;
 227             ib     = 0;
 228             i      = 0;
 229             while (i < ilist->nr)
 230             {
 231                 /* Check if this interaction is perturbed */
 232                 if (ip_q_pert(ftype, iatoms+i, iparams, qA, qB))
 233                 {
 234                     /* Copy to the perturbed buffer */
 235                     if (ib + 1 + nral > iabuf_nalloc)
 236                     {
 237                         iabuf_nalloc = over_alloc_large(ib+1+nral);
 238                         srenew(iabuf, iabuf_nalloc);
 239                     }
 240                     for (a = 0; a < 1+nral; a++)
 241                     {
 242                         iabuf[ib++] = iatoms[i++];
 243                     }
 244                 }
 245                 else
 246                 {
 247                     /* Copy in place */
 248                     for (a = 0; a < 1+nral; a++)
 249                     {
 250                         iatoms[ic++] = iatoms[i++];
 251                     }
 252                 }
 253             }
 254             /* Now we now the number of non-perturbed interactions */
 255             ilist->nr_nonperturbed = ic;
 256
 257             /* Copy the buffer with perturbed interactions to the ilist */
 258             for (a = 0; a < ib; a++)
 259             {
 260                 iatoms[ic++] = iabuf[a];
 261             }
 262
 263             if (debug)
 264             {
 265                 fprintf(debug, "%s non-pert %d pert %d\n",
 266                         interaction_function[ftype].longname,
 267                         ilist->nr_nonperturbed,
 268                         ilist->nr-ilist->nr_nonperturbed);
 269             }
 270         }
 271     }
 272
 273     sfree(iabuf);
 274
 275     idef->ilsort = ilsortFE_SORTED;
 276 }