3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
46 #include "gmx_fatal.h"
57 #include "gmx_parallel_3dfft.h"
67 /* PPPM temporarily disabled while working on 2D PME */
74 /* TODO: fix thread-safety */
76 static void calc_invh(rvec box,int nx,int ny,int nz,rvec invh)
78 invh[XX] = nx/box[XX];
79 invh[YY] = ny/box[YY];
80 invh[ZZ] = nz/box[ZZ];
83 static void calc_weights(int iatom,int nx,int ny,int nz,
84 rvec x,rvec box,rvec invh,ivec ixyz,real WXYZ[])
94 real Wx,Wy,Wzx,Wzy,Wzz;
101 for(m=0; (m<DIM); m++) {
102 /* Put particle in the box... */
105 /* Round to nearest grid point. Do the math in integer! */
116 if ((it < 0) || (it >= nxyz[m]))
117 gmx_fatal(FARGS,"iatom = %d, it = %d, x=%f, ttt=%f",iatom,it,x[m],ttt);
120 /* Fraction (offset) from grid point */
121 abc = ttt - (real)ixyz[m];
123 wxyz[m][0] = sqr((real)(half - abc));
124 wxyz[m][1] = 1.5 - 2.0*sqr(abc);
125 wxyz[m][2] = sqr((real)(half + abc));
130 for(j=m=0; (j<DIM); j++) {
131 Wx = wxyz[XX][j]*fact;
132 for(k=0; (k<DIM); k++,m+=3) {
141 for(j=0; (j<27); j++)
143 fprintf(stderr,"wtot = %g\n",wtot);
146 for(j=0; (j<27); j++)
152 static void calc_nxyz(int nx,int ny,int nz,
153 int **nnx,int **nny,int **nnz)
160 for(i=0; (i<3*nx); i++)
162 for(i=0; (i<3*ny); i++)
164 for(i=0; (i<3*nz); i++)
168 static void spread_q(FILE *log,gmx_bool bVerbose,
170 rvec x[],real charge[],rvec box,
171 t_fftgrid *grid,t_nrnb *nrnb)
173 static gmx_bool bFirst = TRUE;
174 static int *nnx,*nny,*nnz;
182 int i,iX,iY,iZ,index;
183 int jx,jy,jz,jcx,jcy,jcz;
185 int nx,ny,nz,nx2,ny2,nz2,la2,la12;
188 unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,&la2,&la12,TRUE,&ptr);
190 calc_invh(box,nx,ny,nz,invh);
195 "Spreading Charges using Triangle Shaped on %dx%dx%d grid\n",
197 fprintf(log,"invh = %10g,%10g,%10g\n",invh[XX],invh[YY],invh[ZZ]);
200 calc_nxyz(nx,ny,nz,&nnx,&nny,&nnz);
205 for(i=start; (i<start+nr); i++) {
208 /* Each charge is spread over the nearest 27 grid cells,
209 * thus we loop over -1..1 in X,Y and Z direction
210 * We apply the TSC (triangle shaped charge)
211 * see Luty et. al, JCP 103 (1995) 3014
215 calc_weights(i,nx,ny,nz,x[i],box,invh,ixyz,WXYZ);
224 for(jx=-1; (jx<=1); jx++) {
226 for(jy=-1; (jy<=1); jy++) {
228 for(jz=-1; (jz<=1); jz++,nxyz++) {
230 index = INDEX(jcx,jcy,jcz);
232 grid->ptr[index]+=qwt;
236 fprintf(log,"spread %4d %2d %2d %2d %10.3e, weight=%10.3e\n",
237 index,jcx,jcy,jcz,grid->ptr[index],WXYZ[nxyz]);
244 fprintf(log,"q[%3d] = %6.3f, qwt = %6.3f\n",i,qi,qt);
248 inc_nrnb(nrnb,eNR_SPREADQ,9*nr);
249 inc_nrnb(nrnb,eNR_WEIGHTS,3*nr);
252 static real gather_inner(int JCXYZ[],real WXYZ[],int ixw[],int iyw[],int izw[],
254 real c1x,real c1y,real c1z,real c2x,real c2y,real c2z,
255 real qi,rvec f,real ptr[])
257 real pi,fX,fY,fZ,weight;
258 int jxyz,m,jcx,jcy,jcz;
266 /* Now loop over 27 surrounding vectors */
267 for(jxyz=m=0; (jxyz < 27); jxyz++,m+=3) {
277 /* Electrostatic Potential ! */
278 pi += weight * ptr[INDEX(jcx0,jcy0,jcz0)];
281 fX += weight * ((c1x*(ptr[INDEX(ixw[jcx-1],jcy0,jcz0)] -
282 ptr[INDEX(ixw[jcx+1],jcy0,jcz0)] )) +
283 (c2x*(ptr[INDEX(ixw[jcx-2],jcy0,jcz0)] -
284 ptr[INDEX(ixw[jcx+2],jcy0,jcz0)] )));
285 fY += weight * ((c1y*(ptr[INDEX(jcx0,iyw[jcy-1],jcz0)] -
286 ptr[INDEX(jcx0,iyw[jcy+1],jcz0)] )) +
287 (c2y*(ptr[INDEX(jcx0,iyw[jcy-2],jcz0)] -
288 ptr[INDEX(jcx0,iyw[jcy+2],jcz0)] )));
289 fZ += weight * ((c1z*(ptr[INDEX(jcx0,jcy0,izw[jcz-1])] -
290 ptr[INDEX(jcx0,jcy0,izw[jcz+1])] )) +
291 (c2z*(ptr[INDEX(jcx0,jcy0,izw[jcz-2])] -
292 ptr[INDEX(jcx0,jcy0,izw[jcz+2])] )));
301 static real gather_f(FILE *log,gmx_bool bVerbose,
302 int start,int nr,rvec x[],rvec f[],real charge[],rvec box,
303 real pot[],t_fftgrid *grid,rvec beta,t_nrnb *nrnb)
305 static gmx_bool bFirst=TRUE;
306 static int *nnx,*nny,*nnz;
307 static int JCXYZ[81];
314 real c1x,c1y,c1z,c2x,c2y,c2z;
315 int ixw[7],iyw[7],izw[7];
317 int nx,ny,nz,nx2,ny2,nz2,la2,la12;
320 unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,&la2,&la12,TRUE,&ptr);
322 calc_invh(box,nx,ny,nz,invh);
324 for(m=0; (m<DIM); m++) {
325 c1[m] = (beta[m]/2.0)*invh[m];
326 c2[m] = ((1.0-beta[m])/4.0)*invh[m];
337 fprintf(log,"Gathering Forces using Triangle Shaped on %dx%dx%d grid\n",
339 fprintf(log,"beta = %10g,%10g,%10g\n",beta[XX],beta[YY],beta[ZZ]);
340 fprintf(log,"c1 = %10g,%10g,%10g\n",c1[XX],c1[YY],c1[ZZ]);
341 fprintf(log,"c2 = %10g,%10g,%10g\n",c2[XX],c2[YY],c2[ZZ]);
342 fprintf(log,"invh = %10g,%10g,%10g\n",invh[XX],invh[YY],invh[ZZ]);
344 calc_nxyz(nx,ny,nz,&nnx,&nny,&nnz);
346 for(i=0; (i<27); i++) {
347 JCXYZ[3*i] = 2 + (i/9);
348 JCXYZ[3*i+1] = 2 + (i/3) % 3;
349 JCXYZ[3*i+2] = 2 + (i % 3);
356 for(i=start; (i<start+nr); i++) {
357 /* Each charge is spread over the nearest 27 grid cells,
358 * thus we loop over -1..1 in X,Y and Z direction
359 * We apply the TSC (triangle shaped charge)
360 * see Luty et. al, JCP 103 (1995) 3014
363 calc_weights(i,nx,ny,nz,x[i],box,invh,ixyz,WXYZ);
365 for(ll=llim2; (ll<=ulim2); ll++) {
366 ixw[ll-llim2] = nnx[ixyz[XX]+ll+nx];
367 iyw[ll-llim2] = nny[ixyz[YY]+ll+ny];
368 izw[ll-llim2] = nnz[ixyz[ZZ]+ll+nz];
372 pi = gather_inner(JCXYZ,WXYZ,ixw,iyw,izw,la2,la12,
373 c1x,c1y,c1z,c2x,c2y,c2z,
380 inc_nrnb(nrnb,eNR_GATHERF,27*nr);
381 inc_nrnb(nrnb,eNR_WEIGHTS,3*nr);
386 static void convolution(FILE *fp,gmx_bool bVerbose,t_fftgrid *grid,real ***ghat,
391 int nx,ny,nz,nx2,ny2,nz2,la2,la12;
396 unpack_fftgrid(grid,&nx,&ny,&nz,&nx2,&ny2,&nz2,
397 &la2,&la12,FALSE,(real **)&ptr);
398 snew(nTest,grid->nptr);
401 #if (defined GMX_MPI && !defined GMX_WITHOUT_FFTW)
402 jstart=grid->pfft.local_y_start_after_transpose;
403 jend=jstart+grid->pfft.local_ny_after_transpose;
405 for(j=jstart; (j<jend); j++) { /* local cells */
406 for(i=0; (i<nx); i++) {
407 for(k=0;k<(nz/2+1); k++) {
409 index = INDEX(j,i,k);
417 for(j=jstart; (j<jend); j++) {
418 for(i=0; (i<nx); i++) {
419 for(k=0; k<(nz/2+1); k++) {
420 index = INDEX(j,i,k);
421 if (nTest[index] != 1)
422 fprintf(fp,"Index %d sucks, set %d times\n",index,nTest[index]);
428 } else { /* if not running in parallel */
429 for(i=0; (i<nx); i++) {
430 for(j=0; (j<ny); j++) {
431 for(k=0;k<(nz/2+1); k++) {
433 index = INDEX(i,j,k);
441 for(i=0; (i<nx); i++) {
442 for(j=0; (j<ny); j++) {
443 for(k=0; k<(nz/2+1); k++) {
444 index = INDEX(i,j,k);
445 if (nTest[index] != 1)
446 fprintf(fp,"Index %d sucks, set %d times\n",index,nTest[index]);
456 void solve_pppm(FILE *fp,t_commrec *cr,
457 t_fftgrid *grid,real ***ghat,rvec box,
458 gmx_bool bVerbose,t_nrnb *nrnb)
463 print_fftgrid(fp,"Q-Real",grid,grid->nxyz,"qreal.pdb",box,TRUE);*/
465 gmxfft3D(grid,GMX_FFT_REAL_TO_COMPLEX,cr);
468 print_fftgrid(fp,"Q-k",grid,1.0,"qk-re.pdb",box,TRUE);
469 print_fftgrid(fp,"Q-k",grid,1.0,"qk-im.pdb",box,FALSE);
470 fprintf(stderr,"Doing convolution\n");
473 convolution(fp,bVerbose,grid,ghat,cr);
476 print_fftgrid(fp,"Convolution",grid,1.0,
477 "convolute.pdb",box,TRUE);*/
479 gmxfft3D(grid,GMX_FFT_COMPLEX_TO_REAL,cr);
482 print_fftgrid(fp,"Potential",grid,1.0,"pot.pdb",box,TRUE);*/
485 npppm = ntot*log((real)ntot)/log(2.0);
486 inc_nrnb(nrnb,eNR_FFT,2*npppm);
487 inc_nrnb(nrnb,eNR_CONV,ntot);
492 static real ***ghat=NULL;
493 static t_fftgrid *grid=NULL;
498 int gmx_pppm_init(FILE *log, t_commrec *cr,
499 const output_env_t oenv, gmx_bool bVerbose,
500 gmx_bool bOld, matrix box,
501 char *ghatfn, t_inputrec *ir,
502 gmx_bool bReproducible)
504 int nx,ny,nz,m,porder;
507 const real tol = 1e-5;
508 rvec box_diag,spacing;
511 gmx_fatal(FARGS,"PPPM is not functional in the current version, we plan to implement PPPM through a small modification of the PME code.");
515 #ifdef GMX_WITHOUT_FFTW
516 gmx_fatal(FARGS,"PPPM used, but GROMACS was compiled without FFTW support!\n");
522 fprintf(log,"Initializing parallel PPPM.\n");
525 "Will use the PPPM algorithm for long-range electrostatics\n");
529 box_diag[m] = box[m][m];
531 if (!gmx_fexist(ghatfn)) {
532 beta[XX]=beta[YY]=beta[ZZ]= 1.85;
538 fprintf(log,"Generating Ghat function\n");
539 fprintf(log,"Grid size is %d x %d x %d\n",nx,ny,nz);
542 if ((nx < 4) || (ny < 4) || (nz < 4))
543 gmx_fatal(FARGS,"Grid must be at least 4 points in all directions");
545 ghat = mk_rgrid(nx,ny,nz);
546 mk_ghat(NULL,nx,ny,nz,ghat,box_diag,
547 ir->rcoulomb_switch,ir->rcoulomb,TRUE,bOld);
550 pr_scalar_gk("generghat.xvg",oenv,nx,ny,nz,box_diag,ghat);
553 fprintf(stderr,"Reading Ghat function from %s\n",ghatfn);
554 ghat = rd_ghat(log,oenv,ghatfn,grids,spacing,beta,&porder,&r1,&rc);
556 /* Check whether cut-offs correspond */
557 if ((fabs(r1-ir->rcoulomb_switch)>tol) || (fabs(rc-ir->rcoulomb)>tol)) {
559 fprintf(log,"rcoulomb_switch = %10.3e rcoulomb = %10.3e"
560 " r1 = %10.3e rc = %10.3e\n",
561 ir->rcoulomb_switch,ir->rcoulomb,r1,rc);
564 gmx_fatal(FARGS,"Cut-off lengths in tpb file and Ghat file %s "
565 "do not match\nCheck your log file!",ghatfn);
568 /* Check whether boxes correspond */
569 for(m=0; (m<DIM); m++)
570 if (fabs(box_diag[m]-grids[m]*spacing[m]) > tol) {
572 pr_rvec(log,0,"box",box_diag,DIM,TRUE);
573 pr_rvec(log,0,"grid-spacing",spacing,DIM,TRUE);
574 pr_ivec(log,0,"grid size",grids,DIM,TRUE);
577 gmx_fatal(FARGS,"Box sizes in tpb file and Ghat file %s do not match\n"
578 "Check your log file!",ghatfn);
582 gmx_fatal(FARGS,"porder = %d, should be 2 in %s",porder,ghatfn);
589 pr_scalar_gk("optimghat.xvg",oenv,nx,ny,nz,box_diag,ghat);
591 /* Now setup the FFT things */
593 cr->mpi_comm_mygroup=cr->mpi_comm_mysim;
595 grid = mk_fftgrid(nx,ny,nz,NULL,NULL,cr,bReproducible);
601 int gmx_pppm_do(FILE *log, gmx_pme_t pme,
604 real charge[], rvec box,
605 real phi[], t_commrec *cr,
608 int pme_order, real *energy)
611 gmx_fatal(FARGS,"PPPM temporarily disabled while working on 2DPME\n");
615 /* Make the grid empty */
618 /* First step: spreading the charges over the grid. */
619 spread_q(log,bVerbose,start,nr,x,charge,box,grid,nrnb);
621 /* In the parallel code we have to sum the grids from neighbouring nodes */
623 gmx_sum_qgrid(pme,cr,grid,GMX_SUM_QGRID_FORWARD);
625 /* Second step: solving the poisson equation in Fourier space */
626 solve_pppm(log,cr,grid,ghat,box,bVerbose,nrnb);
628 /* In the parallel code we have to sum once again... */
630 gmx_sum_qgrid(pme,cr,grid,GMX_SUM_QGRID_BACKWARD);
632 /* Third and last step: gather the forces, energies and potential
635 *energy = gather_f(log,bVerbose,start,nr,x,f,charge,box,
643 static int gmx_pppm_opt_do(FILE *log, gmx_pme_t pme,
644 t_inputrec *ir, gmx_bool bVerbose,
647 real charge[], rvec box,
648 real phi[], t_commrec *cr,
649 t_nrnb *nrnb, rvec beta,
650 t_fftgrid *grid, gmx_bool bOld,
657 fprintf(log,"Generating Ghat function\n");
661 ghat = mk_rgrid(nx,ny,nz);
662 mk_ghat(NULL,nx,ny,nz,ghat,box,ir->rcoulomb_switch,ir->rcoulomb,TRUE,bOld);
664 /* pr_scalar_gk("generghat.xvg",nx,ny,nz,box,ghat); */
666 /* Now start the actual PPPM procedure.
667 * First step: spreading the charges over the grid.
669 /* Make the grid empty */
672 spread_q(log,bVerbose,0,natoms,x,charge,box,grid,nrnb);
674 /* Second step: solving the poisson equation in Fourier space */
675 solve_pppm(log,cr,grid,ghat,box,bVerbose,nrnb);
677 /* Third and last step: gather the forces, energies and potential
680 *energy = gather_f(log,bVerbose,0,natoms,x,f,charge,box,phi,grid,beta,nrnb);
682 free_rgrid(ghat,nx,ny);