src/kernel/tomorse.c

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