src/kernel/tomorse.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <stdlib.h>
  44 #include <math.h>
  45 #include <ctype.h>
  46 #include "typedefs.h"
  47 #include "string2.h"
  48 #include "grompp.h"
  49 #include "futil.h"
  50 #include "smalloc.h"
  51 #include "toputil.h"
  52 #include "gmx_fatal.h"
  53 #include "gpp_atomtype.h"
  54 #include "gpp_tomorse.h"
  55
  56 typedef struct {
  57   char *ai,*aj;
  58   real e_diss;
  59 } t_2morse;
  60
  61 static t_2morse *read_dissociation_energies(int *n2morse)
  62 {
  63   FILE     *fp;
  64   char     ai[32],aj[32];
  65   double   e_diss;
  66   const char *fn="edissoc.dat";
  67   t_2morse *t2m=NULL;
  68   int      maxn2m=0,n2m=0;
  69   int      nread;
  70
  71   /* Open the file with dissociation energies */
  72   fp = libopen(fn);
  73   do {
  74     /* Try and read two atom names and an energy term from it */
  75     nread = fscanf(fp,"%s%s%lf",ai,aj,&e_diss);
  76     if (nread == 3) {
  77       /* If we got three terms, it probably was OK, no further checking */
  78       if (n2m >= maxn2m) {
  79         /* Increase memory for 16 at once, some mallocs are stupid */
  80         maxn2m += 16;
  81         srenew(t2m,maxn2m);
  82       }
  83       /* Copy the values */
  84       t2m[n2m].ai = strdup(ai);
  85       t2m[n2m].aj = strdup(aj);
  86       t2m[n2m].e_diss = e_diss;
  87       /* Increment counter */
  88       n2m++;
  89     }
  90     /* If we did not read three items, quit reading */
  91   } while (nread == 3);
  92   ffclose(fp);
  93
  94   /* Set the return values */
  95   *n2morse = n2m;
  96
  97   return t2m;
  98 }
  99
 100 static int nequal(char *a1,char *a2)
 101 {
 102   int i;
 103
 104   /* Count the number of (case insensitive) characters that are equal in
 105    * two strings. If they are equally long their respective null characters are
 106    * counted also.
 107    */
 108   for(i=0; (a1[i] != '\0') && (a2[i] != '\0'); i++)
 109     if (toupper(a1[i]) != toupper(a2[i]))
 110       break;
 111   if ((a1[i] == '\0') && (a2[i] == '\0'))
 112     i++;
 113
 114   return i;
 115 }
 116
 117 static real search_e_diss(int n2m,t_2morse t2m[],char *ai,char *aj)
 118 {
 119   int i;
 120   int ibest=-1;
 121   int nii,njj,nbstii=0,nbstjj=0;
 122   real ediss = 400;
 123
 124   /* Do a best match search for dissociation energies */
 125   for(i=0; (i<n2m); i++) {
 126     /* Check for a perfect match */
 127     if (((gmx_strcasecmp(t2m[i].ai,ai) == 0) && (gmx_strcasecmp(t2m[i].aj,aj) == 0)) ||
 128         ((gmx_strcasecmp(t2m[i].aj,ai) == 0) && (gmx_strcasecmp(t2m[i].ai,aj) == 0))) {
 129       ibest = i;
 130       break;
 131     }
 132     else {
 133       /* Otherwise count the number of equal characters in the strings ai and aj
 134        * and the ones from the file
 135        */
 136       nii = nequal(t2m[i].ai,ai);
 137       njj = nequal(t2m[i].aj,aj);
 138       if (((nii >  nbstii) && (njj >= nbstjj)) ||
 139           ((nii >= nbstii) && (njj >  nbstjj))) {
 140         if ((nii > 0) && (njj > 0)) {
 141           ibest  = i;
 142           nbstii = nii;
 143           nbstjj = njj;
 144         }
 145       }
 146       else {
 147         /* Swap ai and aj (at least in counting the number of equal chars) */
 148         nii = nequal(t2m[i].ai,aj);
 149         njj = nequal(t2m[i].aj,ai);
 150         if (((nii >  nbstii) && (njj >= nbstjj)) ||
 151             ((nii >= nbstii) && (njj >  nbstjj))) {
 152           if ((nii > 0) && (njj > 0)) {
 153             ibest  = i;
 154             nbstii = nii;
 155             nbstjj = njj;
 156           }
 157         }
 158       }
 159     }
 160   }
 161   /* Return the dissocation energy corresponding to the best match, if we have
 162    * found one. Do some debug output anyway.
 163    */
 164   if (ibest == -1) {
 165     if (debug)
 166       fprintf(debug,"MORSE: Couldn't find E_diss for bond %s - %s, using default %g\n",ai,aj,ediss);
 167     return ediss;
 168   }
 169   else {
 170     if (debug)
 171       fprintf(debug,"MORSE: Dissoc. E (%10.3f) for bond %4s-%4s taken from bond %4s-%4s\n",
 172               t2m[ibest].e_diss,ai,aj,t2m[ibest].ai,t2m[ibest].aj);
 173     return t2m[ibest].e_diss;
 174   }
 175 }
 176
 177 void convert_harmonics(int nrmols,t_molinfo mols[],gpp_atomtype_t atype)
 178 {
 179   int      n2m;
 180   t_2morse *t2m;
 181
 182   int  i,j,k,last,ni,nj;
 183   int  nrharm,nrmorse,bb;
 184   real edis,kb,b0,beta;
 185   gmx_bool *bRemoveHarm;
 186
 187   /* First get the data */
 188   t2m = read_dissociation_energies(&n2m);
 189   if (debug)
 190     fprintf(debug,"MORSE: read %d dissoc energies\n",n2m);
 191   if (n2m <= 0) {
 192     fprintf(stderr,"No dissocation energies read\n");
 193     return;
 194   }
 195
 196   /* For all the molecule types */
 197   for(i=0; (i<nrmols); i++) {
 198     /* Check how many morse and harmonic BONDSs there are, increase size of
 199      * morse with the number of harmonics
 200      */
 201     nrmorse = mols[i].plist[F_MORSE].nr;
 202
 203     for(bb=0; (bb < F_NRE); bb++) {
 204       if ((interaction_function[bb].flags & IF_BTYPE) && (bb != F_MORSE)) {
 205         nrharm  = mols[i].plist[bb].nr;
 206         pr_alloc(nrharm,&(mols[i].plist[F_MORSE]));
 207         snew(bRemoveHarm,nrharm);
 208
 209         /* Now loop over the harmonics, trying to convert them */
 210         for(j=0; (j<nrharm); j++) {
 211           ni   = mols[i].plist[bb].param[j].AI;
 212           nj   = mols[i].plist[bb].param[j].AJ;
 213           edis =
 214             search_e_diss(n2m,t2m,
 215                           get_atomtype_name(mols[i].atoms.atom[ni].type,atype),
 216                           get_atomtype_name(mols[i].atoms.atom[nj].type,atype));
 217           if (edis != 0) {
 218             bRemoveHarm[j] = TRUE;
 219             b0   = mols[i].plist[bb].param[j].c[0];
 220             kb   = mols[i].plist[bb].param[j].c[1];
 221             beta = sqrt(kb/(2*edis));
 222             mols[i].plist[F_MORSE].param[nrmorse].a[0] = ni;
 223             mols[i].plist[F_MORSE].param[nrmorse].a[1] = nj;
 224             mols[i].plist[F_MORSE].param[nrmorse].c[0] = b0;
 225             mols[i].plist[F_MORSE].param[nrmorse].c[1] = edis;
 226             mols[i].plist[F_MORSE].param[nrmorse].c[2] = beta;
 227             nrmorse++;
 228           }
 229         }
 230         mols[i].plist[F_MORSE].nr = nrmorse;
 231
 232         /* Now remove the harmonics */
 233         for(j=last=0; (j<nrharm); j++) {
 234           if (!bRemoveHarm[j]) {
 235             /* Copy it to the last position */
 236             for(k=0; (k<MAXATOMLIST); k++)
 237               mols[i].plist[bb].param[last].a[k] =
 238                 mols[i].plist[bb].param[j].a[k];
 239             for(k=0; (k<MAXFORCEPARAM); k++)
 240               mols[i].plist[bb].param[last].c[k] =
 241                 mols[i].plist[bb].param[j].c[k];
 242             last++;
 243           }
 244         }
 245         sfree(bRemoveHarm);
 246         fprintf(stderr,"Converted %d out of %d %s to morse bonds for mol %d\n",
 247                 nrharm-last,nrharm,interaction_function[bb].name,i);
 248         mols[i].plist[bb].nr = last;
 249       }
 250     }
 251   }
 252   sfree(t2m);
 253 }