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 3.0
  11  *
  12  * Copyright (c) 1991-2001
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
  32  *
  33  * And Hey:
  34  * GROup of MAchos and Cynical Suckers
  35  */
  36 static char *SRCID_calcfdev_c = "$Id$";
  37 #include "typedefs.h"
  38 #include "main.h"
  39 #include "vec.h"
  40 #include "txtdump.h"
  41
  42 void calc_force(int natom,rvec f[],rvec fff[])
  43 {
  44   int  i,j,m;
  45   int  jindex[] = { 0, 5, 10};
  46   rvec dx,df;
  47   real msf1,msf2;
  48
  49   for(j=0; (j<2); j++) {
  50     clear_rvec(fff[j]);
  51     for(i=jindex[j]; (i<jindex[j+1]); i++) {
  52       for(m=0; (m<DIM); m++) {
  53         fff[j][m] += f[i][m];
  54       }
  55     }
  56   }
  57
  58   msf1 = iprod(fff[0],fff[0]);
  59   msf2 = iprod(fff[1],fff[1]);
  60   if (debug) {
  61     pr_rvecs(debug,0,"force",f,natom);
  62
  63     fprintf(debug,"FMOL:  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f\n",
  64             fff[0][XX],fff[0][YY],fff[0][ZZ],fff[1][XX],fff[1][YY],fff[1][ZZ]);
  65     fprintf(debug,"RMSF:  %10.3e  %10.3e\n",msf1,msf2);
  66   }
  67 }
  68
  69
  70 void calc_f_dev(int natoms,real charge[],rvec x[],rvec f[],
  71                 t_idef *idef,real *xiH,real *xiS)
  72 {
  73   enum { wwwO, wwwH1, wwwH2, wwwS, wwwNR };
  74   int  lj_index[wwwNR] = { 0,    4,     4,     8 };
  75   real rmsf,dFH,dFS,dFO,dr_14;
  76   real q[wwwNR],c6[wwwNR],c12[wwwNR],c12ratio;
  77   rvec fff[2],dx;
  78   int  i,j,aj;
  79
  80   for(i=0; (i<wwwNR); i++) {
  81     q[i]   = charge[i];
  82     c12[i] = idef->iparams[lj_index[i]].lj.c12;
  83     c6[i]  = idef->iparams[lj_index[i]].lj.c6;
  84   }
  85
  86   calc_force(natoms,f,fff);
  87   rmsf = norm(fff[0]);
  88   dFS = 0;
  89   dFH = 0;
  90   dFO = 0;
  91   for(i=0; (i<4); i++) {
  92     for(j=4; (j<8); j++) {
  93       if (c12[i] != 0) {
  94         rvec_sub(x[i],x[j],dx);
  95         aj       = j % 4;
  96         dr_14    = pow(iprod(dx,dx),-7);
  97         c12ratio = -12*sqrt(c12[aj]/c12[i])*dr_14*iprod(fff[0],dx);
  98
  99         switch (i) {
 100         case wwwH1:
 101         case wwwH2:
 102           dFH += c12ratio;
 103           break;
 104         case wwwS:
 105           dFS += c12ratio;
 106           break;
 107         case wwwO:
 108           dFS += c12ratio;
 109           break;
 110         }
 111       }
 112     }
 113   }
 114
 115   if (debug) {
 116     fprintf(debug,"FFF: dFS=%10.3e,  dFH=%10.3e,  dFO=%10.3e, rmsf=%10.3e\n",
 117             dFS,dFH,dFO,rmsf);
 118   }
 119   if (dFH == 0)
 120     *xiH = 1;
 121   else
 122     *xiH=rmsf/(10*c12[wwwH1]*dFH);
 123
 124   if (dFS == 0)
 125     *xiS = 1;
 126   else
 127     *xiS=rmsf/(10*c12[wwwS]*dFS);
 128 }