include/shift_util.h

   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  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29
  30 #ifndef _shift_util_h
  31 #define _shift_util_h
  32
  33 static char *SRCID_shift_util_h = "$Id$";
  34
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "typedefs.h"
  41 #include "complex.h"
  42 #include "physics.h"
  43
  44 /* Routines to set global constants for speeding up the calculation
  45  * of potentials, forces and gk-values.
  46  */
  47 extern void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr);
  48
  49 extern real gk(real k,real rc,real r1);
  50 /* Compute the Ghat function for a single k-value */
  51
  52 extern real gknew(real k,real rc,real r1);
  53 /* Compute the (new!) Ghat function for a single k-value */
  54
  55 extern void pr_scalar_gk(char *fn,int nx,int ny,int nz,rvec box,real ***ghat);
  56
  57 extern real calc_dx2(rvec xi,rvec xj,rvec box);
  58
  59 extern void calc_dx(rvec xi,rvec xj,rvec box,rvec dx);
  60
  61 extern real phi_sr(FILE *log,int nj,rvec x[],real charge[],real rc,real r1,
  62                    rvec box,real phi[],t_block *excl,rvec f_sr[],bool bOld);
  63
  64 extern real shiftfunction(real r1,real rc,real R);
  65
  66 extern real spreadfunction(real r1,real rc,real R);
  67
  68 extern real potential(real r1,real rc,real R);
  69
  70 extern void calc_ener(FILE *fp,char *title,bool bHeader,
  71                       int nmol,int natoms,
  72                       real phi[],real charge[],t_block *excl);
  73
  74 extern real shift_LRcorrection(FILE *fp,t_nsborder *nsb,
  75                                t_commrec *cr,t_forcerec *fr,
  76                                real charge[],t_block *excl,rvec x[],
  77                                bool bOld,matrix box,matrix lrvir);
  78 /* Calculate the self energy and forces
  79  * when using long range electrostatics methods.
  80  * Part of this is a constant, it is computed only once and stored in
  81  * a local variable. The remainder is computed every step.
  82  * PBC is taken into account. (Erik L.)
  83  */
  84
  85 extern void calc_weights(int iatom,int nx,int ny,int nz,
  86                          rvec x,rvec box,rvec invh,ivec ixyz,real WXYZ[]);
  87
  88 static void calc_lll(rvec box,rvec lll)
  89 {
  90   lll[XX] = 2.0*M_PI/box[XX];
  91   lll[YY] = 2.0*M_PI/box[YY];
  92   lll[ZZ] = 2.0*M_PI/box[ZZ];
  93 }
  94
  95 static void calc_k(rvec lll,int ix,int iy,int iz,int nx,int ny,int nz,rvec k)
  96 {
  97 #define IDX(i,n,x)  (i<=n/2) ? (i*x) : ((i-n)*x)
  98   k[XX] = IDX(ix,nx,lll[XX]);
  99   k[YY] = IDX(iy,ny,lll[YY]);
 100   k[ZZ] = IDX(iz,nz,lll[ZZ]);
 101 #undef IDX
 102 }
 103
 104 /******************************************************************
 105  *
 106  *   PLOTTING ROUTINES FOR DEBUGGING
 107  *
 108  ******************************************************************/
 109
 110 extern void plot_phi(char *fn,rvec box,int natoms,rvec x[],real phi[]);
 111 /* Plot potential (or whatever) in a postscript matrix */
 112
 113 extern void print_phi(char *fn,int natoms,rvec x[],real phi[]);
 114 /* Print to a text file in x y phi format */
 115
 116 extern void plot_qtab(char *fn,int nx,int ny,int nz,real ***qtab);
 117 /* Plot a charge table to a postscript matrix */
 118
 119 extern void write_grid_pqr(char *fn,int nx,int ny,int nz,real ***phi);
 120 extern void write_pqr(char *fn,t_atoms *atoms,rvec x[],real phi[],real dx);
 121 /* Write a pdb file where the potential phi is printed as B-factor (for
 122  * viewing with rasmol). All atoms are moved over a distance dx in the X
 123  * direction, to enable viewing of two data sets simultaneously with rasmol
 124  */
 125
 126 /******************************************************************
 127  *
 128  *   ROUTINES FOR GHAT MANIPULATION
 129  *
 130  ******************************************************************/
 131
 132 extern void symmetrize_ghat(int nx,int ny,int nz,real ***ghat);
 133 /* Symmetrize the Ghat function. It is assumed that the
 134  * first octant of the Ghat function is either read or generated
 135  * (all k-vectors from 0..nx/2 0..ny/2 0..nz/2).
 136  * Since Gk depends on the absolute value of k only,
 137  * symmetry operations may shorten the time to generate it.
 138  */
 139
 140 extern void mk_ghat(FILE *fp,int nx,int ny,int nz,real ***ghat,
 141                     rvec box,real r1,real rc,bool bSym,bool bOld);
 142 /* Generate a Ghat function from scratch. The ghat grid should
 143  * be allocated using the mk_rgrid function. When bSym, only
 144  * the first octant of the function is generated by direct calculation
 145  * and the above mentioned function is called for computing the rest.
 146  * When !bOld a new experimental function form will be used.
 147  */
 148
 149 extern real ***rd_ghat(FILE *log,char *fn,ivec igrid,rvec gridspacing,
 150                        rvec beta,int *porder,real *rshort,real *rlong);
 151 /* Read a Ghat function from a file as generated by the program
 152  * mk_ghat. The grid size (number of grid points) is returned in
 153  * igrid, the space between grid points in gridspacing,
 154  * beta is a constant that determines the contribution of first
 155  * and second neighbours in the grid to the force
 156  * (See Luty et al. JCP 103 (1995) 3014)
 157  * porder determines whether 8 (when porder = 1) or 27 (when
 158  * porder = 2) neighbouring grid points are used for spreading
 159  * the charge.
 160  * rshort and rlong are the lengths used for generating the Ghat
 161  * function.
 162  */
 163
 164 extern void wr_ghat(char *fn,int n1max,int n2max,int n3max,real h1,
 165                     real h2,real h3,real ***ghat,int nalias,
 166                     int porder,int niter,bool bSym,rvec beta,
 167                     real r1,real rc,real pval,real zval,real eref,real qopt);
 168 /* Write a ghat file. (see above) */
 169
 170 extern void pr_scalar_gk(char *fn,int nx,int ny,int nz,rvec box,real ***ghat);
 171
 172 extern real analyse_diff(FILE *log,char *label,
 173                          int natom,rvec ffour[],rvec fpppm[],
 174                          real phi_f[],real phi_p[],real phi_sr[],
 175                          char *fcorr,char *pcorr,
 176                          char *ftotcorr,char *ptotcorr);
 177 /* Analyse difference between forces from fourier (_f) and other (_p)
 178  * LR solvers (and potential also).
 179  * If the filenames are given, xvgr files are written.
 180  * returns the root mean square error in the force.
 181  */
 182
 183 #endif