src/kernel/add_par.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 <string.h>
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43 #include "grompp.h"
  44 #include "toputil.h"
  45 #include "hackblock.h"
  46 #include "string2.h"
  47 #include "gmx_fatal.h"
  48
  49 static void clear_atom_list(int i0, atom_id a[])
  50 {
  51   int i;
  52
  53   for (i=i0; i < MAXATOMLIST; i++)
  54     a[i]=-1;
  55 }
  56
  57 static void clear_force_param(int i0, real c[])
  58 {
  59   int i;
  60
  61   for (i=i0; i < MAXFORCEPARAM; i++)
  62     c[i]=NOTSET;
  63 }
  64
  65 void add_param(t_params *ps,int ai,int aj, real *c, char *s)
  66 {
  67   int i;
  68
  69   if ((ai < 0) || (aj < 0))
  70     gmx_fatal(FARGS,"Trying to add impossible atoms: ai=%d, aj=%d",ai,aj);
  71   pr_alloc(1,ps);
  72   ps->param[ps->nr].AI=ai;
  73   ps->param[ps->nr].AJ=aj;
  74   clear_atom_list(2, ps->param[ps->nr].a);
  75   if (c==NULL)
  76     clear_force_param(0, ps->param[ps->nr].c);
  77   else
  78     for(i=0; (i < MAXFORCEPARAM); i++)
  79       ps->param[ps->nr].c[i]=c[i];
  80   set_p_string(&(ps->param[ps->nr]),s);
  81   ps->nr++;
  82 }
  83
  84 void add_imp_param(t_params *ps,int ai,int aj,int ak,int al,real c0, real c1,
  85                    char *s)
  86 {
  87   pr_alloc(1,ps);
  88   ps->param[ps->nr].AI=ai;
  89   ps->param[ps->nr].AJ=aj;
  90   ps->param[ps->nr].AK=ak;
  91   ps->param[ps->nr].AL=al;
  92   clear_atom_list  (4, ps->param[ps->nr].a);
  93   ps->param[ps->nr].C0=c0;
  94   ps->param[ps->nr].C1=c1;
  95   clear_force_param(2, ps->param[ps->nr].c);
  96   set_p_string(&(ps->param[ps->nr]),s);
  97   ps->nr++;
  98 }
  99
 100 void add_dih_param(t_params *ps,int ai,int aj,int ak,int al,real c0, real c1,
 101                    real c2,char *s)
 102 {
 103   pr_alloc(1,ps);
 104   ps->param[ps->nr].AI=ai;
 105   ps->param[ps->nr].AJ=aj;
 106   ps->param[ps->nr].AK=ak;
 107   ps->param[ps->nr].AL=al;
 108   clear_atom_list  (4, ps->param[ps->nr].a);
 109   ps->param[ps->nr].C0=c0;
 110   ps->param[ps->nr].C1=c1;
 111   ps->param[ps->nr].C2=c2;
 112   clear_force_param(3, ps->param[ps->nr].c);
 113   set_p_string(&(ps->param[ps->nr]),s);
 114   ps->nr++;
 115 }
 116
 117 void add_cmap_param(t_params *ps, int ai, int aj, int ak, int al, int am, char *s)
 118 {
 119         pr_alloc(1,ps);
 120         ps->param[ps->nr].AI=ai;
 121         ps->param[ps->nr].AJ=aj;
 122         ps->param[ps->nr].AK=ak;
 123         ps->param[ps->nr].AL=al;
 124         ps->param[ps->nr].AM=am;
 125         clear_atom_list(5,ps->param[ps->nr].a);
 126         clear_force_param(0,ps->param[ps->nr].c);
 127         set_p_string(&(ps->param[ps->nr]),s);
 128         ps->nr++;
 129 }
 130
 131 void add_vsite2_atoms(t_params *ps,int ai,int aj,int ak)
 132 {
 133   pr_alloc(1,ps);
 134   ps->param[ps->nr].AI=ai;
 135   ps->param[ps->nr].AJ=aj;
 136   ps->param[ps->nr].AK=ak;
 137   clear_atom_list  (3, ps->param[ps->nr].a);
 138   clear_force_param(0, ps->param[ps->nr].c);
 139   set_p_string(&(ps->param[ps->nr]),"");
 140   ps->nr++;
 141 }
 142
 143 void add_vsite2_param(t_params *ps,int ai,int aj,int ak,real c0)
 144 {
 145   pr_alloc(1,ps);
 146   ps->param[ps->nr].AI=ai;
 147   ps->param[ps->nr].AJ=aj;
 148   ps->param[ps->nr].AK=ak;
 149   clear_atom_list  (3, ps->param[ps->nr].a);
 150   ps->param[ps->nr].C0=c0;
 151   clear_force_param(1, ps->param[ps->nr].c);
 152   set_p_string(&(ps->param[ps->nr]),"");
 153   ps->nr++;
 154 }
 155
 156 void add_vsite3_param(t_params *ps,int ai,int aj,int ak,int al,
 157                     real c0, real c1)
 158 {
 159   pr_alloc(1,ps);
 160   ps->param[ps->nr].AI=ai;
 161   ps->param[ps->nr].AJ=aj;
 162   ps->param[ps->nr].AK=ak;
 163   ps->param[ps->nr].AL=al;
 164   clear_atom_list  (4, ps->param[ps->nr].a);
 165   ps->param[ps->nr].C0=c0;
 166   ps->param[ps->nr].C1=c1;
 167   clear_force_param(2, ps->param[ps->nr].c);
 168   set_p_string(&(ps->param[ps->nr]),"");
 169   ps->nr++;
 170 }
 171
 172 void add_vsite3_atoms(t_params *ps,int ai,int aj,int ak,int al, bool bSwapParity)
 173 {
 174   pr_alloc(1,ps);
 175   ps->param[ps->nr].AI=ai;
 176   ps->param[ps->nr].AJ=aj;
 177   ps->param[ps->nr].AK=ak;
 178   ps->param[ps->nr].AL=al;
 179   clear_atom_list  (4, ps->param[ps->nr].a);
 180   clear_force_param(0, ps->param[ps->nr].c);
 181   if (bSwapParity)
 182     ps->param[ps->nr].C1=-1;
 183   set_p_string(&(ps->param[ps->nr]),"");
 184   ps->nr++;
 185 }
 186
 187 void add_vsite4_atoms(t_params *ps,int ai,int aj,int ak,int al,int am)
 188 {
 189   pr_alloc(1,ps);
 190   ps->param[ps->nr].AI=ai;
 191   ps->param[ps->nr].AJ=aj;
 192   ps->param[ps->nr].AK=ak;
 193   ps->param[ps->nr].AL=al;
 194   ps->param[ps->nr].AM=am;
 195   clear_atom_list  (5, ps->param[ps->nr].a);
 196   clear_force_param(0, ps->param[ps->nr].c);
 197   set_p_string(&(ps->param[ps->nr]),"");
 198   ps->nr++;
 199 }
 200
 201 int search_jtype(t_restp *rtp,char *name,bool bNterm)
 202 {
 203   int  niter,iter,j,k,kmax,jmax,minstrlen;
 204   char *rtpname,searchname[12];
 205
 206   strcpy(searchname,name);
 207
 208   /* Do a best match comparison */
 209   /* for protein N-terminus, allow renaming of H1, H2 and H3 to H */
 210   if ( bNterm && (strlen(searchname)==2) && (searchname[0] == 'H') &&
 211        ( (searchname[1] == '1') || (searchname[1] == '2') ||
 212          (searchname[1] == '3') ) ) {
 213     niter = 2;
 214   } else {
 215     niter = 1;
 216   }
 217   kmax=0;
 218   jmax=-1;
 219   for(iter=0; (iter<niter && jmax==-1); iter++) {
 220     if (iter == 1) {
 221       /* Try without the hydrogen number in the N-terminus */
 222       searchname[1] = '\0';
 223     }
 224     for(j=0; (j<rtp->natom); j++) {
 225       rtpname=*(rtp->atomname[j]);
 226       if (gmx_strcasecmp(searchname,rtpname) == 0) {
 227         jmax=j;
 228         kmax=strlen(searchname);
 229         break;
 230       }
 231       if (iter == niter - 1) {
 232         minstrlen = min(strlen(searchname), strlen(rtpname));
 233         for(k=0; k < minstrlen; k++)
 234           if (searchname[k] != rtpname[k])
 235             break;
 236         if (k > kmax) {
 237           kmax=k;
 238           jmax=j;
 239         }
 240       }
 241     }
 242   }
 243   if (jmax == -1)
 244     gmx_fatal(FARGS,"Atom %s not found in rtp database in residue %s",
 245               searchname,rtp->resname);
 246   if (kmax != strlen(searchname))
 247     gmx_fatal(FARGS,"Atom %s not found in rtp database in residue %s, "
 248               "it looks a bit like %s",
 249               searchname,rtp->resname,*(rtp->atomname[jmax]));
 250   return jmax;
 251 }
 252