src/mdlib/gmx_qhop_db.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include "typedefs.h"
   5 #include "smalloc.h"
   6 #include "string2.h"
   7 #include "futil.h"
   8 #include "symtab.h"
   9 #include "hackblock.h"
  10 #include "resall.h"
  11 #include "gmx_qhop_parm.h"
  12 #include "gmx_qhop_xml.h"
  13 #include "gmx_qhop_db.h"
  14
  15 #if 0
  16 /* These routines link to invalid dependencies in kernel/resall.c! */
  17
  18 typedef struct {
  19   int id,charge,natom;
  20   int ndonor,*donor,nacceptor,*acceptor;
  21 } qhop_resinfo_t;
  22
  23 typedef struct {
  24   int               nrtp;
  25   t_restp           *rtp;
  26   qhop_resinfo_t    *resinfo;
  27   int               bts[4];
  28   int               nrexcl;
  29   bool              bAllDih,bHH14,bRemoveDih;
  30   t_atomtype        atype;
  31   t_symtab          tab;
  32   int               ngqh;
  33   gmx_qhop          *gqh;
  34   t_qhop_parameters *qhop_param;
  35 } gmx_qhop_db_t;
  36
  37 static void fill_resinfo(t_restp *rtp,qhop_resinfo_t *ri)
  38 {
  39   int    j,k,m;
  40   bool   bDonor;
  41   double qtot;
  42
  43   qtot = 0;
  44   ri->natom = rtp->natom;
  45   ri->ndonor = 0;
  46   snew(ri->donor,ri->natom);
  47   ri->nacceptor = 0;
  48   snew(ri->acceptor,ri->natom);
  49   for(j=0; (j<rtp->natom); j++) {
  50     qtot += rtp->atom[j].q;
  51     /* First check atoms usually involved in Hbonds */
  52     if ((rtp->atom[j].atomnumber == 7) ||
  53         (rtp->atom[j].atomnumber == 8) ||
  54         (rtp->atom[j].atomnumber == 9) ||
  55         (rtp->atom[j].atomnumber == 16)) {
  56       ri->acceptor[ri->nacceptor++] = j;
  57       /* Now check whether this atom has a proton, i.e. is a donor */
  58       bDonor = FALSE;
  59       for(k=0; (k<rtp->rb[ebtsBONDS].nb); k++) {
  60         if (strcasecmp(*(rtp->atomname[j]),rtp->rb[ebtsBONDS].b[k].a[0]) == 0)
  61           for(m=0; (m<rtp->natom); m++)
  62             if ((strcasecmp(*(rtp->atomname[m]),rtp->rb[ebtsBONDS].b[k].a[1]) == 0) &&
  63                 (rtp->atom[m].atomnumber == 1) &&
  64                 (rtp->atom[m].q > 0))
  65               bDonor = TRUE;
  66       }
  67       if (bDonor)
  68         ri->donor[ri->ndonor++] = j;
  69     }
  70   }
  71   ri->charge = qtot;
  72 }
  73
  74 static void lo_fill_qhp(gmx_qhop gqh,char *name,real *xx)
  75 {
  76   double x;
  77   int    test;
  78
  79   test = gmx_qhop_get_value(gqh,name,&x);
  80   if (test == 1)
  81     *xx = x;
  82   else {
  83     fprintf(stderr,"WARNING: Could not extract %s from qhop data for %s-%s\n",
  84             name,gmx_qhop_get_donor(gqh),gmx_qhop_get_acceptor(gqh));
  85     *xx = 0;
  86   }
  87 }
  88
  89 static void fill_qhp(t_qhop_parameters *qhp,gmx_qhop gqh)
  90 {
  91   lo_fill_qhp(gqh,"alpha",&qhp->alpha);
  92   lo_fill_qhp(gqh,"beta",&qhp->beta);
  93   lo_fill_qhp(gqh,"gamma",&qhp->gamma);
  94   lo_fill_qhp(gqh,"k_1",&qhp->k_1);
  95   lo_fill_qhp(gqh,"k_2",&qhp->k_2);
  96   lo_fill_qhp(gqh,"k_3",&qhp->k_3);
  97   lo_fill_qhp(gqh,"m_1",&qhp->m_1);
  98   lo_fill_qhp(gqh,"m_2",&qhp->m_2);
  99   lo_fill_qhp(gqh,"m_3",&qhp->m_3);
 100   lo_fill_qhp(gqh,"s_A",&qhp->s_A);
 101   lo_fill_qhp(gqh,"t_A",&qhp->t_A);
 102   lo_fill_qhp(gqh,"v_A",&qhp->v_A);
 103   lo_fill_qhp(gqh,"s_B",&qhp->s_B);
 104   lo_fill_qhp(gqh,"s_C",&qhp->s_C);
 105   lo_fill_qhp(gqh,"t_C",&qhp->t_C);
 106   lo_fill_qhp(gqh,"v_C",&qhp->v_C);
 107   lo_fill_qhp(gqh,"f",&qhp->f);
 108   lo_fill_qhp(gqh,"g",&qhp->g);
 109   lo_fill_qhp(gqh,"h",&qhp->h);
 110   lo_fill_qhp(gqh,"p_1",&qhp->p_1);
 111   lo_fill_qhp(gqh,"q_1",&qhp->q_1);
 112   lo_fill_qhp(gqh,"q_2",&qhp->q_2);
 113   lo_fill_qhp(gqh,"q_3",&qhp->q_3);
 114   lo_fill_qhp(gqh,"r_1",&qhp->r_1);
 115   lo_fill_qhp(gqh,"r_2",&qhp->r_2);
 116   lo_fill_qhp(gqh,"r_3",&qhp->r_3);
 117 }
 118
 119 static void dump_qhp(FILE *fp,t_qhop_parameters *qhp)
 120 {
 121   fprintf(fp,"alpha = %g\n",qhp->alpha);
 122   fprintf(fp,"beta  = %g\n",qhp->beta);
 123   fprintf(fp,"gamma = %g\n",qhp->gamma);
 124   fprintf(fp,"k_1   = %g\n",qhp->k_1);
 125   fprintf(fp,"k_2   = %g\n",qhp->k_2);
 126   fprintf(fp,"k_3   = %g\n",qhp->k_3);
 127   fprintf(fp,"m_1   = %g\n",qhp->m_1);
 128   fprintf(fp,"m_2   = %g\n",qhp->m_2);
 129   fprintf(fp,"m_3   = %g\n",qhp->m_3);
 130   fprintf(fp,"s_A   = %g\n",qhp->s_A);
 131   fprintf(fp,"t_A   = %g\n",qhp->t_A);
 132   fprintf(fp,"v_A   = %g\n",qhp->v_A);
 133   fprintf(fp,"s_B   = %g\n",qhp->s_B);
 134   fprintf(fp,"s_C   = %g\n",qhp->s_C);
 135   fprintf(fp,"t_C   = %g\n",qhp->t_C);
 136   fprintf(fp,"v_C   = %g\n",qhp->v_C);
 137   fprintf(fp,"f     = %g\n",qhp->f);
 138   fprintf(fp,"g     = %g\n",qhp->g);
 139   fprintf(fp,"h     = %g\n",qhp->h);
 140   fprintf(fp,"p_1   = %g\n",qhp->p_1);
 141   fprintf(fp,"q_1   = %g\n",qhp->q_1);
 142   fprintf(fp,"q_2   = %g\n",qhp->q_2);
 143   fprintf(fp,"q_3   = %g\n",qhp->q_3);
 144   fprintf(fp,"r_1   = %g\n",qhp->r_1);
 145   fprintf(fp,"r_2   = %g\n",qhp->r_2);
 146   fprintf(fp,"r_3   = %g\n",qhp->r_3);
 147 }
 148
 149 gmx_qhop_db gmx_qhop_db_read(char *forcefield)
 150 {
 151   gmx_qhop_db_t *qdb;
 152   char buf[256];
 153   char *fn;
 154   int i,j;
 155   double qtot;
 156
 157   snew(qdb,1);
 158   open_symtab(&(qdb->tab));
 159   sprintf(buf,"%s-qhop",forcefield);
 160   qdb->atype = read_atype(forcefield,&(qdb->tab));
 161   qdb->nrtp = read_resall(buf,qdb->bts,&(qdb->rtp),qdb->atype,
 162                           &(qdb->tab),&(qdb->bAllDih),
 163                           &(qdb->nrexcl),&(qdb->bHH14),&(qdb->bRemoveDih));
 164   snew(qdb->resinfo,qdb->nrtp);
 165   for(i=0; (i<qdb->nrtp); i++)
 166     fill_resinfo(&qdb->rtp[i],&qdb->resinfo[i]);
 167
 168   sprintf(buf,"%s-qhop.dat",forcefield);
 169   fn = (char *)libfn(buf);
 170   /* Read the xml data file */
 171   qdb->gqh = gmx_qhops_read(fn,&qdb->ngqh);
 172   sprintf(buf,"%s-qhop-debug.dat",forcefield);
 173   gmx_qhops_write(buf,qdb->ngqh,qdb->gqh);
 174   snew(qdb->qhop_param,qdb->ngqh);
 175   for(i=0; (i<qdb->ngqh); i++) {
 176     fill_qhp(&(qdb->qhop_param[i]),qdb->gqh[i]);
 177     if (debug) {
 178       fprintf(debug,"Donor: %s  Acceptor: %s\n",
 179               gmx_qhop_get_donor(qdb->gqh[i]),
 180               gmx_qhop_get_acceptor(qdb->gqh[i]));
 181       dump_qhp(debug,&(qdb->qhop_param[i]));
 182     }
 183   }
 184   sfree(fn);
 185   return (gmx_qhop_db) qdb;
 186 }
 187
 188 int gmx_qhop_db_write(char *fn,gmx_qhop_db qdb)
 189 {
 190   gmx_qhop_db_t *db = (gmx_qhop_db_t *) qdb;
 191   FILE *fp;
 192
 193   fp=ffopen(fn,"w");
 194   print_resall(fp,db->bts,db->nrtp,db->rtp,db->atype,db->bAllDih,
 195                db->nrexcl,db->bHH14,db->bRemoveDih);
 196   ffclose(fp);
 197
 198   return 1;
 199 }
 200
 201 int gmx_qhop_db_done(gmx_qhop_db qdb)
 202 {
 203   fprintf(stderr,"gmx_qhop_db_done not implemented yet.\n");
 204
 205   return 1;
 206 }
 207
 208 int gmx_qhop_db_get_nstates(gmx_qhop_db qdb,char *resname)
 209 {
 210   gmx_qhop_db_t *db = (gmx_qhop_db_t *) qdb;
 211   int i,nstate=0;
 212
 213   for(i=0; (i<db->nrtp); i++) {
 214     if (strncmp(db->rtp[i].resname,resname,3) == 0) {
 215       nstate++;
 216     }
 217   }
 218   return nstate;
 219 }
 220
 221 int gmx_qhop_db_get_qstate(gmx_qhop_db qdb,char *resname,int istate)
 222 {
 223   gmx_qhop_db_t *db = (gmx_qhop_db_t *) qdb;
 224   int i,nstate=0;
 225
 226   for(i=0; (i<db->nrtp); i++) {
 227     if (strncmp(db->rtp[i].resname,resname,3) == 0) {
 228       nstate++;
 229       if (nstate == istate)
 230         return db->resinfo[i].charge;
 231     }
 232   }
 233   return nstate;
 234 }
 235
 236 char **gmx_qhop_db_get_donors(gmx_qhop_db qdb,char *resname,int state)
 237 {
 238   return NULL;
 239 }
 240
 241 char **gmx_qhop_db_get_acceptors(gmx_qhop_db qdb,char *resname,int state)
 242 {
 243   return NULL;
 244 }
 245
 246 int gmx_qhop_db_set_charges(gmx_qhop_db qdb,char *resname,int state,
 247                             int natoms,real q[])
 248 {
 249   return 0;
 250 }
 251
 252 int gmx_qhop_db_get_parameters(gmx_qhop_db qdb,
 253                                char *donor,char *acceptor,
 254                                t_qhop_parameters *qp)
 255 {
 256   gmx_qhop_db_t *db = (gmx_qhop_db_t *) qdb;
 257   char *aa,*dd;
 258   int i;
 259
 260   for(i=0; (i<db->ngqh); i++) {
 261     aa = gmx_qhop_get_acceptor(db->gqh[i]);
 262     dd = gmx_qhop_get_donor(db->gqh[i]);
 263     if (strncmp(donor,dd,3) == 0) {
 264       if (strncmp(acceptor,aa,3) == 0) {
 265         memcpy(qp,&(db->qhop_param[i]),sizeof(*qp));
 266
 267         return 1;
 268       }
 269     }
 270   }
 271   return 0;
 272 }
 273
 274 #else
 275 int gmx_qhop_db_dummy=0;
 276 #endif