src/gromacs/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     {
  78         A   = (2*r1-5*rc)/(p3(rc)*p2(rc-r1));
  79         B   = (4*rc-2*r1)/(p3(rc)*p3(rc-r1));
  80         /*C   = (10*rc*rc-5*rc*r1+r1*r1)/(6*rc*rc); Hermans Eq. not correct */
  81     }
  82     else
  83     {
  84         gmx_fatal(FARGS, "r1 (%f) >= rc (%f) in %s, line %d",
  85                   r1, rc, __FILE__, __LINE__);
  86     }
  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     {
  98         fprintf(debug, "Constants for short-range and fourier stuff:\n"
  99                 "r1 = %10.3f,  rc = %10.3f\n"
 100                 "A  = %10.3e,  B  = %10.3e,  C  = %10.3e, FourPi_V = %10.3e\n",
 101                 r1, rc, A, B, C, FourPi_V);
 102     }
 103
 104     /* Constants derived by Mathematica */
 105     c1 = -40*rc*rc    + 50*rc*r1    - 16*r1*r1;
 106     c2 =  60*rc       - 30*r1;
 107     c3 = -10*rc*rc*rc + 20*rc*rc*r1 - 13*rc*r1*r1 + 3*r1*r1*r1;
 108     c4 = -20*rc*rc    + 40*rc*r1    - 14*r1*r1;
 109     c5 = -c2;
 110     c6 = -5*rc*rc*r1  +  7*rc*r1*r1 - 2*r1*r1*r1;
 111
 112     if (debug)
 113     {
 114         fprintf(debug, "c1 = %10.3e,  c2 = %10.3e,  c3 = %10.3e\n"
 115                 "c4 = %10.3e,  c5 = %10.3e,  c6 = %10.3e,  N0 = %10.3e\n",
 116                 c1, c2, c3, c4, c5, c6, N0);
 117     }
 118
 119     One_4pi = 1.0/(4.0*M_PI);
 120 #ifdef GMX_THREAD_MPI
 121     tMPI_Thread_mutex_unlock(&shift_mutex);
 122 #endif
 123 }