src/local/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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_glaasje_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "smalloc.h"
  35 #include "glaasje.h"
  36 #include "macros.h"
  37
  38 void do_glas(FILE *log,int start,int homenr,rvec x[],rvec f[],
  39              t_forcerec *fr,t_mdatoms *md,int atnr,t_inputrec *ir,
  40              real ener[])
  41 {
  42   static bool   bFirst=TRUE,bGlas;
  43   static real   d[2],pi6,pi12,rc9,rc4,rc10,rc3,rc;
  44   static real   *c6,*c12;
  45   real wd,wdd,zi,fz,dd,d10,d4,d9,d3,r9,r3,sign,cc6,cc12;
  46   int  *type;
  47   int  i,j,ti;
  48
  49   type=md->typeA;
  50   if (bFirst) {
  51     pi6  = ir->userreal1;
  52     pi12 = ir->userreal2;
  53     d[0] = ir->userreal3;
  54     d[1] = ir->userreal4;
  55
  56     /* Check whether these constants have been set. */
  57     bGlas = (pi6 != 0) && (pi12 != 0) && (d[0] != 0) && (d[1] != 0);
  58
  59     if (bGlas) {
  60       if (ir->bDispCorr) {
  61         fatal_error(0,"Can not have Long Range C6 corrections and GLASMD");
  62       }
  63       rc   = max(fr->rvdw,fr->rlist);
  64       rc3  = rc*rc*rc;
  65       rc4  = rc3*rc;
  66       rc9  = rc3*rc3*rc3;
  67       rc10 = rc9*rc;
  68
  69       fprintf(log,
  70               "Constants for GLASMD: pi6 = %10g, pi12 = %10g\n"
  71               "                      d1  = %10g, d2   = %10g\n"
  72               "                      rc3 = %10g, rc4  = %10g\n"
  73               "                      rc9 = %10g, rc10 = %10g\n",
  74               pi6,pi12,d[0],d[1],rc3,rc4,rc9,rc10);
  75       if (d[0] > d[1])
  76         fatal_error(0,"d1 > d2 for GLASMD (check log file)");
  77
  78       snew(c6,atnr);
  79       snew(c12,atnr);
  80
  81       for(i=0; (i<atnr); i++) {
  82         c6[i]  = C6 (fr->nbfp,atnr,i,i);
  83         c12[i] = C12(fr->nbfp,atnr,i,i);
  84       }
  85     }
  86     else
  87       fprintf(stderr,"No glasmd!\n");
  88     bFirst = FALSE;
  89   }
  90
  91   if (bGlas) {
  92     wd=0;
  93     for(i=start; (i<start+homenr); i++) {
  94       ti   = type[i];
  95       if ((c6[ti] != 0) || (c12[ti] != 0)) {
  96         zi   = x[i][ZZ];
  97         cc6  = M_PI*sqrt(c6[ti]*pi6);
  98         cc12 = M_PI*sqrt(c12[ti]*pi12);
  99
 100         /* Use a factor for the sign, this saves taking absolute values */
 101         sign = 1;
 102         for(j=0; (j<2); j++) {
 103           dd = sign*(zi-d[j]);
 104           if (dd >= rc) {
 105             d3  = dd*dd*dd;
 106             d9  = d3*d3*d3;
 107             wdd = cc12/(45.0*d9) - cc6/(6.0*d3);
 108             d4  = d3*dd;
 109             d10 = d9*dd;
 110             fz  = sign*(cc12/(5.0*d10) - cc6/(2.0*d4));
 111           }
 112           else {
 113             wdd = cc12*(2.0/(9.0*rc9) - dd/(5.0*rc10)) -
 114               cc6*(2.0/(3.0*rc3) - dd/(2.0*rc4));
 115             fz  = sign*(cc12/(5.0*rc10)-cc6/(2.0*rc4));
 116           }
 117           wd       += wdd;
 118           f[i][ZZ] += fz;
 119           sign      = -sign;
 120         }
 121       }
 122     }
 123     ener[F_LJLR] = wd;
 124   }
 125 }
 126