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 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gromacs Runs On Most of All Computer Systems
  35  */
  36
  37 #ifndef _shift_util_h
  38 #define _shift_util_h
  39
  40 #ifdef HAVE_CONFIG_H
  41 #include <config.h>
  42 #endif
  43
  44 #include <math.h>
  45 #include "typedefs.h"
  46 #include "gmxcomplex.h"
  47 #include "physics.h"
  48
  49 /* Routines to set global constants for speeding up the calculation
  50  * of potentials, forces and gk-values.
  51  */
  52 extern void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr);
  53
  54 extern real gk(real k,real rc,real r1);
  55 /* Compute the Ghat function for a single k-value */
  56
  57 extern real gknew(real k,real rc,real r1);
  58 /* Compute the (new!) Ghat function for a single k-value */
  59
  60 extern void pr_scalar_gk(char *fn,int nx,int ny,int nz,rvec box,real ***ghat);
  61
  62 extern real calc_dx2(rvec xi,rvec xj,rvec box);
  63
  64 extern void calc_dx(rvec xi,rvec xj,rvec box,rvec dx);
  65
  66 extern real phi_sr(FILE *log,int nj,rvec x[],real charge[],real rc,real r1,
  67                    rvec box,real phi[],t_block *excl,rvec f_sr[],bool bOld);
  68
  69 extern real shiftfunction(real r1,real rc,real R);
  70
  71 extern real spreadfunction(real r1,real rc,real R);
  72
  73 extern real potential(real r1,real rc,real R);
  74
  75 extern void calc_ener(FILE *fp,char *title,bool bHeader,
  76                       int nmol,int natoms,
  77                       real phi[],real charge[],t_block *excl);
  78
  79 extern real shift_LRcorrection(FILE *fp,t_nsborder *nsb,
  80                                t_commrec *cr,t_forcerec *fr,
  81                                real charge[],t_block *excl,rvec x[],
  82                                bool bOld,matrix box,matrix lrvir);
  83 /* Calculate the self energy and forces
  84  * when using long range electrostatics methods.
  85  * Part of this is a constant, it is computed only once and stored in
  86  * a local variable. The remainder is computed every step.
  87  * PBC is taken into account. (Erik L.)
  88  */
  89
  90 extern void calc_weights(int iatom,int nx,int ny,int nz,
  91                          rvec x,rvec box,rvec invh,ivec ixyz,real WXYZ[]);
  92
  93 static void calc_lll(rvec box,rvec lll)
  94 {
  95   lll[XX] = 2.0*M_PI/box[XX];
  96   lll[YY] = 2.0*M_PI/box[YY];
  97   lll[ZZ] = 2.0*M_PI/box[ZZ];
  98 }
  99
 100 static void calc_k(rvec lll,int ix,int iy,int iz,int nx,int ny,int nz,rvec k)
 101 {
 102 #define IDX(i,n,x)  (i<=n/2) ? (i*x) : ((i-n)*x)
 103   k[XX] = IDX(ix,nx,lll[XX]);
 104   k[YY] = IDX(iy,ny,lll[YY]);
 105   k[ZZ] = IDX(iz,nz,lll[ZZ]);
 106 #undef IDX
 107 }
 108
 109 /******************************************************************
 110  *
 111  *   PLOTTING ROUTINES FOR DEBUGGING
 112  *
 113  ******************************************************************/
 114
 115 extern void plot_phi(char *fn,rvec box,int natoms,rvec x[],real phi[]);
 116 /* Plot potential (or whatever) in a postscript matrix */
 117
 118 extern void print_phi(char *fn,int natoms,rvec x[],real phi[]);
 119 /* Print to a text file in x y phi format */
 120
 121 extern void plot_qtab(char *fn,int nx,int ny,int nz,real ***qtab);
 122 /* Plot a charge table to a postscript matrix */
 123
 124 extern void write_grid_pqr(char *fn,int nx,int ny,int nz,real ***phi);
 125 extern void write_pqr(char *fn,t_atoms *atoms,rvec x[],real phi[],real dx);
 126 /* Write a pdb file where the potential phi is printed as B-factor (for
 127  * viewing with rasmol). All atoms are moved over a distance dx in the X
 128  * direction, to enable viewing of two data sets simultaneously with rasmol
 129  */
 130
 131 /******************************************************************
 132  *
 133  *   ROUTINES FOR GHAT MANIPULATION
 134  *
 135  ******************************************************************/
 136
 137 extern void symmetrize_ghat(int nx,int ny,int nz,real ***ghat);
 138 /* Symmetrize the Ghat function. It is assumed that the
 139  * first octant of the Ghat function is either read or generated
 140  * (all k-vectors from 0..nx/2 0..ny/2 0..nz/2).
 141  * Since Gk depends on the absolute value of k only,
 142  * symmetry operations may shorten the time to generate it.
 143  */
 144
 145 extern void mk_ghat(FILE *fp,int nx,int ny,int nz,real ***ghat,
 146                     rvec box,real r1,real rc,bool bSym,bool bOld);
 147 /* Generate a Ghat function from scratch. The ghat grid should
 148  * be allocated using the mk_rgrid function. When bSym, only
 149  * the first octant of the function is generated by direct calculation
 150  * and the above mentioned function is called for computing the rest.
 151  * When !bOld a new experimental function form will be used.
 152  */
 153
 154 extern real ***rd_ghat(FILE *log,char *fn,ivec igrid,rvec gridspacing,
 155                        rvec beta,int *porder,real *rshort,real *rlong);
 156 /* Read a Ghat function from a file as generated by the program
 157  * mk_ghat. The grid size (number of grid points) is returned in
 158  * igrid, the space between grid points in gridspacing,
 159  * beta is a constant that determines the contribution of first
 160  * and second neighbours in the grid to the force
 161  * (See Luty et al. JCP 103 (1995) 3014)
 162  * porder determines whether 8 (when porder = 1) or 27 (when
 163  * porder = 2) neighbouring grid points are used for spreading
 164  * the charge.
 165  * rshort and rlong are the lengths used for generating the Ghat
 166  * function.
 167  */
 168
 169 extern void wr_ghat(char *fn,int n1max,int n2max,int n3max,real h1,
 170                     real h2,real h3,real ***ghat,int nalias,
 171                     int porder,int niter,bool bSym,rvec beta,
 172                     real r1,real rc,real pval,real zval,real eref,real qopt);
 173 /* Write a ghat file. (see above) */
 174
 175 extern void pr_scalar_gk(char *fn,int nx,int ny,int nz,rvec box,real ***ghat);
 176
 177 extern real analyse_diff(FILE *log,char *label,
 178                          int natom,rvec ffour[],rvec fpppm[],
 179                          real phi_f[],real phi_p[],real phi_sr[],
 180                          char *fcorr,char *pcorr,
 181                          char *ftotcorr,char *ptotcorr);
 182 /* Analyse difference between forces from fourier (_f) and other (_p)
 183  * LR solvers (and potential also).
 184  * If the filenames are given, xvgr files are written.
 185  * returns the root mean square error in the force.
 186  */
 187
 188 #endif