src/contrib/calcfdev.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2006, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "typedefs.h"
  43 #include "main.h"
  44 #include "vec.h"
  45 #include "txtdump.h"
  46
  47 void calc_force(int natom,rvec f[],rvec fff[])
  48 {
  49   int  i,j,m;
  50   int  jindex[] = { 0, 5, 10};
  51   rvec dx,df;
  52   real msf1,msf2;
  53
  54   for(j=0; (j<2); j++) {
  55     clear_rvec(fff[j]);
  56     for(i=jindex[j]; (i<jindex[j+1]); i++) {
  57       for(m=0; (m<DIM); m++) {
  58         fff[j][m] += f[i][m];
  59       }
  60     }
  61   }
  62
  63   msf1 = iprod(fff[0],fff[0]);
  64   msf2 = iprod(fff[1],fff[1]);
  65   if (debug) {
  66     pr_rvecs(debug,0,"force",f,natom);
  67
  68     fprintf(debug,"FMOL:  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f  %10.3f\n",
  69             fff[0][XX],fff[0][YY],fff[0][ZZ],fff[1][XX],fff[1][YY],fff[1][ZZ]);
  70     fprintf(debug,"RMSF:  %10.3e  %10.3e\n",msf1,msf2);
  71   }
  72 }
  73
  74
  75 void calc_f_dev(int natoms,real charge[],rvec x[],rvec f[],
  76                 t_idef *idef,real *xiH,real *xiS)
  77 {
  78   enum { wwwO, wwwH1, wwwH2, wwwS, wwwNR };
  79   int  lj_index[wwwNR] = { 0,    4,     4,     8 };
  80   real rmsf,dFH,dFS,dFO,dr_14;
  81   real q[wwwNR],c6[wwwNR],c12[wwwNR],c12ratio;
  82   rvec fff[2],dx;
  83   int  i,j,aj;
  84
  85   for(i=0; (i<wwwNR); i++) {
  86     q[i]   = charge[i];
  87     c12[i] = idef->iparams[lj_index[i]].lj.c12;
  88     c6[i]  = idef->iparams[lj_index[i]].lj.c6;
  89   }
  90
  91   calc_force(natoms,f,fff);
  92   rmsf = norm(fff[0]);
  93   dFS = 0;
  94   dFH = 0;
  95   dFO = 0;
  96   for(i=0; (i<4); i++) {
  97     for(j=4; (j<8); j++) {
  98       if (c12[i] != 0) {
  99         rvec_sub(x[i],x[j],dx);
 100         aj       = j % 4;
 101         dr_14    = pow(iprod(dx,dx),-7);
 102         c12ratio = -12*sqrt(c12[aj]/c12[i])*dr_14*iprod(fff[0],dx);
 103
 104         switch (i) {
 105         case wwwH1:
 106         case wwwH2:
 107           dFH += c12ratio;
 108           break;
 109         case wwwS:
 110           dFS += c12ratio;
 111           break;
 112         case wwwO:
 113           dFS += c12ratio;
 114           break;
 115         }
 116       }
 117     }
 118   }
 119
 120   if (debug) {
 121     fprintf(debug,"FFF: dFS=%10.3e,  dFH=%10.3e,  dFO=%10.3e, rmsf=%10.3e\n",
 122             dFS,dFH,dFO,rmsf);
 123   }
 124   if (dFH == 0)
 125     *xiH = 1;
 126   else
 127     *xiH=rmsf/(10*c12[wwwH1]*dFH);
 128
 129   if (dFS == 0)
 130     *xiS = 1;
 131   else
 132     *xiS=rmsf/(10*c12[wwwS]*dFS);
 133 }