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 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  * GROningen Mixture of Alchemy and Childrens' Stories
  31  */
  32 static char *SRCID_convparm_c = "$Id$";
  33 #include <math.h>
  34 #include "sysstuff.h"
  35 #include "physics.h"
  36 #include "vec.h"
  37 #include "assert.h"
  38 #include "smalloc.h"
  39 #include "typedefs.h"
  40 #include "fatal.h"
  41 #include "topio.h"
  42 #include "toputil.h"
  43 #include "convparm.h"
  44 #include "names.h"
  45
  46 static void assign_param(t_functype ftype,t_iparams *new,
  47                          real old[MAXFORCEPARAM])
  48 {
  49   int i,j;
  50
  51   /* Set to zero */
  52   for(j=0; (j<MAXFORCEPARAM); j++)
  53     new->generic.buf[j]=0.0;
  54
  55   switch (ftype) {
  56   case F_G96ANGLES:
  57     /* Post processing of input data: store cosine iso angle itself */
  58     new->harmonic.rA =cos(old[0]*DEG2RAD);
  59     new->harmonic.krA=old[1];
  60     new->harmonic.rB =cos(old[2]*DEG2RAD);
  61     new->harmonic.krB=old[3];
  62     break;
  63   case F_G96BONDS:
  64     /* Post processing of input data: store square of length itself */
  65     new->harmonic.rA =sqr(old[0]);
  66     new->harmonic.krA=old[1];
  67     new->harmonic.rB =sqr(old[2]);
  68     new->harmonic.krB=old[3];
  69     break;
  70   case F_ANGLES:
  71   case F_BONDS:
  72   case F_HARMONIC:
  73   case F_IDIHS:
  74     new->harmonic.rA =old[0];
  75     new->harmonic.krA=old[1];
  76     new->harmonic.rB =old[2];
  77     new->harmonic.krB=old[3];
  78     break;
  79   case F_MORSE:
  80     new->morse.b0    =old[0];
  81     new->morse.cb    =old[1];
  82     new->morse.beta  =old[2];
  83     break;
  84   case F_CUBICBONDS:
  85     new->cubic.b0    =old[0];
  86     new->cubic.kb    =old[1];
  87     new->cubic.kcub  =old[2];
  88     break;
  89   case F_CONNBONDS:
  90     break;
  91   case F_WPOL:
  92     new->wpol.kx     =old[0];
  93     new->wpol.ky     =old[1];
  94     new->wpol.kz     =old[2];
  95     new->wpol.rOH    =old[3];
  96     new->wpol.rHH    =old[4];
  97     new->wpol.rOD    =old[5];
  98     break;
  99   case F_BHAM:
 100     new->bham.a = old[0];
 101     new->bham.b = old[1];
 102     new->bham.c = old[2];
 103     break;
 104   case F_LJ14:
 105     new->lj14.c6A  = old[0];
 106     new->lj14.c12A = old[1];
 107     new->lj14.c6B  = old[2];
 108     new->lj14.c12B = old[3];
 109     break;
 110   case F_LJ:
 111     new->lj.c6  = old[0];
 112     new->lj.c12 = old[1];
 113     break;
 114   case F_PDIHS:
 115   case F_ANGRES:
 116   case F_ANGRESZ:
 117     new->pdihs.phiA=old[0];
 118     new->pdihs.cpA =old[1];
 119     new->pdihs.mult=old[2];
 120     new->pdihs.phiB=old[3];
 121     new->pdihs.cpB =old[4];
 122     break;
 123   case F_POSRES:
 124     new->posres.fc[XX]   = old[0];
 125     new->posres.fc[YY]   = old[1];
 126     new->posres.fc[ZZ]   = old[2];
 127     new->posres.pos0[XX] = old[3];
 128     new->posres.pos0[YY] = old[4];
 129     new->posres.pos0[ZZ] = old[5];
 130     break;
 131   case F_DISRES:
 132     new->disres.label = old[0];
 133     new->disres.type  = old[1];
 134     new->disres.low   = old[2];
 135     new->disres.up1   = old[3];
 136     new->disres.up2   = old[4];
 137     new->disres.kfac  = old[5];
 138     break;
 139   case F_ORIRES:
 140     if (old[0] < 0)
 141       fatal_error(0,"Found experiment number for orientation restraints which is smaller than 1 (%d)",old[0]);
 142     new->orires.ex    = old[0] - 1;
 143     new->orires.label = old[1];
 144     new->orires.power = old[2];
 145     new->orires.c     = old[3];
 146     new->orires.obs   = old[4];
 147     new->orires.kfac  = old[5];
 148     break;
 149   case F_RBDIHS:
 150     for (i=0; (i<NR_RBDIHS); i++) {
 151       new->rbdihs.rbcA[i]=old[i];
 152       new->rbdihs.rbcB[i]=old[NR_RBDIHS+i];
 153     }
 154     break;
 155   case F_SHAKE:
 156   case F_SHAKENC:
 157     new->shake.dA = old[0];
 158     new->shake.dB = old[1];
 159     break;
 160   case F_SETTLE:
 161     new->settle.doh=old[0];
 162     new->settle.dhh=old[1];
 163     break;
 164   case F_DUMMY2:
 165   case F_DUMMY3:
 166   case F_DUMMY3FD:
 167   case F_DUMMY3OUT:
 168   case F_DUMMY4FD:
 169     new->dummy.a=old[0];
 170     new->dummy.b=old[1];
 171     new->dummy.c=old[2];
 172     new->dummy.d=old[3];
 173     new->dummy.e=old[4];
 174     new->dummy.f=old[5];
 175     break;
 176   case F_DUMMY3FAD:
 177     new->dummy.a=old[1] * cos(DEG2RAD * old[0]);
 178     new->dummy.b=old[1] * sin(DEG2RAD * old[0]);
 179     new->dummy.c=old[2];
 180     new->dummy.d=old[3];
 181     new->dummy.e=old[4];
 182     new->dummy.f=old[5];
 183     break;
 184   default:
 185     fatal_error(0,"unknown function type %d in %s line %d",
 186                 ftype,__FILE__,__LINE__);
 187   }
 188 }
 189
 190 static int enter_params(t_idef *idef, t_functype ftype,
 191                         real forceparams[MAXFORCEPARAM],int start,bool bAppend)
 192 {
 193   t_iparams new;
 194   int       type;
 195
 196   assign_param(ftype,&new,forceparams);
 197   if (!bAppend) {
 198     for (type=start; (type<idef->ntypes); type++) {
 199       if (idef->functype[type]==ftype) {
 200         if (memcmp(&new,&idef->iparams[type],(size_t)sizeof(new)) == 0)
 201           return type;
 202       }
 203     }
 204   }
 205   else
 206     type=idef->ntypes;
 207   if (debug)
 208     fprintf(debug,"copying new to idef->iparams[%d] (ntypes=%d)\n",
 209             type,idef->ntypes);
 210   memcpy(&idef->iparams[type],&new,(size_t)sizeof(new));
 211
 212   idef->ntypes++;
 213   idef->functype[type]=ftype;
 214
 215   return type;
 216 }
 217
 218 static void append_interaction(t_ilist *ilist,
 219                                int type,int nral,atom_id a[MAXATOMLIST])
 220 {
 221   int i,where1;
 222
 223   where1     = ilist->nr;
 224   ilist->nr += nral+1;
 225
 226   ilist->iatoms[where1++]=type;
 227   for (i=0; (i<nral); i++)
 228     ilist->iatoms[where1++]=a[i];
 229 }
 230
 231 static void enter_function(t_params *p,t_functype ftype,
 232                            t_idef *idef,int *maxtypes,bool bNB,bool bAppend)
 233 {
 234   int     k,type,nr,nral,delta,start;
 235   t_ilist *il;
 236
 237   il    = &(idef->il[ftype]);
 238   start = idef->ntypes;
 239   nr    = p->nr;
 240   nral  = NRAL(ftype);
 241   delta = nr*(nral+1);
 242   srenew(il->iatoms,il->nr+delta);
 243
 244   for (k=0; k<nr; k++) {
 245     if (*maxtypes <= idef->ntypes) {
 246       *maxtypes += 1000;
 247       srenew(idef->functype,*maxtypes);
 248       srenew(idef->iparams, *maxtypes);
 249       if (debug)
 250         fprintf(debug,"%s, line %d: srenewed idef->functype and idef->iparams to %d\n",
 251                 __FILE__,__LINE__,*maxtypes);
 252     }
 253     type = enter_params(idef,ftype,p->param[k].c,start,bAppend);
 254     if (!bNB)
 255       append_interaction(il,type,nral,p->param[k].a);
 256   }
 257 }
 258
 259 static void new_interaction_list(t_ilist *ilist)
 260 {
 261   int i;
 262
 263   ilist->nr=0;
 264   for(i=0; (i<MAXNODES); i++)
 265     ilist->multinr[i]=0;
 266   ilist->iatoms=NULL;
 267 }
 268
 269 void convert_params(int atnr,t_params nbtypes[],
 270                     t_params plist[],t_idef *idef)
 271 {
 272   int    i,j,maxtypes;
 273   unsigned long  flags;
 274
 275   maxtypes=0;
 276
 277   idef->ntypes   = 0;
 278   idef->atnr     = atnr;
 279   idef->nodeid   = 0;
 280   idef->functype = NULL;
 281   idef->iparams  = NULL;
 282   for(i=0; (i<F_NRE); i++) {
 283     idef->il[i].nr=0;
 284     idef->il[i].iatoms=NULL;
 285   }
 286   enter_function(&(nbtypes[F_LJ]),  (t_functype)F_LJ,  idef,
 287                  &maxtypes,TRUE,TRUE);
 288   enter_function(&(nbtypes[F_BHAM]),(t_functype)F_BHAM,idef,
 289                  &maxtypes,TRUE,TRUE);
 290   enter_function(&(plist[F_POSRES]),(t_functype)F_POSRES,idef,
 291                  &maxtypes,FALSE,TRUE);
 292
 293   for(i=0; (i<F_NRE); i++) {
 294     flags = interaction_function[i].flags;
 295     if ((i != F_LJ) && (i != F_BHAM) && (i != F_POSRES) &&
 296         ((flags & IF_BOND) || (flags & IF_DUMMY) || (flags & IF_CONSTRAINT)))
 297       enter_function(&(plist[i]),(t_functype)i,idef,&maxtypes,FALSE,FALSE);
 298   }
 299   if (debug)
 300     fprintf(debug,"%s, line %d: There are %d functypes in idef\n",
 301             __FILE__,__LINE__,idef->ntypes);
 302   for(j=0; (j<F_NRE); j++) {
 303     for (i=0; (i<MAXNODES); i++)
 304       idef->il[j].multinr[i]=idef->il[j].nr;
 305
 306     if (idef->il[j].nr > 0)
 307       printf("# %10s:   %d\n",interaction_function[j].name,idef->il[j].nr);
 308   }
 309 }