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