src/kernel/tomorse.c

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