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