set_shift_consts(FILE *log,real r1,real rc,rvec box,
t_forcerec *fr);
-real
-shift_LRcorrection(FILE *fp,int start,int natoms,
- t_commrec *cr,t_forcerec *fr,
- real charge[],t_blocka *excl,rvec x[],
- gmx_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.)
- */
-
#ifdef __cplusplus
}
#endif
+++ /dev/null
-/*
- *
- * 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 _pppm_h
-#define _pppm_h
-
-#include <stdio.h>
-#include "typedefs.h"
-#include "gmxcomplex.h"
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-int gmx_pppm_init(FILE *log, t_commrec *cr,
- const output_env_t oenv, gmx_bool bVerbose,
- gmx_bool bOld, matrix box,
- char *ghatfn, t_inputrec *ir,
- gmx_bool bReproducible);
-/* Setup stuff for PPPM.
- * Either reads a ghat function from file (when the file exists)
- * or generate a ghat function from scratch.
- */
-
-int gmx_pppm_do(FILE *log, gmx_pme_t pme,
- gmx_bool bVerbose,
- rvec x[], rvec f[],
- real charge[], rvec box,
- real phi[], t_commrec *cr,
- int start, int nt,
- t_nrnb *nrnb,
- int pme_order, real *energy);
-/* Do a PPPM calculation for the long range electrostatics. */
-
-/******************************************************************
- *
- * ROUTINES FOR GHAT MANIPULATION
- *
- ******************************************************************/
-
-real gk(real k,real rc,real r1);
-/* Compute the Ghat function for a single k-value */
-
-real gknew(real k,real rc,real r1);
-/* Compute the (new!) Ghat function for a single k-value */
-
-void pr_scalar_gk(const char *fn,const output_env_t oenv,
- int nx,int ny,int nz, rvec box,real ***ghat);
-
-void mk_ghat(FILE *fp,int nx,int ny,int nz,
- real ***ghat, rvec box,real r1,real rc,gmx_bool bSym,gmx_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.
- */
-
-real ***rd_ghat(FILE *log,const output_env_t oenv,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.
- */
-
-void wr_ghat(const char *fn,const output_env_t oenv,
- int n1max,int n2max, int n3max,real h1,
- real h2,real h3,real ***ghat,int nalias,
- int porder,int niter,gmx_bool bSym,rvec beta,
- real r1,real rc,real pval,real zval,real eref,real qopt);
-/* Write a ghat file. (see above) */
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
* separately (through the implicit_solvent option).
*/
enum {
- eelCUT, eelRF, eelGRF, eelPME, eelEWALD, eelPPPM,
+ eelCUT, eelRF, eelGRF, eelPME, eelEWALD, eelP3M_AD,
eelPOISSON, eelSWITCH, eelSHIFT, eelUSER, eelGB_NOTUSED, eelRF_NEC, eelENCADSHIFT,
eelPMEUSER, eelPMESWITCH, eelPMEUSERSWITCH, eelRF_ZERO, eelNR
};
#define EEL_RF(e) ((e) == eelRF || (e) == eelGRF || (e) == eelRF_NEC || (e) == eelRF_ZERO )
-#define EEL_PME(e) ((e) == eelPME || (e) == eelPMESWITCH || (e) == eelPMEUSER || (e) == eelPMEUSERSWITCH)
-#define EEL_FULL(e) (EEL_PME(e) || (e) == eelPPPM || (e) == eelPOISSON || (e) == eelEWALD)
+#define EEL_PME(e) ((e) == eelPME || (e) == eelPMESWITCH || (e) == eelPMEUSER || (e) == eelPMEUSERSWITCH || (e) == eelP3M_AD)
+#define EEL_FULL(e) (EEL_PME(e) || (e) == eelPOISSON || (e) == eelEWALD)
#define EEL_SWITCHED(e) ((e) == eelSWITCH || (e) == eelSHIFT || (e) == eelENCADSHIFT || (e) == eelPMESWITCH || (e) == eelPMEUSERSWITCH)
real calc_temp(real ekin,real nrdf);
/* Calculate the temperature */
-real calc_pres(int ePBC,int nwall,matrix box,
- tensor ekin,tensor vir,tensor pres,real Elr);
+real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
+ tensor pres);
/* Calculate the pressure tensor, returns the scalar pressure.
- * The unit of pressure is bar, If Elr != 0
- * a long range correction based on Ewald/PPPM is made (see c-code)
+ * The unit of pressure is bar.
*/
void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
reference - in most cases PME will perform much better.</dd>
<dt><b><!--Idx-->PME<!--EIdx--></b></dt>
-<dd>Fast Particle-Mesh Ewald electrostatics. Direct space is similar
+<dd>Fast smooth Particle-Mesh Ewald (SPME) electrostatics. Direct space is similar
to the Ewald sum, while the reciprocal part is performed with
FFTs. Grid dimensions are controlled with <b>fourierspacing</b> and the
interpolation order with <b>pme-order</b>. With a grid spacing of 0.1
parallelizes better than the FFT, so try decreasing grid dimensions
while increasing interpolation.</dd>
-<dt><b><!--Idx-->PPPM<!--EIdx--></b></dt>
-<dd>Particle-Particle Particle-Mesh algorithm for long range
-electrostatic interactions.
-Use for example <b>rlist</b><tt>=0.9</tt>, <b>rcoulomb</b><tt>=0.9</TT>.
-The grid dimensions are controlled by <b>fourierspacing</b>.
-Reasonable grid spacing for PPPM is 0.05-0.1 nm.
-See <tt>Shift</tt> for the details of the particle-particle potential.
-<br>
-NOTE: PPPM is not functional in the current version, but we plan to implement
-PPPM through a small modification of the PME code.</dd>
+<dt><b><!--Idx-->P3M-AD<!--EIdx--></b></dt>
+<dd>Particle-Particle Particle-Mesh algorithm with analytical derivative
+for for long range electrostatic interactions. The method and code
+is identical to SPME, except that the influence function is optimized
+for the grid. This gives a slight increase in accuracy.</dd>
<dt><b><!--Idx-->Reaction-Field<!--EIdx--></b></dt>
<dd>Reaction field with Coulomb cut-off <b>rcoulomb</b>,
<h3>Ewald</h3>
<dl>
<dt><b>fourierspacing: (0.12) [nm]</b></dt>
-<dd>The maximum grid spacing for the FFT grid when using PPPM or PME.
+<dd>The maximum grid spacing for the FFT grid when using PME or P3M.
For ordinary Ewald the spacing times the box dimensions determines the
highest magnitude to use in each direction. In all cases
each direction can be overridden by entering a non-zero value for
<dt><b>fourier-nx (0) ; fourier-ny (0) ; fourier-nz: (0)</b></dt>
<dd>Highest magnitude of wave vectors in reciprocal space when using Ewald.</dd>
-<dd>Grid size when using PPPM or PME. These values override
+<dd>Grid size when using PME or P3M. These values override
<b>fourierspacing</b> per direction. The best choice is powers of
2, 3, 5 and 7. Avoid large primes.</dd>
"M. Hoefling, N. Lima, D. Haenni, C.A.M. Seidel, B. Schuler, H. Grubmuller",
"Structural Heterogeneity and Quantitative FRET Efficiency Distributions of Polyprolines through a Hybrid Atomistic Simulation and Monte Carlo Approach",
"PLoS ONE",
- 6, 2011, "e19791"
- }
+ 6, 2011, "e19791" },
+ { "Hockney1988",
+ "R. W. Hockney and J. W. Eastwood",
+ "Computer simulation using particles",
+ "IOP, Bristol",
+ 1, 1988, "1" },
+ { "Ballenegger2012",
+ "V. Ballenegger, J.J. Cerda, and C. Holm",
+ "How to Convert SPME to P3M: Influence Functions and Error Estimates",
+ "J. Chem. Theory Comput.",
+ 8, 2012, "936-947" }
};
#define NSTR (int)asize(citedb)
const char *eel_names[eelNR+1] = {
"Cut-off", "Reaction-Field", "Generalized-Reaction-Field",
- "PME", "Ewald", "PPPM", "Poisson", "Switch", "Shift", "User",
+ "PME", "Ewald", "P3M-AD", "Poisson", "Switch", "Shift", "User",
"Generalized-Born", "Reaction-Field-nec", "Encad-shift",
"PME-User", "PME-Switch", "PME-User-Switch",
"Reaction-Field-zero", NULL
#include "writeps.h"
#include "macros.h"
#include "xvgr.h"
-#include "pppm.h"
#include "gmxfio.h"
#include "thread_mpi.h"
tMPI_Thread_mutex_unlock(&shift_mutex);
#endif
}
-
-real gk(real k,real rc,real r1)
-/* Spread function in Fourier space. Eqs. 56-64 */
-{
- real N,gg;
-
- N = N0*p4(k);
-
- /* c1 thru c6 consts are global variables! */
- gg = (FourPi_V / (N*k)) *
- ( c1*k*cos(k*rc) + (c2+c3*k*k)*sin(k*rc) +
- c4*k*cos(k*r1) + (c5+c6*k*k)*sin(k*r1) );
-
- return gg;
-}
-
-real gknew(real k,real rc,real r1)
-{
- real rck,rck2;
-
- rck = rc*k;
- rck2= rck*rck;
-
- return -15.0*((rck2-3.0)*sin(rck) + 3*rck*cos(rck))/(Vol*rck2*rck2*rck);
-}
-
-real calc_dx2dx(rvec xi,rvec xj,rvec box,rvec dx)
-{
- int m;
- real ddx,dx2;
-
- dx2=0;
- for(m=0; (m<DIM); m++) {
- ddx=xj[m]-xi[m];
- if (ddx < -0.5*box[m])
- ddx+=box[m];
- else if (ddx >= 0.5*box[m])
- ddx-=box[m];
- dx[m]=ddx;
- dx2 += ddx*ddx;
- }
- return dx2;
-}
-
-real shiftfunction(real r1,real rc,real R)
-{
- real dr;
-
- if (R <= r1)
- return 0.0;
- else if (R >= rc)
- return -1.0/(R*R);
-
- dr=R-r1;
-
- return A*dr*dr+B*dr*dr*dr;
-}
-
-real potential(real r1,real rc,real R)
-{
- if (R < r1)
- return (1.0/R-C);
- else if (R <= rc)
- return (1.0/R - A_3*p3(R-r1) - B_4*p4(R-r1) - C);
- else
- return 0.0;
-}
-
-
-
-real shift_LRcorrection(FILE *fp,int start,int natoms,
- t_commrec *cr,t_forcerec *fr,
- real charge[],t_blocka *excl,rvec x[],
- gmx_bool bOld,matrix box,matrix lr_vir)
-{
- static gmx_bool bFirst=TRUE;
- static real Vself;
- int i,i1,i2,j,k,m,iv,jv;
- int *AA;
- double qq; /* Necessary for precision */
- double isp=0.564189583547756;
- real qi,dr,dr2,dr_1,dr_3,fscal,Vexcl,qtot=0;
- rvec df,dx;
- real r1=fr->rcoulomb_switch;
- real rc=fr->rcoulomb;
- ivec shift;
-
- if (bFirst) {
- qq =0;
- for(i=start; (i<start+natoms); i++)
- qq += charge[i]*charge[i];
-
- /* Obsolete, until we implement dipole and charge corrections.
- qtot=0;
- for(i=0;i<nsb->natoms;i++)
- qtot+=charge[i];
- */
-
- Vself = 0.5*C*ONE_4PI_EPS0*qq;
- if(debug) {
- fprintf(fp,"Long Range corrections for shifted interactions:\n");
- fprintf(fp,"r1 = %g, rc=%g\n",r1,rc);
- fprintf(fp,"start=%d,natoms=%d\n",start,natoms);
- fprintf(fp,"qq = %g, Vself=%g\n",qq,Vself);
- }
-
- }
- AA = excl->a;
- Vexcl = 0;
-
- for(i=start; (i<start+natoms); i++) {
- /* Initiate local variables (for this i-particle) to 0 */
- i1 = excl->index[i];
- i2 = excl->index[i+1];
- qi = charge[i]*ONE_4PI_EPS0;
-
- /* Loop over excluded neighbours */
- for(j=i1; (j<i2); j++) {
- k = AA[j];
- /*
- * First we must test whether k <> i, and then, because the
- * exclusions are all listed twice i->k and k->i we must select
- * just one of the two.
- * As a minor optimization we only compute forces when the charges
- * are non-zero.
- */
- if (k > i) {
- qq = qi*charge[k];
- if (qq != 0.0) {
- dr2 = 0;
- rvec_sub(x[i],x[k],dx);
- for(m=DIM-1; m>=0; m--) {
- if (dx[m] > 0.5*box[m][m])
- rvec_dec(dx,box[m]);
- else if (dx[m] < -0.5*box[m][m])
- rvec_inc(dx,box[m]);
-
- dr2 += dx[m]*dx[m];
- }
- dr_1 = gmx_invsqrt(dr2);
- dr = 1.0/dr_1;
- dr_3 = dr_1*dr_1*dr_1;
- /* Compute exclusion energy and scalar force */
-
- Vexcl += qq*(dr_1-potential(r1,rc,dr));
- fscal = qq*(-shiftfunction(r1,rc,dr))*dr_3;
-
- if ((fscal != 0) && debug)
- fprintf(debug,"i: %d, k: %d, dr: %.3f fscal: %.3f\n",i,k,dr,fscal);
-
- /* The force vector is obtained by multiplication with the
- * distance vector
- */
- svmul(fscal,dx,df);
- rvec_inc(fr->f_novirsum[k],df);
- rvec_dec(fr->f_novirsum[i],df);
- for(iv=0;iv<DIM;iv++)
- for(jv=0;jv<DIM;jv++)
- lr_vir[iv][jv]+=0.5*dx[iv]*df[jv];
- }
- }
- }
- }
- if (bFirst && debug)
- fprintf(fp,"Long Range correction: Vexcl=%g\n",Vexcl);
-
- bFirst = FALSE;
- /* Return the correction to the energy */
- return (-(Vself+Vexcl));
-}
-
-real phi_aver(int natoms,real phi[])
-{
- real phitot;
- int i;
-
- phitot=0;
- for(i=0; (i<natoms); i++)
- phitot=phitot+phi[i];
-
- return (phitot/natoms);
-}
-
svmul(ir->wall_ewald_zfac,box[ZZ],box[ZZ]);
max_spacing = calc_grid(stdout,box,opts->fourierspacing,
&(ir->nkx),&(ir->nky),&(ir->nkz));
- if ((ir->coulombtype == eelPPPM) && (max_spacing > 0.1)) {
- set_warning_line(wi,mdparin,-1);
- warning_note(wi,"Grid spacing larger then 0.1 while using PPPM.");
- }
}
if (ir->ePull != epullNO)
#include "disre.h"
#include "orires.h"
#include "dihre.h"
-#include "pppm.h"
#include "pme.h"
#include "mdatoms.h"
#include "repl_ex.h"
#include "disre.h"
#include "orires.h"
#include "dihre.h"
-#include "pppm.h"
#include "pme.h"
#include "mdatoms.h"
#include "qmmm.h"
warning_note(wi,"Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
}
- sprintf(err_buf,"Free-energy not implemented for Ewald and PPPM");
- CHECK((ir->coulombtype==eelEWALD || ir->coulombtype==eelPPPM)
- && (ir->efep!=efepNO));
+ sprintf(err_buf,"Free-energy not implemented for Ewald");
+ CHECK((ir->coulombtype==eelEWALD) && (ir->efep!=efepNO));
sprintf(err_buf,"Twin-range neighbour searching (NS) with simple NS"
" algorithm not implemented");
(trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
- sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
- CHECK(ir->coulombtype == eelPPPM);
-
if (epcPARRINELLORAHMAN == ir->epct && opts->bGenVel)
{
sprintf(warn_buf,
warning(wi,warn_buf);
}
}
- else if (ir->coulombtype == eelPPPM)
- {
- sprintf(warn_buf,"The pressure with PPPM is incorrect, if you need the pressure use PME");
- warning(wi,warn_buf);
- }
if (EI_VV(ir->eI))
{
/* ELECTROSTATICS */
/* More checks are in triple check (grompp.c) */
- if (ir->coulombtype == eelPPPM)
- {
- warning_error(wi,"PPPM is not functional in the current version, we plan to implement PPPM through a small modification of the PME code");
- }
if (ir->coulombtype == eelSWITCH) {
sprintf(warn_buf,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
eel_names[ir->coulombtype]);
CHECK(ir->rcoulomb > ir->rlist);
} else {
- if (ir->coulombtype == eelPME) {
+ if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD) {
sprintf(err_buf,
"With coulombtype = %s, rcoulomb must be equal to rlist\n"
"If you want optimal energy conservation or exact integration use %s",
#include "disre.h"
#include "orires.h"
#include "dihre.h"
-#include "pppm.h"
#include "pme.h"
#include "mdatoms.h"
#include "repl_ex.h"
}
}
- /* Initiate PPPM if necessary */
- if (fr->eeltype == eelPPPM)
- {
- if (mdatoms->nChargePerturbed)
- {
- gmx_fatal(FARGS,"Free energy with %s is not implemented",
- eel_names[fr->eeltype]);
- }
- status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
- getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
- if (status != 0)
- {
- gmx_fatal(FARGS,"Error %d initializing PPPM",status);
- }
- }
-
if (EEL_PME(fr->eeltype))
{
ewaldcoeff = fr->ewaldcoeff;
alpha *= ekind->tcstat[0].ekinscalef_nhc;
msmul(ekind->ekin,alpha,ekinmod);
- /* for now, we use Elr = 0, because if you want to get it right, you
- really should be using PME. Maybe print a warning? */
-
- pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres,0.0) + pcorr;
+ pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres) + pcorr;
vol = det(box);
GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
*/
real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
- tensor pres,real Elr)
+ tensor pres)
{
int n,m;
- real fac,Plr;
+ real fac;
if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
clear_mat(pres);
* het systeem...
*/
- /* Long range correction for periodic systems, see
- * Neumann et al. JCP
- * divide by 6 because it is multiplied by fac later on.
- * If Elr = 0, no correction is made.
- */
-
- /* This formula should not be used with Ewald or PME,
- * where the full long-range virial is calculated. EL 990823
- */
- Plr = Elr/6.0;
-
fac=PRESFAC*2.0/det(box);
for(n=0; (n<DIM); n++)
for(m=0; (m<DIM); m++)
- pres[n][m]=(ekin[n][m]-vir[n][m]+Plr)*fac;
+ pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
if (debug) {
pr_rvecs(debug,0,"PC: pres",pres,DIM);
#include "mshift.h"
#include "txtdump.h"
#include "coulomb.h"
-#include "pppm.h"
#include "pme.h"
#include "mdrun.h"
#include "domdec.h"
Vcorr = 0;
}
}
- else
- {
- Vcorr = shift_LRcorrection(fplog,md->start,md->homenr,cr,fr,
- md->chargeA,excl,x,TRUE,box,
- fr->vir_el_recip);
- }
dvdlambda = 0;
status = 0;
switch (fr->eeltype)
{
- case eelPPPM:
- status = gmx_pppm_do(fplog,fr->pmedata,FALSE,x,fr->f_novirsum,
- md->chargeA,
- box_size,fr->phi,cr,md->start,md->homenr,
- nrnb,ir->pme_order,&Vlr);
- break;
case eelPME:
case eelPMESWITCH:
case eelPMEUSER:
case eelPMEUSERSWITCH:
+ case eelP3M_AD:
if (cr->duty & DUTY_PME)
{
if (fr->n_tpi == 0 || (flags & GMX_FORCE_STATECHANGED))
{
if (fp)
fprintf(fp,"Will do PME sum in reciprocal space.\n");
- please_cite(fp,"Essmann95a");
+ if (ir->coulombtype == eelP3M_AD)
+ {
+ please_cite(fp,"Hockney1988");
+ please_cite(fp,"Ballenegger2012");
+ }
+ else
+ {
+ please_cite(fp,"Essmann95a");
+ }
if (ir->ewald_geometry == eewg3DC)
{
{
init_generalized_rf(fp,mtop,ir,fr);
}
- else if (EEL_FULL(fr->eeltype) || (fr->eeltype == eelSHIFT) ||
- (fr->eeltype == eelUSER) || (fr->eeltype == eelSWITCH))
+ else if (fr->eeltype == eelSHIFT)
{
- /* We must use the long range cut-off for neighboursearching...
- * An extra range of e.g. 0.1 nm (half the size of a charge group)
- * is necessary for neighboursearching. This allows diffusion
- * into the cut-off range (between neighborlist updates),
- * and gives more accurate forces because all atoms within the short-range
- * cut-off rc must be taken into account, while the ns criterium takes
- * only those with the center of geometry within the cut-off.
- * (therefore we have to add half the size of a charge group, plus
- * something to account for diffusion if we have nstlist > 1)
- */
for(m=0; (m<DIM); m++)
box_size[m]=box[m][m];
- if (fr->eeltype == eelPPPM && fr->phi == NULL)
- snew(fr->phi,natoms);
-
- if ((fr->eeltype==eelPPPM) || (fr->eeltype==eelPOISSON) ||
- (fr->eeltype == eelSHIFT && fr->rcoulomb > fr->rcoulomb_switch))
+ if ((fr->eeltype == eelSHIFT && fr->rcoulomb > fr->rcoulomb_switch))
set_shift_consts(fp,fr->rcoulomb_switch,fr->rcoulomb,box_size,fr);
}
+++ /dev/null
-/*
- *
- * 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:
- * GROwing Monsters And Cloning Shrimps
- */
-/* This file is completely threadsafe - keep it that way! */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include "typedefs.h"
-#include "futil.h"
-#include "vec.h"
-#include "physics.h"
-#include "coulomb.h"
-#include "pppm.h"
-#include "xvgr.h"
-#include "gmxfio.h"
-#include "pppm.h"
-#include "smalloc.h"
-
-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
-}
-
-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.
- */
-
-{
- int i,j,k;
- int iip,jjp,kkp;
- real ggg;
-
- /* fprintf(stderr,"Symmetrizing Ghat function\n"); */
- /* Only the lower octant of the rectangle has been saved,
- * so we must construct the other 7 octants by symmetry operations.
- */
- for(i=0; (i<=nx/2); i++) {
- iip = (nx-i) % nx;
- for(j=0; (j<=ny/2); j++) {
- jjp = (ny-j) % ny;
- for(k=0; (k<=nz/2); k++) {
- kkp = (nz-k) % nz;
- ggg = ghat[i][j][k];
- ghat[i] [jjp][k] = ggg;
- ghat[i] [j] [kkp] = ggg;
- ghat[i] [jjp][kkp] = ggg;
- ghat[iip][j] [k] = ggg;
- ghat[iip][jjp][k] = ggg;
- ghat[iip][j] [kkp] = ggg;
- ghat[iip][jjp][kkp] = ggg;
- }
- }
- }
-}
-
-void mk_ghat(FILE *fp,int nx,int ny,int nz,real ***ghat,
- rvec box,real r1,real rc,gmx_bool bSym,gmx_bool bOld)
-{
- int ix,iy,iz;
- int ixmax,iymax,izmax;
- real k2,ggg;
- rvec k,lll;
-
- calc_lll(box,lll);
-
- if (bSym) {
- ixmax=nx/2+1;
- iymax=ny/2+1;
- izmax=nz/2+1;
- }
- else {
- ixmax=nx;
- iymax=ny;
- izmax=nz;
- }
- /* Loop over lattice vectors in fourier space */
- for(ix=0; (ix < ixmax); ix++) {
- for(iy=0; (iy < iymax); iy++) {
- for(iz=0; (iz < izmax); iz++) {
- calc_k(lll,ix,iy,iz,nx,ny,nz,k);
- k2 = iprod(k,k);
- if ((ix == 0) && (iy == 0) && (iz == 0))
- ggg = 0.0;
- else {
- if (bOld)
- ggg = gk(sqrt(k2),rc,r1)/(k2*EPSILON0);
- else
- ggg = gknew(sqrt(k2),rc,r1)/(k2*EPSILON0);
- }
- ghat[ix][iy][iz]=ggg;
- }
- }
- }
- if (bSym)
- symmetrize_ghat(nx,ny,nz,ghat);
-}
-
-void wr_ghat(const char *fn,const output_env_t oenv,
- int n1max,int n2max,int n3max,real h1,real h2,real h3,
- real ***ghat,int nalias,int porder,int niter,gmx_bool bSym,rvec beta,
- real r1,real rc,real pval,real zval,real eref,real qopt)
-{
- FILE *fp;
- int N1MAX,N2MAX,N3MAX;
- gmx_bool bNL=FALSE;
- real rx,ry,rz;
- int ii,jj,kk,nn;
-
- fp = gmx_fio_fopen(fn,"w");
- fprintf(fp,"%8d %8d %8d %15.10e %15.10e %15.10e\n",
- n1max,n2max,n3max,h1,h2,h3);
- fprintf(fp,"%8d %8d %8d %8d %15.10e %15.10e %15.10e\n",
- nalias,porder,niter,bSym,beta[XX],beta[YY],beta[ZZ]);
- fprintf(fp,"%10g %10g %10g %10g %10g %10g\n",
- rc,r1,pval,zval,eref,qopt);
-
- if (bSym) {
- N1MAX = n1max/2+1;
- N2MAX = n2max/2+1;
- N3MAX = n3max/2+1;
- }
- else {
- N1MAX = n1max;
- N2MAX = n2max;
- N3MAX = n3max;
- }
- for(ii=0; (ii<N1MAX); ii++) {
- for(jj=0; (jj<N2MAX); jj++) {
- for(kk=0,nn=1; (kk<N3MAX); kk++,nn++) {
- bNL=FALSE;
- fprintf(fp," %12.5e",ghat[ii][jj][kk]);
- if ((nn % 6) == 0) {
- fprintf(fp,"\n");
- bNL=TRUE;
- }
- }
- if (!bNL)
- fprintf(fp,"\n");
- }
- }
- gmx_fio_fclose(fp);
-
- fp=xvgropen("ghat.xvg","G-Hat","k","gk",oenv);
- for(ii=0; (ii<N1MAX); ii++) {
- rx=sqr((real)(ii*h1));
- for(jj=0; (jj<N2MAX); jj++) {
- ry=rx+sqr((real)(jj*h2));
- for(kk=0; (kk<N3MAX); kk++) {
- rz=ry+sqr((real)(kk*h3));
- fprintf(fp,"%10g %10g\n",sqrt(rz),ghat[ii][jj][kk]);
- }
- }
- }
- gmx_fio_fclose(fp);
-}
-
-void pr_scalar_gk(const char *fn,const output_env_t oenv,int nx,int ny,int nz,
- rvec box,real ***ghat)
-{
- FILE *fp;
- int ii,jj,kk;
- real k1;
- rvec k,lll;
-
- calc_lll(box,lll);
-
- fp=xvgropen(fn,"G-Hat","k","gk",oenv);
- for(ii=0; (ii<nx); ii++) {
- for(jj=0; (jj<ny); jj++) {
- for(kk=0; (kk<nz); kk++) {
- calc_k(lll,ii,jj,kk,nx,ny,nz,k);
- k1 = norm(k);
- fprintf(fp,"%10g %10g\n",k1,ghat[ii][jj][kk]);
- }
- }
- }
- gmx_fio_fclose(fp);
-}
-
-static real ***mk_rgrid(int nx,int ny,int nz)
-{
- real *ptr1;
- real **ptr2;
- real ***ptr3;
- int i,j,n2,n3;
-
- snew(ptr1,nx*ny*nz);
- snew(ptr2,nx*ny);
- snew(ptr3,nx);
-
- n2=n3=0;
- for(i=0; (i<nx); i++) {
- ptr3[i]=&(ptr2[n2]);
- for(j=0; (j<ny); j++,n2++) {
- ptr2[n2] = &(ptr1[n3]);
- n3 += nz;
- }
- }
- return ptr3;
-}
-
-real ***rd_ghat(FILE *log,const output_env_t oenv,char *fn,ivec igrid,
- rvec gridspace, rvec beta,int *porder,real *r1,real *rc)
-{
- FILE *in;
- real ***gh;
- double gx,gy,gz,alX,alY,alZ,ddd;
- double acut,pval,zval,eref,qopt,r11;
- int nalias,niter,bSym;
- int ix,iy,iz,ixmax,iymax,izmax;
-
- in=gmx_fio_fopen(fn,"r");
- if(6 != fscanf(in,"%d%d%d%lf%lf%lf",&ix,&iy,&iz,&gx,&gy,&gz))
- {
- gmx_fatal(FARGS,"Error reading from file %s",fn);
- }
-
-
- igrid[XX]=ix, igrid[YY]=iy, igrid[ZZ]=iz;
- gridspace[XX]=gx, gridspace[YY]=gy, gridspace[ZZ]=gz;
- if(7 != fscanf(in,"%d%d%d%d%lf%lf%lf",&nalias,porder,&niter,&bSym,&alX,&alY,&alZ))
- {
- gmx_fatal(FARGS,"Error reading from file %s",fn);
- }
-
- if(6 != fscanf(in,"%lf%lf%lf%lf%lf%lf",&acut,&r11,&pval,&zval,&eref,&qopt))
- {
- gmx_fatal(FARGS,"Error reading from file %s",fn);
- }
-
- *r1 = r11;
- *rc = acut;
-
- fprintf(log,"\nOpened %s for reading ghat function\n",fn);
- fprintf(log,"gridsize: %10d %10d %10d\n",ix,iy,iz);
- fprintf(log,"spacing: %10g %10g %10g\n",gx,gy,gz);
- fprintf(log," nalias porder niter bSym beta[X-Z]\n"
- "%10d%10d%10d%10d%10g%10g%10g\n",
- nalias,*porder,niter,bSym,alX,alY,alZ);
- fprintf(log," acut r1 pval zval eref qopt\n"
- "%10g%10g%10g%10g%10g%10g\n",acut,*r1,pval,zval,eref,qopt);
- fflush(log);
-
- beta[XX] = alX;
- beta[YY] = alY;
- beta[ZZ] = alZ;
-
- gh = mk_rgrid(ix,iy,iz);
- if (bSym) {
- ixmax=igrid[XX]/2+1;
- iymax=igrid[YY]/2+1;
- izmax=igrid[ZZ]/2+1;
- }
- else {
- ixmax=igrid[XX];
- iymax=igrid[YY];
- izmax=igrid[ZZ];
- }
- fprintf(log,"Reading ghat of %d %d %d\n",ixmax,iymax,izmax);
- for(ix=0; (ix<ixmax); ix++)
- for(iy=0; (iy<iymax); iy++)
- for(iz=0; (iz<izmax); iz++) {
- if( 1 != fscanf(in,"%lf",&ddd))
- {
- gmx_fatal(FARGS,"Error reading from file %s",fn);
- }
-
- gh[ix][iy][iz] = ddd;
- }
- gmx_fio_fclose(in);
-
- wr_ghat("output.hat",oenv,igrid[XX],igrid[YY],igrid[ZZ],gx,gy,gz,gh,
- nalias,*porder,niter,bSym,beta,
- *r1,*rc,pval,zval,eref,qopt);
-
- if (bSym)
- symmetrize_ghat(igrid[XX],igrid[YY],igrid[ZZ],gh);
-
- fprintf(log,"\nSuccessfully read ghat function!\n");
-
-
- return gh;
-}
-
* Use the box from last timestep since we already called update().
*/
- enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres,
- (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
+ enerd->term[F_PRES] = calc_pres(fr->ePBC,ir->nwall,box,ekind->ekin,total_vir,pres);
/* Calculate long range corrections to pressure and energy */
/* this adds to enerd->term[F_PRES] and enerd->term[F_ETOT],
clear_mat(ekin);
enerd->term[F_PRES] =
- calc_pres(fr->ePBC,inputrec->nwall,ems->s.box,ekin,vir,pres,
- (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
+ calc_pres(fr->ePBC,inputrec->nwall,ems->s.box,ekin,vir,pres);
sum_dhdl(enerd,ems->s.lambda,inputrec);
gmx_bool bPPnode; /* Node also does particle-particle forces */
gmx_bool bFEP; /* Compute Free energy contribution */
int nkx,nky,nkz; /* Grid dimensions */
+ gmx_bool bP3M; /* Do P3M: optimize the influence function */
int pme_order;
real epsilon_r;
}
-
-void make_bspline_moduli(splinevec bsp_mod,int nx,int ny,int nz,int order)
+static void make_bspline_moduli(splinevec bsp_mod,
+ int nx,int ny,int nz,int order)
{
int nmax=max(nx,max(ny,nz));
real *data,*ddata,*bsp_data;
sfree(bsp_data);
}
+
+/* Return the P3M optimal influence function */
+static double do_p3m_influence(double z, int order)
+{
+ double z2,z4;
+
+ z2 = z*z;
+ z4 = z2*z2;
+
+ /* The formula and most constants can be found in:
+ * Ballenegger et al., JCTC 8, 936 (2012)
+ */
+ switch(order)
+ {
+ case 2:
+ return 1.0 - 2.0*z2/3.0;
+ break;
+ case 3:
+ return 1.0 - z2 + 2.0*z4/15.0;
+ break;
+ case 4:
+ return 1.0 - 4.0*z2/3.0 + 2.0*z4/5.0 + 4.0*z2*z4/315.0;
+ break;
+ case 5:
+ return 1.0 - 5.0*z2/3.0 + 7.0*z4/9.0 - 17.0*z2*z4/189.0 + 2.0*z4*z4/2835.0;
+ break;
+ case 6:
+ return 1.0 - 2.0*z2 + 19.0*z4/15.0 - 256.0*z2*z4/945.0 + 62.0*z4*z4/4725.0 + 4.0*z2*z4*z4/155925.0;
+ break;
+ case 7:
+ return 1.0 - 7.0*z2/3.0 + 28.0*z4/15.0 - 16.0*z2*z4/27.0 + 26.0*z4*z4/405.0 - 2.0*z2*z4*z4/1485.0 + 4.0*z4*z4*z4/6081075.0;
+ case 8:
+ return 1.0 - 8.0*z2/3.0 + 116.0*z4/45.0 - 344.0*z2*z4/315.0 + 914.0*z4*z4/4725.0 - 248.0*z4*z4*z2/22275.0 + 21844.0*z4*z4*z4/212837625.0 - 8.0*z4*z4*z4*z2/638512875.0;
+ break;
+ }
+
+ return 0.0;
+}
+
+/* Calculate the P3M B-spline moduli for one dimension */
+static void make_p3m_bspline_moduli_dim(real *bsp_mod,int n,int order)
+{
+ double zarg,zai,sinzai,infl;
+ int maxk,i;
+
+ if (order > 8)
+ {
+ gmx_fatal(FARGS,"The current P3M code only supports orders up to 8");
+ }
+
+ zarg = M_PI/n;
+
+ maxk = (n + 1)/2;
+
+ for(i=-maxk; i<0; i++)
+ {
+ zai = zarg*i;
+ sinzai = sin(zai);
+ infl = do_p3m_influence(sinzai,order);
+ bsp_mod[n+i] = infl*infl*pow(sinzai/zai,-2.0*order);
+ }
+ bsp_mod[0] = 1.0;
+ for(i=1; i<maxk; i++)
+ {
+ zai = zarg*i;
+ sinzai = sin(zai);
+ infl = do_p3m_influence(sinzai,order);
+ bsp_mod[i] = infl*infl*pow(sinzai/zai,-2.0*order);
+ }
+}
+
+/* Calculate the P3M B-spline moduli */
+static void make_p3m_bspline_moduli(splinevec bsp_mod,
+ int nx,int ny,int nz,int order)
+{
+ make_p3m_bspline_moduli_dim(bsp_mod[XX],nx,order);
+ make_p3m_bspline_moduli_dim(bsp_mod[YY],ny,order);
+ make_p3m_bspline_moduli_dim(bsp_mod[ZZ],nz,order);
+}
+
+
static void setup_coordinate_communication(pme_atomcomm_t *atc)
{
int nslab,n,i;
pme->nkx = ir->nkx;
pme->nky = ir->nky;
pme->nkz = ir->nkz;
+ pme->bP3M = (ir->coulombtype == eelP3M_AD || getenv("GMX_PME_P3M") != NULL);
pme->pme_order = ir->pme_order;
pme->epsilon_r = ir->epsilon_r;
pme->cfftgridB = NULL;
}
- make_bspline_moduli(pme->bsp_mod,pme->nkx,pme->nky,pme->nkz,pme->pme_order);
+ if (!pme->bP3M)
+ {
+ /* Use plain SPME B-spline interpolation */
+ make_bspline_moduli(pme->bsp_mod,pme->nkx,pme->nky,pme->nkz,pme->pme_order);
+ }
+ else
+ {
+ /* Use the P3M grid-optimized influence function */
+ make_p3m_bspline_moduli(pme->bsp_mod,pme->nkx,pme->nky,pme->nkz,pme->pme_order);
+ }
/* Use atc[0] for spreading */
init_atomcomm(pme,&pme->atc[0],cr,nnodes_major > 1 ? 0 : 1,TRUE);
}
}
}
+ else
+ {
+ *energy = 0;
+ }
if (debug)
{
+++ /dev/null
-/*
- *
- * 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:
- * GROwing Monsters And Cloning Shrimps
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include <math.h>
-#include "physics.h"
-#include "typedefs.h"
-#include "smalloc.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "gmx_fatal.h"
-#include "txtdump.h"
-#include "network.h"
-#include "nrnb.h"
-#include "pppm.h"
-#include "coulomb.h"
-#include "mdrun.h"
-#include "gmx_fft.h"
-#include "pme.h"
-
-#ifdef GMX_MPI
-#include "gmx_parallel_3dfft.h"
-#endif
-
-#define llim (-1)
-#define ulim (1)
-#define llim2 (-3)
-#define ulim2 (3)
-
-
-
-/* PPPM temporarily disabled while working on 2D PME */
-#define DISABLE_PPPM
-
-
-
-#ifndef DISABLE_PPPM
-
-/* TODO: fix thread-safety */
-
-static void calc_invh(rvec box,int nx,int ny,int nz,rvec invh)
-{
- invh[XX] = nx/box[XX];
- invh[YY] = ny/box[YY];
- invh[ZZ] = nz/box[ZZ];
-}
-
-static void calc_weights(int iatom,int nx,int ny,int nz,
- rvec x,rvec box,rvec invh,ivec ixyz,real WXYZ[])
-{
- const real half=0.5;
- tensor wxyz;
- real abc,ttt,fact;
-#ifdef DEBUG
- real wtot;
-#endif
- ivec nxyz;
- int it,j,k,m,nm;
- real Wx,Wy,Wzx,Wzy,Wzz;
-
- fact = 0.125;
-
- nxyz[XX] = nx;
- nxyz[YY] = ny;
- nxyz[ZZ] = nz;
- for(m=0; (m<DIM); m++) {
- /* Put particle in the box... */
- ttt = x[m]*invh[m];
-
- /* Round to nearest grid point. Do the math in integer! */
- it = (ttt+0.5);
- nm = nxyz[m];
- if (it < 0) {
- ttt+=nm;
- it +=nm;
- }
- else if (it >= nm) {
- ttt-=nm;
- it -=nm;
- }
- if ((it < 0) || (it >= nxyz[m]))
- gmx_fatal(FARGS,"iatom = %d, it = %d, x=%f, ttt=%f",iatom,it,x[m],ttt);
- ixyz[m] = it;
-
- /* Fraction (offset) from grid point */
- abc = ttt - (real)ixyz[m];
-
- wxyz[m][0] = sqr((real)(half - abc));
- wxyz[m][1] = 1.5 - 2.0*sqr(abc);
- wxyz[m][2] = sqr((real)(half + abc));
- }
- Wzx=wxyz[ZZ][XX];
- Wzy=wxyz[ZZ][YY];
- Wzz=wxyz[ZZ][ZZ];
- for(j=m=0; (j<DIM); j++) {
- Wx = wxyz[XX][j]*fact;
- for(k=0; (k<DIM); k++,m+=3) {
- Wy = Wx*wxyz[YY][k];
- WXYZ[m] = Wy*Wzx;
- WXYZ[m+1] = Wy*Wzy;
- WXYZ[m+2] = Wy*Wzz;
- }
- }
-#ifdef DEBUG
- wtot = 0;
- for(j=0; (j<27); j++)
- wtot+=WXYZ[j];
- fprintf(stderr,"wtot = %g\n",wtot);
-#endif
-#ifdef HACK
- for(j=0; (j<27); j++)
- WXYZ[j] = 0;
- WXYZ[13] = 1.0;
-#endif
-}
-
-static void calc_nxyz(int nx,int ny,int nz,
- int **nnx,int **nny,int **nnz)
-{
- int i;
-
- snew(*nnx,3*nx);
- snew(*nny,3*ny);
- snew(*nnz,3*nz);
- for(i=0; (i<3*nx); i++)
- (*nnx)[i] = i % nx;
- for(i=0; (i<3*ny); i++)
- (*nny)[i] = i % ny;
- for(i=0; (i<3*nz); i++)
- (*nnz)[i] = i % nz;
-}
-
-static void spread_q(FILE *log,gmx_bool bVerbose,
- int start,int nr,
- rvec x[],real charge[],rvec box,
- t_fftgrid *grid,t_nrnb *nrnb)
-{
- static gmx_bool bFirst = TRUE;
- static int *nnx,*nny,*nnz;
- rvec invh;
- real qi,qwt;
-#ifdef DEBUG
- real qt;
-#endif
- real WXYZ[27];
- ivec ixyz;
- int i,iX,iY,iZ,index;
- int jx,jy,jz,jcx,jcy,jcz;
- int nxyz;
- int nx,ny,nz,nx2,ny2,nz2,la2,la12;
- real *ptr;
-
- unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,&la2,&la12,TRUE,&ptr);
-
- calc_invh(box,nx,ny,nz,invh);
-
- if (bFirst) {
- if (log) {
- fprintf(log,
- "Spreading Charges using Triangle Shaped on %dx%dx%d grid\n",
- nx,ny,nz);
- fprintf(log,"invh = %10g,%10g,%10g\n",invh[XX],invh[YY],invh[ZZ]);
- }
-
- calc_nxyz(nx,ny,nz,&nnx,&nny,&nnz);
-
- bFirst = FALSE;
- }
-
- for(i=start; (i<start+nr); i++) {
- qi=charge[i];
-
- /* Each charge is spread over the nearest 27 grid cells,
- * thus we loop over -1..1 in X,Y and Z direction
- * We apply the TSC (triangle shaped charge)
- * see Luty et. al, JCP 103 (1995) 3014
- */
-
- if (qi != 0.0) {
- calc_weights(i,nx,ny,nz,x[i],box,invh,ixyz,WXYZ);
- iX = ixyz[XX] + nx;
- iY = ixyz[YY] + ny;
- iZ = ixyz[ZZ] + nz;
-
-#ifdef DEBUG
- qt=0;
-#endif
- nxyz = 0;
- for(jx=-1; (jx<=1); jx++) {
- jcx = nnx[iX + jx];
- for(jy=-1; (jy<=1); jy++) {
- jcy = nny[iY + jy];
- for(jz=-1; (jz<=1); jz++,nxyz++) {
- jcz = nnz[iZ + jz];
- index = INDEX(jcx,jcy,jcz);
- qwt = qi*WXYZ[nxyz];
- grid->ptr[index]+=qwt;
-#ifdef DEBUG
- qt += qwt;
- if (bVerbose && log)
- fprintf(log,"spread %4d %2d %2d %2d %10.3e, weight=%10.3e\n",
- index,jcx,jcy,jcz,grid->ptr[index],WXYZ[nxyz]);
-#endif
- }
- }
- }
-#ifdef DEBUG
- if (log)
- fprintf(log,"q[%3d] = %6.3f, qwt = %6.3f\n",i,qi,qt);
-#endif
- }
- }
- inc_nrnb(nrnb,eNR_SPREADQ,9*nr);
- inc_nrnb(nrnb,eNR_WEIGHTS,3*nr);
-}
-
-static real gather_inner(int JCXYZ[],real WXYZ[],int ixw[],int iyw[],int izw[],
- int la2,int la12,
- real c1x,real c1y,real c1z,real c2x,real c2y,real c2z,
- real qi,rvec f,real ptr[])
-{
- real pi,fX,fY,fZ,weight;
- int jxyz,m,jcx,jcy,jcz;
- int jcx0,jcy0,jcz0;
-
- pi = 0.0;
- fX = 0.0;
- fY = 0.0;
- fZ = 0.0;
-
- /* Now loop over 27 surrounding vectors */
- for(jxyz=m=0; (jxyz < 27); jxyz++,m+=3) {
- jcx = JCXYZ[m];
- jcy = JCXYZ[m+1];
- jcz = JCXYZ[m+2];
- weight = WXYZ[jxyz];
-
- jcx0 = ixw[jcx];
- jcy0 = iyw[jcy];
- jcz0 = izw[jcz];
-
- /* Electrostatic Potential ! */
- pi += weight * ptr[INDEX(jcx0,jcy0,jcz0)];
-
- /* Forces */
- fX += weight * ((c1x*(ptr[INDEX(ixw[jcx-1],jcy0,jcz0)] -
- ptr[INDEX(ixw[jcx+1],jcy0,jcz0)] )) +
- (c2x*(ptr[INDEX(ixw[jcx-2],jcy0,jcz0)] -
- ptr[INDEX(ixw[jcx+2],jcy0,jcz0)] )));
- fY += weight * ((c1y*(ptr[INDEX(jcx0,iyw[jcy-1],jcz0)] -
- ptr[INDEX(jcx0,iyw[jcy+1],jcz0)] )) +
- (c2y*(ptr[INDEX(jcx0,iyw[jcy-2],jcz0)] -
- ptr[INDEX(jcx0,iyw[jcy+2],jcz0)] )));
- fZ += weight * ((c1z*(ptr[INDEX(jcx0,jcy0,izw[jcz-1])] -
- ptr[INDEX(jcx0,jcy0,izw[jcz+1])] )) +
- (c2z*(ptr[INDEX(jcx0,jcy0,izw[jcz-2])] -
- ptr[INDEX(jcx0,jcy0,izw[jcz+2])] )));
- }
- f[XX] += qi*fX;
- f[YY] += qi*fY;
- f[ZZ] += qi*fZ;
-
- return pi;
-}
-
-static real gather_f(FILE *log,gmx_bool bVerbose,
- int start,int nr,rvec x[],rvec f[],real charge[],rvec box,
- real pot[],t_fftgrid *grid,rvec beta,t_nrnb *nrnb)
-{
- static gmx_bool bFirst=TRUE;
- static int *nnx,*nny,*nnz;
- static int JCXYZ[81];
- int i,m;
- real energy;
- real qi,pi;
- ivec ixyz;
- rvec invh,c1,c2;
- real WXYZ[27];
- real c1x,c1y,c1z,c2x,c2y,c2z;
- int ixw[7],iyw[7],izw[7];
- int ll;
- int nx,ny,nz,nx2,ny2,nz2,la2,la12;
- real *ptr;
-
- unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,&la2,&la12,TRUE,&ptr);
-
- calc_invh(box,nx,ny,nz,invh);
-
- for(m=0; (m<DIM); m++) {
- c1[m] = (beta[m]/2.0)*invh[m];
- c2[m] = ((1.0-beta[m])/4.0)*invh[m];
- }
- c1x = c1[XX];
- c1y = c1[YY];
- c1z = c1[ZZ];
- c2x = c2[XX];
- c2y = c2[YY];
- c2z = c2[ZZ];
-
- if (bFirst) {
- if (log) {
- fprintf(log,"Gathering Forces using Triangle Shaped on %dx%dx%d grid\n",
- nx,ny,nz);
- fprintf(log,"beta = %10g,%10g,%10g\n",beta[XX],beta[YY],beta[ZZ]);
- fprintf(log,"c1 = %10g,%10g,%10g\n",c1[XX],c1[YY],c1[ZZ]);
- fprintf(log,"c2 = %10g,%10g,%10g\n",c2[XX],c2[YY],c2[ZZ]);
- fprintf(log,"invh = %10g,%10g,%10g\n",invh[XX],invh[YY],invh[ZZ]);
- }
- calc_nxyz(nx,ny,nz,&nnx,&nny,&nnz);
-
- for(i=0; (i<27); i++) {
- JCXYZ[3*i] = 2 + (i/9);
- JCXYZ[3*i+1] = 2 + (i/3) % 3;
- JCXYZ[3*i+2] = 2 + (i % 3);
- }
-
- bFirst = FALSE;
- }
-
- energy=0.0;
- for(i=start; (i<start+nr); i++) {
- /* Each charge is spread over the nearest 27 grid cells,
- * thus we loop over -1..1 in X,Y and Z direction
- * We apply the TSC (triangle shaped charge)
- * see Luty et. al, JCP 103 (1995) 3014
- */
-
- calc_weights(i,nx,ny,nz,x[i],box,invh,ixyz,WXYZ);
-
- for(ll=llim2; (ll<=ulim2); ll++) {
- ixw[ll-llim2] = nnx[ixyz[XX]+ll+nx];
- iyw[ll-llim2] = nny[ixyz[YY]+ll+ny];
- izw[ll-llim2] = nnz[ixyz[ZZ]+ll+nz];
- }
-
- qi = charge[i];
- pi = gather_inner(JCXYZ,WXYZ,ixw,iyw,izw,la2,la12,
- c1x,c1y,c1z,c2x,c2y,c2z,
- qi,f[i],ptr);
-
- energy += pi*qi;
- pot[i] = pi;
- }
-
- inc_nrnb(nrnb,eNR_GATHERF,27*nr);
- inc_nrnb(nrnb,eNR_WEIGHTS,3*nr);
-
- return energy*0.5;
-}
-
-static void convolution(FILE *fp,gmx_bool bVerbose,t_fftgrid *grid,real ***ghat,
- t_commrec *cr)
-{
- int i,j,k,index;
- real gk;
- int nx,ny,nz,nx2,ny2,nz2,la2,la12;
- t_complex *ptr;
- int *nTest;
- int jstart,jend;
-
- unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,
- &la2,&la12,FALSE,(real **)&ptr);
- snew(nTest,grid->nptr);
-
- if(PAR(cr)) {
-#if (defined GMX_MPI && !defined GMX_WITHOUT_FFTW)
- jstart=grid->pfft.local_y_start_after_transpose;
- jend=jstart+grid->pfft.local_ny_after_transpose;
-
- for(j=jstart; (j<jend); j++) { /* local cells */
- for(i=0; (i<nx); i++) {
- for(k=0;k<(nz/2+1); k++) {
- gk = ghat[i][j][k];
- index = INDEX(j,i,k);
- ptr[index].re *= gk;
- ptr[index].im *= gk;
- nTest[index]++;
- }
- }
- }
-#ifdef DEBUG
- for(j=jstart; (j<jend); j++) {
- for(i=0; (i<nx); i++) {
- for(k=0; k<(nz/2+1); k++) {
- index = INDEX(j,i,k);
- if (nTest[index] != 1)
- fprintf(fp,"Index %d sucks, set %d times\n",index,nTest[index]);
- }
- }
- }
-#endif /* DEBUG */
-#endif /* GMX_MPI */
- } else { /* if not running in parallel */
- for(i=0; (i<nx); i++) {
- for(j=0; (j<ny); j++) {
- for(k=0;k<(nz/2+1); k++) {
- gk = ghat[i][j][k];
- index = INDEX(i,j,k);
- ptr[index].re *= gk;
- ptr[index].im *= gk;
- nTest[index]++;
- }
- }
- }
-#ifdef DEBUG
- for(i=0; (i<nx); i++) {
- for(j=0; (j<ny); j++) {
- for(k=0; k<(nz/2+1); k++) {
- index = INDEX(i,j,k);
- if (nTest[index] != 1)
- fprintf(fp,"Index %d sucks, set %d times\n",index,nTest[index]);
- }
- }
- }
-#endif
- }
- sfree(nTest);
-}
-
-
-void solve_pppm(FILE *fp,t_commrec *cr,
- t_fftgrid *grid,real ***ghat,rvec box,
- gmx_bool bVerbose,t_nrnb *nrnb)
-{
- int ntot,npppm;
-
- /* if (bVerbose)
- print_fftgrid(fp,"Q-Real",grid,grid->nxyz,"qreal.pdb",box,TRUE);*/
-
- gmxfft3D(grid,GMX_FFT_REAL_TO_COMPLEX,cr);
-
- /* if (bVerbose) {
- print_fftgrid(fp,"Q-k",grid,1.0,"qk-re.pdb",box,TRUE);
- print_fftgrid(fp,"Q-k",grid,1.0,"qk-im.pdb",box,FALSE);
- fprintf(stderr,"Doing convolution\n");
- }*/
-
- convolution(fp,bVerbose,grid,ghat,cr);
-
- /* if (bVerbose)
- print_fftgrid(fp,"Convolution",grid,1.0,
- "convolute.pdb",box,TRUE);*/
-
- gmxfft3D(grid,GMX_FFT_COMPLEX_TO_REAL,cr);
-
- /* if (bVerbose)
- print_fftgrid(fp,"Potential",grid,1.0,"pot.pdb",box,TRUE);*/
-
- ntot = grid->nxyz;
- npppm = ntot*log((real)ntot)/log(2.0);
- inc_nrnb(nrnb,eNR_FFT,2*npppm);
- inc_nrnb(nrnb,eNR_CONV,ntot);
-}
-
-
-static rvec beta;
-static real ***ghat=NULL;
-static t_fftgrid *grid=NULL;
-
-#endif
-
-
-int gmx_pppm_init(FILE *log, t_commrec *cr,
- const output_env_t oenv, gmx_bool bVerbose,
- gmx_bool bOld, matrix box,
- char *ghatfn, t_inputrec *ir,
- gmx_bool bReproducible)
-{
- int nx,ny,nz,m,porder;
- ivec grids;
- real r1,rc;
- const real tol = 1e-5;
- rvec box_diag,spacing;
-
-#ifdef DISABLE_PPPM
- gmx_fatal(FARGS,"PPPM is not functional in the current version, we plan to implement PPPM through a small modification of the PME code.");
- return -1;
-#else
-
-#ifdef GMX_WITHOUT_FFTW
- gmx_fatal(FARGS,"PPPM used, but GROMACS was compiled without FFTW support!\n");
-#endif
-
- if (log) {
- if (cr != NULL) {
- if (cr->nnodes > 1)
- fprintf(log,"Initializing parallel PPPM.\n");
- }
- fprintf(log,
- "Will use the PPPM algorithm for long-range electrostatics\n");
- }
-
- for(m=0; m<DIM; m++)
- box_diag[m] = box[m][m];
-
- if (!gmx_fexist(ghatfn)) {
- beta[XX]=beta[YY]=beta[ZZ]= 1.85;
- nx = ir->nkx;
- ny = ir->nky;
- nz = ir->nkz;
-
- if (log) {
- fprintf(log,"Generating Ghat function\n");
- fprintf(log,"Grid size is %d x %d x %d\n",nx,ny,nz);
- }
-
- if ((nx < 4) || (ny < 4) || (nz < 4))
- gmx_fatal(FARGS,"Grid must be at least 4 points in all directions");
-
- ghat = mk_rgrid(nx,ny,nz);
- mk_ghat(NULL,nx,ny,nz,ghat,box_diag,
- ir->rcoulomb_switch,ir->rcoulomb,TRUE,bOld);
-
- if (bVerbose)
- pr_scalar_gk("generghat.xvg",oenv,nx,ny,nz,box_diag,ghat);
- }
- else {
- fprintf(stderr,"Reading Ghat function from %s\n",ghatfn);
- ghat = rd_ghat(log,oenv,ghatfn,grids,spacing,beta,&porder,&r1,&rc);
-
- /* Check whether cut-offs correspond */
- if ((fabs(r1-ir->rcoulomb_switch)>tol) || (fabs(rc-ir->rcoulomb)>tol)) {
- if (log) {
- fprintf(log,"rcoulomb_switch = %10.3e rcoulomb = %10.3e"
- " r1 = %10.3e rc = %10.3e\n",
- ir->rcoulomb_switch,ir->rcoulomb,r1,rc);
- fflush(log);
- }
- gmx_fatal(FARGS,"Cut-off lengths in tpb file and Ghat file %s "
- "do not match\nCheck your log file!",ghatfn);
- }
-
- /* Check whether boxes correspond */
- for(m=0; (m<DIM); m++)
- if (fabs(box_diag[m]-grids[m]*spacing[m]) > tol) {
- if (log) {
- pr_rvec(log,0,"box",box_diag,DIM,TRUE);
- pr_rvec(log,0,"grid-spacing",spacing,DIM,TRUE);
- pr_ivec(log,0,"grid size",grids,DIM,TRUE);
- fflush(log);
- }
- gmx_fatal(FARGS,"Box sizes in tpb file and Ghat file %s do not match\n"
- "Check your log file!",ghatfn);
- }
-
- if (porder != 2)
- gmx_fatal(FARGS,"porder = %d, should be 2 in %s",porder,ghatfn);
-
- nx = grids[XX];
- ny = grids[YY];
- nz = grids[ZZ];
-
- if (bVerbose)
- pr_scalar_gk("optimghat.xvg",oenv,nx,ny,nz,box_diag,ghat);
- }
- /* Now setup the FFT things */
-#ifdef GMX_MPI
- cr->mpi_comm_mygroup=cr->mpi_comm_mysim;
-#endif
- grid = mk_fftgrid(nx,ny,nz,NULL,NULL,cr,bReproducible);
-
- return 0;
-#endif
-}
-
-int gmx_pppm_do(FILE *log, gmx_pme_t pme,
- gmx_bool bVerbose,
- rvec x[], rvec f[],
- real charge[], rvec box,
- real phi[], t_commrec *cr,
- int start, int nr,
- t_nrnb *nrnb,
- int pme_order, real *energy)
-{
-#ifdef DISABLE_PPPM
- gmx_fatal(FARGS,"PPPM temporarily disabled while working on 2DPME\n");
- return -1;
-#else
-
- /* Make the grid empty */
- clear_fftgrid(grid);
-
- /* First step: spreading the charges over the grid. */
- spread_q(log,bVerbose,start,nr,x,charge,box,grid,nrnb);
-
- /* In the parallel code we have to sum the grids from neighbouring nodes */
- if (PAR(cr))
- gmx_sum_qgrid(pme,cr,grid,GMX_SUM_QGRID_FORWARD);
-
- /* Second step: solving the poisson equation in Fourier space */
- solve_pppm(log,cr,grid,ghat,box,bVerbose,nrnb);
-
- /* In the parallel code we have to sum once again... */
- if (PAR(cr))
- gmx_sum_qgrid(pme,cr,grid,GMX_SUM_QGRID_BACKWARD);
-
- /* Third and last step: gather the forces, energies and potential
- * from the grid.
- */
- *energy = gather_f(log,bVerbose,start,nr,x,f,charge,box,
- phi,grid,beta,nrnb);
-
- return 0;
-#endif
-}
-
-#ifndef DISABLE_PPPM
-static int gmx_pppm_opt_do(FILE *log, gmx_pme_t pme,
- t_inputrec *ir, gmx_bool bVerbose,
- int natoms,
- rvec x[], rvec f[],
- real charge[], rvec box,
- real phi[], t_commrec *cr,
- t_nrnb *nrnb, rvec beta,
- t_fftgrid *grid, gmx_bool bOld,
- real *energy)
-{
- real ***ghat;
- int nx,ny,nz;
-
- if (log)
- fprintf(log,"Generating Ghat function\n");
- nx = ir->nkx;
- ny = ir->nky;
- nz = ir->nkz;
- ghat = mk_rgrid(nx,ny,nz);
- mk_ghat(NULL,nx,ny,nz,ghat,box,ir->rcoulomb_switch,ir->rcoulomb,TRUE,bOld);
-
- /* pr_scalar_gk("generghat.xvg",nx,ny,nz,box,ghat); */
-
- /* Now start the actual PPPM procedure.
- * First step: spreading the charges over the grid.
- */
- /* Make the grid empty */
- clear_fftgrid(grid);
-
- spread_q(log,bVerbose,0,natoms,x,charge,box,grid,nrnb);
-
- /* Second step: solving the poisson equation in Fourier space */
- solve_pppm(log,cr,grid,ghat,box,bVerbose,nrnb);
-
- /* Third and last step: gather the forces, energies and potential
- * from the grid.
- */
- *energy = gather_f(log,bVerbose,0,natoms,x,f,charge,box,phi,grid,beta,nrnb);
-
- free_rgrid(ghat,nx,ny);
-
- return 0;
-}
-
-#endif
#include "force.h"
#include "bondf.h"
#include "pme.h"
-#include "pppm.h"
#include "disre.h"
#include "orires.h"
#include "network.h"
case eelCUT:
tabsel[etiCOUL] = etabCOUL;
break;
- case eelPPPM:
case eelPOISSON:
tabsel[etiCOUL] = etabShift;
break;
break;
case eelEWALD:
case eelPME:
+ case eelP3M_AD:
tabsel[etiCOUL] = etabEwald;
break;
case eelPMESWITCH:
#include "disre.h"
#include "orires.h"
#include "dihre.h"
-#include "pppm.h"
#include "pme.h"
#include "mdatoms.h"
#include "qmmm.h"
}
}
- /* Initiate PPPM if necessary */
- if (fr->eeltype == eelPPPM)
- {
- if (mdatoms->nChargePerturbed)
- {
- gmx_fatal(FARGS,"Free energy with %s is not implemented",
- eel_names[fr->eeltype]);
- }
- status = gmx_pppm_init(fplog,cr,oenv,FALSE,TRUE,box,
- getenv("GMXGHAT"),inputrec, (Flags & MD_REPRODUCIBLE));
- if (status != 0)
- {
- gmx_fatal(FARGS,"Error %d initializing PPPM",status);
- }
- }
-
if (EEL_PME(fr->eeltype))
{
ewaldcoeff = fr->ewaldcoeff;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff