src/gmxlib/shift_util.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.2.0
  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-2004, 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 Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <stdio.h>
  40 #include <math.h>
  41 #include "typedefs.h"
  42 #include "vec.h"
  43 #include "coulomb.h"
  44 #include "smalloc.h"
  45 #include "physics.h"
  46 #include "txtdump.h"
  47 #include "futil.h"
  48 #include "names.h"
  49 #include "writeps.h"
  50 #include "macros.h"
  51 #include "xvgr.h"
  52 #include "gmxfio.h"
  53
  54 #ifdef GMX_THREAD_MPI
  55 #include "thread_mpi/threads.h"
  56 #endif
  57
  58 #define p2(x) ((x)*(x))
  59 #define p3(x) ((x)*(x)*(x))
  60 #define p4(x) ((x)*(x)*(x)*(x))
  61
  62 static real A,A_3,B,B_4,C,c1,c2,c3,c4,c5,c6,One_4pi,FourPi_V,Vol,N0;
  63 #ifdef GMX_THREAD_MPI
  64 static tMPI_Thread_mutex_t shift_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
  65 #endif
  66
  67
  68 void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr)
  69 {
  70 #ifdef GMX_THREAD_MPI
  71   /* at the very least we shouldn't allow multiple threads to set these
  72      simulataneously */
  73   tMPI_Thread_mutex_lock(&shift_mutex);
  74 #endif
  75   /* A, B and C are recalculated in tables.c */
  76   if (r1 < rc) {
  77     A   = (2*r1-5*rc)/(p3(rc)*p2(rc-r1));
  78     B   = (4*rc-2*r1)/(p3(rc)*p3(rc-r1));
  79     /*C   = (10*rc*rc-5*rc*r1+r1*r1)/(6*rc*rc); Hermans Eq. not correct */
  80   }
  81   else
  82     gmx_fatal(FARGS,"r1 (%f) >= rc (%f) in %s, line %d",
  83               r1,rc,__FILE__,__LINE__);
  84
  85   A_3 = A/3.0;
  86   B_4 = B/4.0;
  87   C   = 1/rc-A_3*p3(rc-r1)-B_4*p4(rc-r1);
  88   N0  = 2.0*M_PI*p3(rc)*p3(rc-r1);
  89
  90   Vol     =(box[XX]*box[YY]*box[ZZ]);
  91   FourPi_V=4.0*M_PI/Vol;
  92
  93   if (debug) {
  94       fprintf(debug,"Constants for short-range and fourier stuff:\n"
  95               "r1 = %10.3f,  rc = %10.3f\n"
  96               "A  = %10.3e,  B  = %10.3e,  C  = %10.3e, FourPi_V = %10.3e\n",
  97               r1,rc,A,B,C,FourPi_V);
  98   }
  99
 100   /* Constants derived by Mathematica */
 101   c1 = -40*rc*rc    + 50*rc*r1    - 16*r1*r1;
 102   c2 =  60*rc       - 30*r1;
 103   c3 = -10*rc*rc*rc + 20*rc*rc*r1 - 13*rc*r1*r1 + 3*r1*r1*r1;
 104   c4 = -20*rc*rc    + 40*rc*r1    - 14*r1*r1;
 105   c5 = -c2;
 106   c6 = -5*rc*rc*r1  +  7*rc*r1*r1 - 2*r1*r1*r1;
 107
 108   if (debug) {
 109       fprintf(debug,"c1 = %10.3e,  c2 = %10.3e,  c3 = %10.3e\n"
 110               "c4 = %10.3e,  c5 = %10.3e,  c6 = %10.3e,  N0 = %10.3e\n",
 111               c1,c2,c3,c4,c5,c6,N0);
 112   }
 113
 114   One_4pi = 1.0/(4.0*M_PI);
 115 #ifdef GMX_THREAD_MPI
 116   tMPI_Thread_mutex_unlock(&shift_mutex);
 117 #endif
 118 }