src/contrib/calcfdev.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_calcfdev_c = "$Id$";
  30
  31 #include "typedefs.h"
  32 #include "main.h"
  33 #include "vec.h"
  34 #include "txtdump.h"
  35
  36 void calc_force(int natom,rvec f[],rvec fff[])
  37 {
  38   int  i,j,m;
  39   int  jindex[] = { 0, 5, 10};
  40   rvec dx,df;
  41   real msf1,msf2;
  42
  43   for(j=0; (j<2); j++) {
  44     clear_rvec(fff[j]);
  45     for(i=jindex[j]; (i<jindex[j+1]); i++) {
  46       for(m=0; (m<DIM); m++) {
  47         fff[j][m] += f[i][m];
  48       }
  49     }
  50   }
  51
  52   msf1 = iprod(fff[0],fff[0]);
  53   msf2 = iprod(fff[1],fff[1]);
  54   if (debug) {
  55     pr_rvecs(debug,0,"force",f,natom);
  56
  57     fprintf(debug,"FMOL:  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f\n",
  58             fff[0][XX],fff[0][YY],fff[0][ZZ],fff[1][XX],fff[1][YY],fff[1][ZZ]);
  59     fprintf(debug,"RMSF:  %10.3e  %10.3e\n",msf1,msf2);
  60   }
  61 }
  62
  63
  64 void calc_f_dev(int natoms,real charge[],rvec x[],rvec f[],
  65                 t_idef *idef,real *xiH,real *xiS)
  66 {
  67   enum { wwwO, wwwH1, wwwH2, wwwS, wwwNR };
  68   int  lj_index[wwwNR] = { 0,    4,     4,     8 };
  69   real rmsf,dFH,dFS,dFO,dr_14;
  70   real q[wwwNR],c6[wwwNR],c12[wwwNR],c12ratio;
  71   rvec fff[2],dx;
  72   int  i,j,aj;
  73
  74   for(i=0; (i<wwwNR); i++) {
  75     q[i]   = charge[i];
  76     c12[i] = idef->iparams[lj_index[i]].lj.c12;
  77     c6[i]  = idef->iparams[lj_index[i]].lj.c6;
  78   }
  79
  80   calc_force(natoms,f,fff);
  81   rmsf = norm(fff[0]);
  82   dFS = 0;
  83   dFH = 0;
  84   dFO = 0;
  85   for(i=0; (i<4); i++) {
  86     for(j=4; (j<8); j++) {
  87       if (c12[i] != 0) {
  88         rvec_sub(x[i],x[j],dx);
  89         aj       = j % 4;
  90         dr_14    = pow(iprod(dx,dx),-7);
  91         c12ratio = -12*sqrt(c12[aj]/c12[i])*dr_14*iprod(fff[0],dx);
  92
  93         switch (i) {
  94         case wwwH1:
  95         case wwwH2:
  96           dFH += c12ratio;
  97           break;
  98         case wwwS:
  99           dFS += c12ratio;
 100           break;
 101         case wwwO:
 102           dFS += c12ratio;
 103           break;
 104         }
 105       }
 106     }
 107   }
 108
 109   if (debug) {
 110     fprintf(debug,"FFF: dFS=%10.3e,  dFH=%10.3e,  dFO=%10.3e, rmsf=%10.3e\n",
 111             dFS,dFH,dFO,rmsf);
 112   }
 113   if (dFH == 0)
 114     *xiH = 1;
 115   else
 116     *xiH=rmsf/(10*c12[wwwH1]*dFH);
 117
 118   if (dFS == 0)
 119     *xiS = 1;
 120   else
 121     *xiS=rmsf/(10*c12[wwwS]*dFS);
 122 }