src/kernel/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.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * GROningen Mixture of Alchemy and Childrens' Stories
  31  */
  32
  33 #include <stdio.h>
  34 #include <math.h>
  35 #include "typedefs.h"
  36 #include "smalloc.h"
  37 #include "maths.h"
  38 #include "glaasje.h"
  39 #include "macros.h"
  40
  41 void do_glas(FILE *log,int start,int homenr,rvec x[],rvec f[],
  42              t_forcerec *fr,t_mdatoms *md,int atnr,t_inputrec *ir,
  43              real ener[])
  44 {
  45   static bool   bFirst=TRUE,bGlas;
  46   static real   d[2],pi6,pi12,rc9,rc4,rc10,rc3,rc;
  47   static real   *c6,*c12;
  48   real wd,wdd,zi,fz,dd,d10,d4,d9,d3,r9,r3,sign,cc6,cc12;
  49   int  *type;
  50   int  i,j,ti;
  51
  52   type=md->typeA;
  53   if (bFirst) {
  54     pi6  = ir->userreal1;
  55     pi12 = ir->userreal2;
  56     d[0] = ir->userreal3;
  57     d[1] = ir->userreal4;
  58
  59     /* Check whether these constants have been set. */
  60     bGlas = (pi6 != 0) && (pi12 != 0) && (d[0] != 0) && (d[1] != 0);
  61
  62     if (bGlas) {
  63       if (ir->eDispCorr != edispcNO) {
  64         fatal_error(0,"Can not have Long Range C6 corrections and GLASMD");
  65       }
  66       rc   = max(fr->rvdw,fr->rlist);
  67       rc3  = rc*rc*rc;
  68       rc4  = rc3*rc;
  69       rc9  = rc3*rc3*rc3;
  70       rc10 = rc9*rc;
  71
  72       fprintf(log,
  73               "Constants for GLASMD: pi6 = %10g, pi12 = %10g\n"
  74               "                      d1  = %10g, d2   = %10g\n"
  75               "                      rc3 = %10g, rc4  = %10g\n"
  76               "                      rc9 = %10g, rc10 = %10g\n",
  77               pi6,pi12,d[0],d[1],rc3,rc4,rc9,rc10);
  78       if (d[0] > d[1])
  79         fatal_error(0,"d1 > d2 for GLASMD (check log file)");
  80
  81       snew(c6,atnr);
  82       snew(c12,atnr);
  83
  84       for(i=0; (i<atnr); i++) {
  85         c6[i]  = C6 (fr->nbfp,atnr,i,i);
  86         c12[i] = C12(fr->nbfp,atnr,i,i);
  87       }
  88     }
  89     else
  90       fprintf(stderr,"No glasmd!\n");
  91     bFirst = FALSE;
  92   }
  93
  94   if (bGlas) {
  95     wd=0;
  96     for(i=start; (i<start+homenr); i++) {
  97       ti   = type[i];
  98       if ((c6[ti] != 0) || (c12[ti] != 0)) {
  99         zi   = x[i][ZZ];
 100         cc6  = M_PI*sqrt(c6[ti]*pi6);
 101         cc12 = M_PI*sqrt(c12[ti]*pi12);
 102
 103         /* Use a factor for the sign, this saves taking absolute values */
 104         sign = 1;
 105         for(j=0; (j<2); j++) {
 106           dd = sign*(zi-d[j]);
 107           if (dd >= rc) {
 108             d3  = dd*dd*dd;
 109             d9  = d3*d3*d3;
 110             wdd = cc12/(45.0*d9) - cc6/(6.0*d3);
 111             d4  = d3*dd;
 112             d10 = d9*dd;
 113             fz  = sign*(cc12/(5.0*d10) - cc6/(2.0*d4));
 114           }
 115           else {
 116             wdd = cc12*(2.0/(9.0*rc9) - dd/(5.0*rc10)) -
 117               cc6*(2.0/(3.0*rc3) - dd/(2.0*rc4));
 118             fz  = sign*(cc12/(5.0*rc10)-cc6/(2.0*rc4));
 119           }
 120           wd       += wdd;
 121           f[i][ZZ] += fz;
 122           sign      = -sign;
 123         }
 124       }
 125     }
 126     ener[F_LJLR] = wd;
 127   }
 128 }
 129