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