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