src/gmxlib/shift_util.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdio.h>
  43 #include <math.h>
  44 #include "typedefs.h"
  45 #include "vec.h"
  46 #include "coulomb.h"
  47 #include "smalloc.h"
  48 #include "physics.h"
  49 #include "txtdump.h"
  50 #include "futil.h"
  51 #include "names.h"
  52 #include "writeps.h"
  53 #include "macros.h"
  54 #include "xvgr.h"
  55 #include "gmxfio.h"
  56
  57 #ifdef GMX_THREAD_MPI
  58 #include "thread_mpi/threads.h"
  59 #endif
  60
  61 #define p2(x) ((x)*(x))
  62 #define p3(x) ((x)*(x)*(x))
  63 #define p4(x) ((x)*(x)*(x)*(x))
  64
  65 static real A,A_3,B,B_4,C,c1,c2,c3,c4,c5,c6,One_4pi,FourPi_V,Vol,N0;
  66 #ifdef GMX_THREAD_MPI
  67 static tMPI_Thread_mutex_t shift_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
  68 #endif
  69
  70
  71 void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr)
  72 {
  73 #ifdef GMX_THREAD_MPI
  74   /* at the very least we shouldn't allow multiple threads to set these
  75      simulataneously */
  76   tMPI_Thread_mutex_lock(&shift_mutex);
  77 #endif
  78   /* A, B and C are recalculated in tables.c */
  79   if (r1 < rc) {
  80     A   = (2*r1-5*rc)/(p3(rc)*p2(rc-r1));
  81     B   = (4*rc-2*r1)/(p3(rc)*p3(rc-r1));
  82     /*C   = (10*rc*rc-5*rc*r1+r1*r1)/(6*rc*rc); Hermans Eq. not correct */
  83   }
  84   else
  85     gmx_fatal(FARGS,"r1 (%f) >= rc (%f) in %s, line %d",
  86               r1,rc,__FILE__,__LINE__);
  87
  88   A_3 = A/3.0;
  89   B_4 = B/4.0;
  90   C   = 1/rc-A_3*p3(rc-r1)-B_4*p4(rc-r1);
  91   N0  = 2.0*M_PI*p3(rc)*p3(rc-r1);
  92
  93   Vol     =(box[XX]*box[YY]*box[ZZ]);
  94   FourPi_V=4.0*M_PI/Vol;
  95
  96   if (debug) {
  97       fprintf(debug,"Constants for short-range and fourier stuff:\n"
  98               "r1 = %10.3f,  rc = %10.3f\n"
  99               "A  = %10.3e,  B  = %10.3e,  C  = %10.3e, FourPi_V = %10.3e\n",
 100               r1,rc,A,B,C,FourPi_V);
 101   }
 102
 103   /* Constants derived by Mathematica */
 104   c1 = -40*rc*rc    + 50*rc*r1    - 16*r1*r1;
 105   c2 =  60*rc       - 30*r1;
 106   c3 = -10*rc*rc*rc + 20*rc*rc*r1 - 13*rc*r1*r1 + 3*r1*r1*r1;
 107   c4 = -20*rc*rc    + 40*rc*r1    - 14*r1*r1;
 108   c5 = -c2;
 109   c6 = -5*rc*rc*r1  +  7*rc*r1*r1 - 2*r1*r1*r1;
 110
 111   if (debug) {
 112       fprintf(debug,"c1 = %10.3e,  c2 = %10.3e,  c3 = %10.3e\n"
 113               "c4 = %10.3e,  c5 = %10.3e,  c6 = %10.3e,  N0 = %10.3e\n",
 114               c1,c2,c3,c4,c5,c6,N0);
 115   }
 116
 117   One_4pi = 1.0/(4.0*M_PI);
 118 #ifdef GMX_THREAD_MPI
 119   tMPI_Thread_mutex_unlock(&shift_mutex);
 120 #endif
 121 }