src/kernel/convparm.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_convparm_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "sysstuff.h"
  33 #include "physics.h"
  34 #include "vec.h"
  35 #include "assert.h"
  36 #include "smalloc.h"
  37 #include "typedefs.h"
  38 #include "fatal.h"
  39 #include "topio.h"
  40 #include "toputil.h"
  41 #include "convparm.h"
  42 #include "names.h"
  43
  44 static void assign_param(t_functype ftype,t_iparams *new,
  45                          real old[MAXFORCEPARAM])
  46 {
  47   int i,j;
  48
  49   /* Set to zero */
  50   for(j=0; (j<MAXFORCEPARAM); j++)
  51     new->generic.buf[j]=0.0;
  52
  53   switch (ftype) {
  54   case F_G96ANGLES:
  55     /* Post processing of input data: store cosine iso angle itself */
  56     new->harmonic.rA =cos(old[0]*DEG2RAD);
  57     new->harmonic.krA=old[1];
  58     new->harmonic.rB =cos(old[2]*DEG2RAD);
  59     new->harmonic.krB=old[3];
  60     break;
  61   case F_G96BONDS:
  62     /* Post processing of input data: store square of length itself */
  63     new->harmonic.rA =sqr(old[0]);
  64     new->harmonic.krA=old[1];
  65     new->harmonic.rB =sqr(old[2]);
  66     new->harmonic.krB=old[3];
  67     break;
  68   case F_ANGLES:
  69   case F_BONDS:
  70   case F_IDIHS:
  71     new->harmonic.rA =old[0];
  72     new->harmonic.krA=old[1];
  73     new->harmonic.rB =old[2];
  74     new->harmonic.krB=old[3];
  75     break;
  76   case F_MORSE:
  77     new->morse.b0    =old[0];
  78     new->morse.cb    =old[1];
  79     new->morse.beta  =old[2];
  80     break;
  81   case F_WPOL:
  82     new->wpol.kx     =old[0];
  83     new->wpol.ky     =old[1];
  84     new->wpol.kz     =old[2];
  85     new->wpol.rOH    =old[3];
  86     new->wpol.rHH    =old[4];
  87     new->wpol.rOD    =old[5];
  88     break;
  89   case F_BHAM:
  90     new->bham.a = old[0];
  91     new->bham.b = old[1];
  92     new->bham.c = old[2];
  93     break;
  94   case F_LJ14:
  95     new->lj14.c6A  = old[0];
  96     new->lj14.c12A = old[1];
  97     new->lj14.c6B  = old[2];
  98     new->lj14.c12B = old[3];
  99     break;
 100   case F_LJ:
 101     new->lj.c6  = old[0];
 102     new->lj.c12 = old[1];
 103     break;
 104   case F_PDIHS:
 105   case F_ANGRES:
 106   case F_ANGRESZ:
 107     new->pdihs.phiA=old[0];
 108     new->pdihs.cpA =old[1];
 109     new->pdihs.mult=old[2];
 110     new->pdihs.phiB=old[3];
 111     new->pdihs.cpB =old[4];
 112     break;
 113   case F_POSRES:
 114     new->posres.fc[XX]   = old[0];
 115     new->posres.fc[YY]   = old[1];
 116     new->posres.fc[ZZ]   = old[2];
 117     new->posres.pos0[XX] = old[3];
 118     new->posres.pos0[YY] = old[4];
 119     new->posres.pos0[ZZ] = old[5];
 120     break;
 121   case F_DISRES:
 122     new->disres.index=old[0];
 123     new->disres.type=old[1];
 124     new->disres.low=old[2];
 125     new->disres.up1=old[3];
 126     new->disres.up2=old[4];
 127     new->disres.fac=old[5];
 128     break;
 129   case F_RBDIHS:
 130     for (i=0; (i<NR_RBDIHS); i++)
 131       new->rbdihs.rbc[i]=old[i];
 132     break;
 133   case F_SHAKE:
 134   case F_SHAKENC:
 135     new->shake.dA = old[0];
 136     new->shake.dB = old[1];
 137     break;
 138   case F_SETTLE:
 139     new->settle.doh=old[0];
 140     new->settle.dhh=old[1];
 141     break;
 142   case F_DUMMY2:
 143   case F_DUMMY3:
 144   case F_DUMMY3FD:
 145   case F_DUMMY3OUT:
 146   case F_DUMMY4FD:
 147     new->dummy.a=old[0];
 148     new->dummy.b=old[1];
 149     new->dummy.c=old[2];
 150     new->dummy.d=old[3];
 151     new->dummy.e=old[4];
 152     new->dummy.f=old[5];
 153     break;
 154   case F_DUMMY3FAD:
 155     new->dummy.a=old[1] * cos(DEG2RAD * old[0]);
 156     new->dummy.b=old[1] * sin(DEG2RAD * old[0]);
 157     new->dummy.c=old[2];
 158     new->dummy.d=old[3];
 159     new->dummy.e=old[4];
 160     new->dummy.f=old[5];
 161     break;
 162   default:
 163     fatal_error(0,"unknown function type %d in %s line %d",
 164                 ftype,__FILE__,__LINE__);
 165   }
 166 }
 167
 168 static int enter_params(t_idef *idef, t_functype ftype,
 169                         real forceparams[MAXFORCEPARAM],int start,bool bAppend)
 170 {
 171   t_iparams new;
 172   int       type;
 173
 174   assign_param(ftype,&new,forceparams);
 175   if (!bAppend) {
 176     for (type=start; (type<idef->ntypes); type++) {
 177       if (idef->functype[type]==ftype) {
 178         if (memcmp(&new,&idef->iparams[type],(size_t)sizeof(new)) == 0)
 179           return type;
 180       }
 181     }
 182   }
 183   else
 184     type=idef->ntypes;
 185   if (debug)
 186     fprintf(debug,"copying new to idef->iparams[%d] (ntypes=%d)\n",
 187             type,idef->ntypes);
 188   memcpy(&idef->iparams[type],&new,(size_t)sizeof(new));
 189
 190   idef->ntypes++;
 191   idef->functype[type]=ftype;
 192
 193   return type;
 194 }
 195
 196 static void append_interaction(t_ilist *ilist,
 197                                int type,int nral,atom_id a[MAXATOMLIST])
 198 {
 199   int i,where1;
 200
 201   where1     = ilist->nr;
 202   ilist->nr += nral+1;
 203
 204   ilist->iatoms[where1++]=type;
 205   for (i=0; (i<nral); i++)
 206     ilist->iatoms[where1++]=a[i];
 207 }
 208
 209 static void enter_function(t_params *p,t_functype ftype,
 210                            t_idef *idef,int *maxtypes,bool bNB,bool bAppend)
 211 {
 212   int     k,type,nr,nral,delta,start;
 213   t_ilist *il;
 214
 215   il    = &(idef->il[ftype]);
 216   start = idef->ntypes;
 217   nr    = p->nr;
 218   nral  = NRAL(ftype);
 219   delta = nr*(nral+1);
 220   srenew(il->iatoms,il->nr+delta);
 221
 222   for (k=0; k<nr; k++) {
 223     if (*maxtypes <= idef->ntypes) {
 224       *maxtypes += 1000;
 225       srenew(idef->functype,*maxtypes);
 226       srenew(idef->iparams, *maxtypes);
 227       if (debug)
 228         fprintf(debug,"%s, line %d: srenewed idef->functype and idef->iparams to %d\n",
 229                 __FILE__,__LINE__,*maxtypes);
 230     }
 231     type = enter_params(idef,ftype,p->param[k].c,start,bAppend);
 232     if (!bNB)
 233       append_interaction(il,type,nral,p->param[k].a);
 234   }
 235 }
 236
 237 static void new_interaction_list(t_ilist *ilist)
 238 {
 239   int i;
 240
 241   ilist->nr=0;
 242   for(i=0; (i<MAXPROC); i++)
 243     ilist->multinr[i]=0;
 244   ilist->iatoms=NULL;
 245 }
 246
 247 void convert_params(int atnr,t_params nbtypes[],
 248                     t_params plist[],t_idef *idef)
 249 {
 250   int    i,j,maxtypes;
 251   unsigned long  flags;
 252
 253   maxtypes=0;
 254
 255   idef->ntypes   = 0;
 256   idef->atnr     = atnr;
 257   idef->pid      = 0;
 258   idef->functype = NULL;
 259   idef->iparams  = NULL;
 260   for(i=0; (i<F_NRE); i++) {
 261     idef->il[i].nr=0;
 262     idef->il[i].iatoms=NULL;
 263   }
 264   enter_function(&(nbtypes[F_LJ]),  (t_functype)F_LJ,  idef,
 265                  &maxtypes,TRUE,TRUE);
 266   enter_function(&(nbtypes[F_BHAM]),(t_functype)F_BHAM,idef,
 267                  &maxtypes,TRUE,TRUE);
 268   enter_function(&(plist[F_POSRES]),(t_functype)F_POSRES,idef,
 269                  &maxtypes,FALSE,TRUE);
 270
 271   for(i=0; (i<F_NRE); i++) {
 272     flags = interaction_function[i].flags;
 273     if ((i != F_LJ) && (i != F_BHAM) && (i != F_POSRES) &&
 274         ((flags & IF_BOND) || (flags & IF_DUMMY) || (flags & IF_CONSTRAINT)))
 275       enter_function(&(plist[i]),(t_functype)i,idef,&maxtypes,FALSE,FALSE);
 276   }
 277   if (debug)
 278     fprintf(debug,"%s, line %d: There are %d functypes in idef\n",
 279             __FILE__,__LINE__,idef->ntypes);
 280   for(j=0; (j<F_NRE); j++) {
 281     for (i=0; (i<MAXPROC); i++)
 282       idef->il[j].multinr[i]=idef->il[j].nr;
 283
 284     if (idef->il[j].nr > 0)
 285       printf("# %10s:   %d\n",interaction_function[j].name,idef->il[j].nr);
 286   }
 287 }