src/contrib/glaasje.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  * Gromacs Runs On Most of All Computer Systems
  35  */
  36 static char *SRCID_glaasje_c = "$Id$";
  37 #include <stdio.h>
  38 #include <math.h>
  39 #include "typedefs.h"
  40 #include "smalloc.h"
  41 #include "glaasje.h"
  42 #include "macros.h"
  43
  44 void do_glas(FILE *log,int start,int homenr,rvec x[],rvec f[],
  45              t_forcerec *fr,t_mdatoms *md,int atnr,t_inputrec *ir,
  46              real ener[])
  47 {
  48   static bool   bFirst=TRUE,bGlas;
  49   static real   d[2],pi6,pi12,rc9,rc4,rc10,rc3,rc;
  50   static real   *c6,*c12;
  51   real wd,wdd,zi,fz,dd,d10,d4,d9,d3,r9,r3,sign,cc6,cc12;
  52   int  *type;
  53   int  i,j,ti;
  54
  55   type=md->typeA;
  56   if (bFirst) {
  57     pi6  = ir->userreal1;
  58     pi12 = ir->userreal2;
  59     d[0] = ir->userreal3;
  60     d[1] = ir->userreal4;
  61
  62     /* Check whether these constants have been set. */
  63     bGlas = (pi6 != 0) && (pi12 != 0) && (d[0] != 0) && (d[1] != 0);
  64
  65     if (bGlas) {
  66       if (ir->eDispCorr != edispcNO) {
  67         fatal_error(0,"Can not have Long Range C6 corrections and GLASMD");
  68       }
  69       rc   = max(fr->rvdw,fr->rlist);
  70       rc3  = rc*rc*rc;
  71       rc4  = rc3*rc;
  72       rc9  = rc3*rc3*rc3;
  73       rc10 = rc9*rc;
  74
  75       fprintf(log,
  76               "Constants for GLASMD: pi6 = %10g, pi12 = %10g\n"
  77               "                      d1  = %10g, d2   = %10g\n"
  78               "                      rc3 = %10g, rc4  = %10g\n"
  79               "                      rc9 = %10g, rc10 = %10g\n",
  80               pi6,pi12,d[0],d[1],rc3,rc4,rc9,rc10);
  81       if (d[0] > d[1])
  82         fatal_error(0,"d1 > d2 for GLASMD (check log file)");
  83
  84       snew(c6,atnr);
  85       snew(c12,atnr);
  86
  87       for(i=0; (i<atnr); i++) {
  88         c6[i]  = C6 (fr->nbfp,atnr,i,i);
  89         c12[i] = C12(fr->nbfp,atnr,i,i);
  90       }
  91     }
  92     else
  93       fprintf(stderr,"No glasmd!\n");
  94     bFirst = FALSE;
  95   }
  96
  97   if (bGlas) {
  98     wd=0;
  99     for(i=start; (i<start+homenr); i++) {
 100       ti   = type[i];
 101       if ((c6[ti] != 0) || (c12[ti] != 0)) {
 102         zi   = x[i][ZZ];
 103         cc6  = M_PI*sqrt(c6[ti]*pi6);
 104         cc12 = M_PI*sqrt(c12[ti]*pi12);
 105
 106         /* Use a factor for the sign, this saves taking absolute values */
 107         sign = 1;
 108         for(j=0; (j<2); j++) {
 109           dd = sign*(zi-d[j]);
 110           if (dd >= rc) {
 111             d3  = dd*dd*dd;
 112             d9  = d3*d3*d3;
 113             wdd = cc12/(45.0*d9) - cc6/(6.0*d3);
 114             d4  = d3*dd;
 115             d10 = d9*dd;
 116             fz  = sign*(cc12/(5.0*d10) - cc6/(2.0*d4));
 117           }
 118           else {
 119             wdd = cc12*(2.0/(9.0*rc9) - dd/(5.0*rc10)) -
 120               cc6*(2.0/(3.0*rc3) - dd/(2.0*rc4));
 121             fz  = sign*(cc12/(5.0*rc10)-cc6/(2.0*rc4));
 122           }
 123           wd       += wdd;
 124           f[i][ZZ] += fz;
 125           sign      = -sign;
 126         }
 127       }
 128     }
 129     ener[F_LJLR] = wd;
 130   }
 131 }
 132