diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 9bd0a9f3fc..a484c0b862 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -321,7 +321,7 @@ void PPPM::init() // calculate the final accuracy double estimated_accuracy = final_accuracy(); - + // print stats int ngrid_max,nfft_both_max,nbuf_max; @@ -456,7 +456,7 @@ void PPPM::setup() } } } - + if (differentiation_flag == 1) { compute_gf_ad(); } else { @@ -631,7 +631,7 @@ void PPPM::allocate() memory->create3d_offset(vdz_brick,nzlo_out,nzhi_out,nylo_out,nyhi_out, nxlo_out,nxhi_out,"pppm:vdz_brick"); } - + // summation coeffs memory->create(gf_b,order,"pppm:gf_b"); @@ -751,11 +751,11 @@ void PPPM::deallocate_peratom() if (differentiation_flag != 1) memory->destroy3d_offset(u_brick,nzlo_out,nylo_out,nxlo_out); - + memory->destroy(buf3); memory->destroy(buf4); } - + /* ---------------------------------------------------------------------- set size of FFT grid (nx,ny,nz_pppm) ------------------------------------------------------------------------- */ @@ -792,43 +792,43 @@ void PPPM::set_grid() // reduce it until accuracy target is met if (!gridflag) { - + if (differentiation_flag == 1) { h = h_x = h_y = h_z = 4.0/g_ewald; int count = 0; while (1) { - + // set grid dimension nx_pppm = static_cast (xprd/h_x); ny_pppm = static_cast (yprd/h_y); nz_pppm = static_cast (zprd_slab/h_z); - + if (nx_pppm <= 1) nx_pppm = 2; if (ny_pppm <= 1) ny_pppm = 2; if (nz_pppm <= 1) nz_pppm = 2; - + //set local grid dimension int npey_fft,npez_fft; if (nz_pppm >= nprocs) { npey_fft = 1; npez_fft = nprocs; } else procs2grid2d(nprocs,ny_pppm,nz_pppm,&npey_fft,&npez_fft); - + int me_y = me % npey_fft; int me_z = me / npey_fft; - + nxlo_fft = 0; nxhi_fft = nx_pppm - 1; nylo_fft = me_y*ny_pppm/npey_fft; nyhi_fft = (me_y+1)*ny_pppm/npey_fft - 1; nzlo_fft = me_z*nz_pppm/npez_fft; nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1; - + double df_kspace = compute_df_kspace(); - + count++; - + // break loop if the accuracy has been reached or // too many loops have been performed @@ -1039,20 +1039,20 @@ double PPPM::estimate_ik_error(double h, double prd, bigint natoms) void PPPM::adjust_gewald() { double dx; - + for (int i = 0; i < LARGE; i++) { dx = newton_raphson_f() / derivf(); - g_ewald -= dx; + g_ewald -= dx; if (fabs(newton_raphson_f()) < SMALL) return; } - + char str[128]; sprintf(str, "Could not compute g_ewald"); error->all(FLERR, str); } /* ---------------------------------------------------------------------- - Calculate f(x) using Newton–Raphson solver + Calculate f(x) using Newton-Raphson solver ------------------------------------------------------------------------- */ double PPPM::newton_raphson_f() @@ -1102,12 +1102,12 @@ double PPPM::final_accuracy() double zprd = domain->zprd; double zprd_slab = zprd*slab_volfactor; bigint natoms = atom->natoms; - + double df_kspace = compute_df_kspace(); double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd_slab); double df_rspace = 2.0 * q2_over_sqrt * exp(-g_ewald*g_ewald*cutoff*cutoff); double df_table = estimate_table_accuracy(q2_over_sqrt,df_rspace); - double estimated_accuracy = sqrt(df_kspace*df_kspace + df_rspace*df_rspace + + double estimated_accuracy = sqrt(df_kspace*df_kspace + df_rspace*df_rspace + df_table*df_table); return estimated_accuracy; @@ -1123,7 +1123,7 @@ void PPPM::set_fft_parameters() // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of // global PPPM grid that I own without ghost cells // for slab PPPM, assign z grid as if it were not extended - + nxlo_in = static_cast (comm->xsplit[comm->myloc[0]] * nx_pppm); nxhi_in = static_cast (comm->xsplit[comm->myloc[0]+1] * nx_pppm) - 1; @@ -1133,7 +1133,7 @@ void PPPM::set_fft_parameters() nzlo_in = static_cast (comm->zsplit[comm->myloc[2]] * nz_pppm/slab_volfactor); nzhi_in = static_cast - (comm->zsplit[comm->myloc[2]+1] * nz_pppm/slab_volfactor) - 1; + (comm->zsplit[comm->myloc[2]+1] * nz_pppm/slab_volfactor) - 1; // nlower,nupper = stencil size for mapping particles to PPPM grid @@ -1147,7 +1147,7 @@ void PPPM::set_fft_parameters() else shift = OFFSET; if (order % 2) shiftone = 0.0; else shiftone = 0.5; - + // nlo_out,nhi_out = lower/upper limits of the 3d sub-brick of // global PPPM grid that my particles can contribute charge to // effectively nlo_in,nhi_in + ghost cells @@ -1172,7 +1172,7 @@ void PPPM::set_fft_parameters() sublo = domain->sublo_lamda; subhi = domain->subhi_lamda; } - + double xprd = prd[0]; double yprd = prd[1]; double zprd = prd[2]; @@ -1186,7 +1186,7 @@ void PPPM::set_fft_parameters() dist[1] = cuthalf/domain->prd[1]; dist[2] = cuthalf/domain->prd[2]; } - + int nlo,nhi; nlo = static_cast ((sublo[0]-dist[0]-boxlo[0]) * @@ -1226,7 +1226,7 @@ void PPPM::set_fft_parameters() // that overlay domain I own // proc in that direction tells me via sendrecv() // if no neighbor proc, value is from self since I have ghosts regardless - + int nplanes; MPI_Status status; @@ -1271,7 +1271,7 @@ void PPPM::set_fft_parameters() &nzlo_ghost,1,MPI_INT,comm->procneigh[2][0],0, world,&status); else nzlo_ghost = nplanes; - + // decomposition of FFT mesh // global indices range from 0 to N-1 // proc owns entire x-dimension, clump of columns in y,z dimensions @@ -1299,7 +1299,7 @@ void PPPM::set_fft_parameters() nzhi_fft = (me_z+1)*nz_pppm/npez_fft - 1; // PPPM grid for this proc, including ghosts - + ngrid = (nxhi_out-nxlo_out+1) * (nyhi_out-nylo_out+1) * (nzhi_out-nzlo_out+1); @@ -1393,7 +1393,7 @@ void PPPM::compute_gf_ik() if (triclinic == 0) prd = domain->prd; else prd = domain->prd_lamda; - + double xprd = prd[0]; double yprd = prd[1]; double zprd = prd[2]; @@ -1401,7 +1401,7 @@ void PPPM::compute_gf_ik() double unitkx = (MY_2PI/xprd); double unitky = (MY_2PI/yprd); double unitkz = (MY_2PI/zprd_slab); - + int nbx = static_cast ((g_ewald*xprd/(MY_PI*nx_pppm)) * pow(-log(EPS_HOC),0.25)); int nby = static_cast ((g_ewald*yprd/(MY_PI*ny_pppm)) * @@ -1643,7 +1643,7 @@ void PPPM::compute_sf_precoeff() if (argz != 0.0) wz0[i+2] = pow(sin(argz)/argz,order); argz = 0.5*qz1/nz_pppm; if (argz != 0.0) wz1[i+2] = pow(sin(argz)/argz,order); - argz = 0.5*qz2/nz_pppm; + argz = 0.5*qz2/nz_pppm; if (argz != 0.0) wz2[i+2] = pow(sin(argz)/argz,order); } @@ -1681,8 +1681,6 @@ void PPPM::compute_sf_precoeff() } } - - /* ---------------------------------------------------------------------- ghost-swap to accumulate full density in brick decomposition remap density from 3d brick decomposition to FFT decomposition